blob: a8f1644284e82af3cdc653548a4e74259bbc81a4 [file] [log] [blame]
 namespace Eigen { /** \eigenManualPage TopicAliasing Aliasing In %Eigen, aliasing refers to assignment statement in which the same matrix (or array or vector) appears on the left and on the right of the assignment operators. Statements like mat = 2 * mat; or mat = mat.transpose(); exhibit aliasing. The aliasing in the first example is harmless, but the aliasing in the second example leads to unexpected results. This page explains what aliasing is, when it is harmful, and what to do about it. \eigenAutoToc \section TopicAliasingExamples Examples Here is a simple example exhibiting aliasing:
ExampleOutput
\include TopicAliasing_block.cpp \verbinclude TopicAliasing_block.out
The output is not what one would expect. The problem is the assignment \code mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2); \endcode This assignment exhibits aliasing: the coefficient \c mat(1,1) appears both in the block mat.bottomRightCorner(2,2) on the left-hand side of the assignment and the block mat.topLeftCorner(2,2) on the right-hand side. After the assignment, the (2,2) entry in the bottom right corner should have the value of \c mat(1,1) before the assignment, which is 5. However, the output shows that \c mat(2,2) is actually 1. The problem is that %Eigen uses lazy evaluation (see \ref TopicEigenExpressionTemplates) for mat.topLeftCorner(2,2). The result is similar to \code mat(1,1) = mat(0,0); mat(1,2) = mat(0,1); mat(2,1) = mat(1,0); mat(2,2) = mat(1,1); \endcode Thus, \c mat(2,2) is assigned the \e new value of \c mat(1,1) instead of the old value. The next section explains how to solve this problem by calling \link DenseBase::eval() eval()\endlink. Aliasing occurs more naturally when trying to shrink a matrix. For example, the expressions vec = vec.head(n) and mat = mat.block(i,j,r,c) exhibit aliasing. In general, aliasing cannot be detected at compile time: if \c mat in the first example were a bit bigger, then the blocks would not overlap, and there would be no aliasing problem. However, %Eigen does detect some instances of aliasing, albeit at run time. The following example exhibiting aliasing was mentioned in \ref TutorialMatrixArithmetic :
ExampleOutput
\include tut_arithmetic_transpose_aliasing.cpp \verbinclude tut_arithmetic_transpose_aliasing.out
Again, the output shows the aliasing issue. However, by default %Eigen uses a run-time assertion to detect this and exits with a message like \verbatim void Eigen::DenseBase::checkTransposeAliasing(const OtherDerived&) const [with OtherDerived = Eigen::Transpose >, Derived = Eigen::Matrix]: Assertion `(!internal::check_transpose_aliasing_selector::IsTransposed,OtherDerived>::run(internal::extract_data(derived()), other)) && "aliasing detected during transposition, use transposeInPlace() or evaluate the rhs into a temporary using .eval()"' failed. \endverbatim The user can turn %Eigen's run-time assertions like the one to detect this aliasing problem off by defining the EIGEN_NO_DEBUG macro, and the above program was compiled with this macro turned off in order to illustrate the aliasing problem. See \ref TopicAssertions for more information about %Eigen's run-time assertions. \section TopicAliasingSolution Resolving aliasing issues If you understand the cause of the aliasing issue, then it is obvious what must happen to solve it: %Eigen has to evaluate the right-hand side fully into a temporary matrix/array and then assign it to the left-hand side. The function \link DenseBase::eval() eval() \endlink does precisely that. For example, here is the corrected version of the first example above:
ExampleOutput
\include TopicAliasing_block_correct.cpp \verbinclude TopicAliasing_block_correct.out
Now, \c mat(2,2) equals 5 after the assignment, as it should be. The same solution also works for the second example, with the transpose: simply replace the line a = a.transpose(); with a = a.transpose().eval();. However, in this common case there is a better solution. %Eigen provides the special-purpose function \link DenseBase::transposeInPlace() transposeInPlace() \endlink which replaces a matrix by its transpose. This is shown below:
ExampleOutput
\include tut_arithmetic_transpose_inplace.cpp \verbinclude tut_arithmetic_transpose_inplace.out
If an xxxInPlace() function is available, then it is best to use it, because it indicates more clearly what you are doing. This may also allow %Eigen to optimize more aggressively. These are some of the xxxInPlace() functions provided:
Original functionIn-place function
DenseBase::reverse() DenseBase::reverseInPlace()
LDLT::solve() LDLT::solveInPlace()
LLT::solve() LLT::solveInPlace()
TriangularView::solve() TriangularView::solveInPlace()
DenseBase::transpose() DenseBase::transposeInPlace()
In the special case where a matrix or vector is shrunk using an expression like vec = vec.head(n), you can use \link PlainObjectBase::conservativeResize() conservativeResize() \endlink. \section TopicAliasingCwise Aliasing and component-wise operations As explained above, it may be dangerous if the same matrix or array occurs on both the left-hand side and the right-hand side of an assignment operator, and it is then often necessary to evaluate the right-hand side explicitly. However, applying component-wise operations (such as matrix addition, scalar multiplication and array multiplication) is safe. The following example has only component-wise operations. Thus, there is no need for \link DenseBase::eval() eval() \endlink even though the same matrix appears on both sides of the assignments.
ExampleOutput
\include TopicAliasing_cwise.cpp \verbinclude TopicAliasing_cwise.out
In general, an assignment is safe if the (i,j) entry of the expression on the right-hand side depends only on the (i,j) entry of the matrix or array on the left-hand side and not on any other entries. In that case it is not necessary to evaluate the right-hand side explicitly. \section TopicAliasingMatrixMult Aliasing and matrix multiplication Matrix multiplication is the only operation in %Eigen that assumes aliasing by default, under the condition that the destination matrix is not resized. Thus, if \c matA is a \b squared matrix, then the statement matA = matA * matA; is safe. All other operations in %Eigen assume that there are no aliasing problems, either because the result is assigned to a different matrix or because it is a component-wise operation.
ExampleOutput
\include TopicAliasing_mult1.cpp \verbinclude TopicAliasing_mult1.out
However, this comes at a price. When executing the expression matA = matA * matA, %Eigen evaluates the product in a temporary matrix which is assigned to \c matA after the computation. This is fine. But %Eigen does the same when the product is assigned to a different matrix (e.g., matB = matA * matA). In that case, it is more efficient to evaluate the product directly into \c matB instead of evaluating it first into a temporary matrix and copying that matrix to \c matB. The user can indicate with the \link MatrixBase::noalias() noalias()\endlink function that there is no aliasing, as follows: matB.noalias() = matA * matA. This allows %Eigen to evaluate the matrix product matA * matA directly into \c matB.
ExampleOutput
\include TopicAliasing_mult2.cpp \verbinclude TopicAliasing_mult2.out
Of course, you should not use \c noalias() when there is in fact aliasing taking place. If you do, then you may get wrong results:
ExampleOutput
\include TopicAliasing_mult3.cpp \verbinclude TopicAliasing_mult3.out
Moreover, 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. Therefore, the following example is also wrong:
ExampleOutput
\include TopicAliasing_mult4.cpp \verbinclude TopicAliasing_mult4.out
As for any aliasing issue, you can resolve it by explicitly evaluating the expression prior to assignment:
ExampleOutput
\include TopicAliasing_mult5.cpp \verbinclude TopicAliasing_mult5.out
\section TopicAliasingSummary Summary Aliasing occurs when the same matrix or array coefficients appear both on the left- and the right-hand side of an assignment operator. - Aliasing is harmless with coefficient-wise computations; this includes scalar multiplication and matrix or array addition. - When you multiply two matrices, %Eigen assumes that aliasing occurs. If you know that there is no aliasing, then you can use \link MatrixBase::noalias() noalias()\endlink. - In all other situations, %Eigen assumes that there is no aliasing issue and thus gives the wrong result if aliasing does in fact occur. To prevent this, you have to use \link DenseBase::eval() eval() \endlink or one of the xxxInPlace() functions. */ }