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