1namespace Eigen { 2 3/** \page TopicAliasing Aliasing 4 5In Eigen, aliasing refers to assignment statement in which the same matrix (or array or vector) appears on the 6left and on the right of the assignment operators. Statements like <tt>mat = 2 * mat;</tt> or <tt>mat = 7mat.transpose();</tt> exhibit aliasing. The aliasing in the first example is harmless, but the aliasing in the 8second example leads to unexpected results. This page explains what aliasing is, when it is harmful, and what 9to do about it. 10 11<b>Table of contents</b> 12 - \ref TopicAliasingExamples 13 - \ref TopicAliasingSolution 14 - \ref TopicAliasingCwise 15 - \ref TopicAliasingMatrixMult 16 - \ref TopicAliasingSummary 17 18 19\section TopicAliasingExamples Examples 20 21Here is a simple example exhibiting aliasing: 22 23<table class="example"> 24<tr><th>Example</th><th>Output</th></tr> 25<tr><td> 26\include TopicAliasing_block.cpp 27</td> 28<td> 29\verbinclude TopicAliasing_block.out 30</td></tr></table> 31 32The output is not what one would expect. The problem is the assignment 33\code 34mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2); 35\endcode 36This assignment exhibits aliasing: the coefficient \c mat(1,1) appears both in the block 37<tt>mat.bottomRightCorner(2,2)</tt> on the left-hand side of the assignment and the block 38<tt>mat.topLeftCorner(2,2)</tt> on the right-hand side. After the assignment, the (2,2) entry in the bottom 39right corner should have the value of \c mat(1,1) before the assignment, which is 5. However, the output shows 40that \c mat(2,2) is actually 1. The problem is that Eigen uses lazy evaluation (see 41\ref TopicEigenExpressionTemplates) for <tt>mat.topLeftCorner(2,2)</tt>. The result is similar to 42\code 43mat(1,1) = mat(0,0); 44mat(1,2) = mat(0,1); 45mat(2,1) = mat(1,0); 46mat(2,2) = mat(1,1); 47\endcode 48Thus, \c mat(2,2) is assigned the \e new value of \c mat(1,1) instead of the old value. The next section 49explains how to solve this problem by calling \link DenseBase::eval() eval()\endlink. 50 51Note that if \c mat were a bigger, then the blocks would not overlap, and there would be no aliasing 52problem. This means that in general aliasing cannot be detected at compile time. However, Eigen does detect 53some instances of aliasing, albeit at run time. The following example exhibiting aliasing was mentioned in 54\ref TutorialMatrixArithmetic : 55 56<table class="example"> 57<tr><th>Example</th><th>Output</th></tr> 58<tr><td> 59\include tut_arithmetic_transpose_aliasing.cpp 60</td> 61<td> 62\verbinclude tut_arithmetic_transpose_aliasing.out 63</td></tr></table> 64 65Again, the output shows the aliasing issue. However, by default Eigen uses a run-time assertion to detect this 66and exits with a message like 67 68\verbatim 69void Eigen::DenseBase<Derived>::checkTransposeAliasing(const OtherDerived&) const 70[with OtherDerived = Eigen::Transpose<Eigen::Matrix<int, 2, 2, 0, 2, 2> >, Derived = Eigen::Matrix<int, 2, 2, 0, 2, 2>]: 71Assertion `(!internal::check_transpose_aliasing_selector<Scalar,internal::blas_traits<Derived>::IsTransposed,OtherDerived>::run(internal::extract_data(derived()), other)) 72&& "aliasing detected during tranposition, use transposeInPlace() or evaluate the rhs into a temporary using .eval()"' failed. 73\endverbatim 74 75The user can turn Eigen's run-time assertions like the one to detect this aliasing problem off by defining the 76EIGEN_NO_DEBUG macro, and the above program was compiled with this macro turned off in order to illustrate the 77aliasing problem. See \ref TopicAssertions for more information about Eigen's run-time assertions. 78 79 80\section TopicAliasingSolution Resolving aliasing issues 81 82If you understand the cause of the aliasing issue, then it is obvious what must happen to solve it: Eigen has 83to evaluate the right-hand side fully into a temporary matrix/array and then assign it to the left-hand 84side. The function \link DenseBase::eval() eval() \endlink does precisely that. 85 86For example, here is the corrected version of the first example above: 87 88<table class="example"> 89<tr><th>Example</th><th>Output</th></tr> 90<tr><td> 91\include TopicAliasing_block_correct.cpp 92</td> 93<td> 94\verbinclude TopicAliasing_block_correct.out 95</td></tr></table> 96 97Now, \c mat(2,2) equals 5 after the assignment, as it should be. 98 99The same solution also works for the second example, with the transpose: simply replace the line 100<tt>a = a.transpose();</tt> with <tt>a = a.transpose().eval();</tt>. However, in this common case there is a 101better solution. Eigen provides the special-purpose function 102\link DenseBase::transposeInPlace() transposeInPlace() \endlink which replaces a matrix by its transpose. 103This is shown below: 104 105<table class="example"> 106<tr><th>Example</th><th>Output</th></tr> 107<tr><td> 108\include tut_arithmetic_transpose_inplace.cpp 109</td> 110<td> 111\verbinclude tut_arithmetic_transpose_inplace.out 112</td></tr></table> 113 114If an xxxInPlace() function is available, then it is best to use it, because it indicates more clearly what you 115are doing. This may also allow Eigen to optimize more aggressively. These are some of the xxxInPlace() 116functions provided: 117 118<table class="manual"> 119<tr><th>Original function</th><th>In-place function</th></tr> 120<tr> <td> MatrixBase::adjoint() </td> <td> MatrixBase::adjointInPlace() </td> </tr> 121<tr class="alt"> <td> DenseBase::reverse() </td> <td> DenseBase::reverseInPlace() </td> </tr> 122<tr> <td> LDLT::solve() </td> <td> LDLT::solveInPlace() </td> </tr> 123<tr class="alt"> <td> LLT::solve() </td> <td> LLT::solveInPlace() </td> </tr> 124<tr> <td> TriangularView::solve() </td> <td> TriangularView::solveInPlace() </td> </tr> 125<tr class="alt"> <td> DenseBase::transpose() </td> <td> DenseBase::transposeInPlace() </td> </tr> 126</table> 127 128 129\section TopicAliasingCwise Aliasing and component-wise operations 130 131As explained above, it may be dangerous if the same matrix or array occurs on both the left-hand side and the 132right-hand side of an assignment operator, and it is then often necessary to evaluate the right-hand side 133explicitly. However, applying component-wise operations (such as matrix addition, scalar multiplication and 134array multiplication) is safe. 135 136The following example has only component-wise operations. Thus, there is no need for .eval() even though 137the same matrix appears on both sides of the assignments. 138 139<table class="example"> 140<tr><th>Example</th><th>Output</th></tr> 141<tr><td> 142\include TopicAliasing_cwise.cpp 143</td> 144<td> 145\verbinclude TopicAliasing_cwise.out 146</td></tr></table> 147 148In general, an assignment is safe if the (i,j) entry of the expression on the right-hand side depends only on 149the (i,j) entry of the matrix or array on the left-hand side and not on any other entries. In that case it is 150not necessary to evaluate the right-hand side explicitly. 151 152 153\section TopicAliasingMatrixMult Aliasing and matrix multiplication 154 155Matrix multiplication is the only operation in Eigen that assumes aliasing by default. Thus, if \c matA is a 156matrix, then the statement <tt>matA = matA * matA;</tt> is safe. All other operations in Eigen assume that 157there are no aliasing problems, either because the result is assigned to a different matrix or because it is a 158component-wise operation. 159 160<table class="example"> 161<tr><th>Example</th><th>Output</th></tr> 162<tr><td> 163\include TopicAliasing_mult1.cpp 164</td> 165<td> 166\verbinclude TopicAliasing_mult1.out 167</td></tr></table> 168 169However, this comes at a price. When executing the expression <tt>matA = matA * matA</tt>, Eigen evaluates the 170product in a temporary matrix which is assigned to \c matA after the computation. This is fine. But Eigen does 171the same when the product is assigned to a different matrix (e.g., <tt>matB = matA * matA</tt>). In that case, 172it is more efficient to evaluate the product directly into \c matB instead of evaluating it first into a 173temporary matrix and copying that matrix to \c matB. 174 175The user can indicate with the \link MatrixBase::noalias() noalias()\endlink function that there is no 176aliasing, as follows: <tt>matB.noalias() = matA * matA</tt>. This allows Eigen to evaluate the matrix product 177<tt>matA * matA</tt> directly into \c matB. 178 179<table class="example"> 180<tr><th>Example</th><th>Output</th></tr> 181<tr><td> 182\include TopicAliasing_mult2.cpp 183</td> 184<td> 185\verbinclude TopicAliasing_mult2.out 186</td></tr></table> 187 188Of course, you should not use \c noalias() when there is in fact aliasing taking place. If you do, then you 189may get wrong results: 190 191<table class="example"> 192<tr><th>Example</th><th>Output</th></tr> 193<tr><td> 194\include TopicAliasing_mult3.cpp 195</td> 196<td> 197\verbinclude TopicAliasing_mult3.out 198</td></tr></table> 199 200 201\section TopicAliasingSummary Summary 202 203Aliasing occurs when the same matrix or array coefficients appear both on the left- and the right-hand side of 204an assignment operator. 205 - Aliasing is harmless with coefficient-wise computations; this includes scalar multiplication and matrix or 206 array addition. 207 - When you multiply two matrices, Eigen assumes that aliasing occurs. If you know that there is no aliasing, 208 then you can use \link MatrixBase::noalias() noalias()\endlink. 209 - In all other situations, Eigen assumes that there is no aliasing issue and thus gives the wrong result if 210 aliasing does in fact occur. To prevent this, you have to use \link DenseBase::eval() eval() \endlink or 211 one of the xxxInPlace() functions. 212 213*/ 214} 215