123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237 |
- 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 <tt>mat = 2 * mat;</tt> or <tt>mat =
- mat.transpose();</tt> 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:
- <table class="example">
- <tr><th>Example</th><th>Output</th></tr>
- <tr><td>
- \include TopicAliasing_block.cpp
- </td>
- <td>
- \verbinclude TopicAliasing_block.out
- </td></tr></table>
- 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
- <tt>mat.bottomRightCorner(2,2)</tt> on the left-hand side of the assignment and the block
- <tt>mat.topLeftCorner(2,2)</tt> 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 <tt>mat.topLeftCorner(2,2)</tt>. 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 <tt>vec =
- vec.head(n)</tt> and <tt>mat = mat.block(i,j,r,c)</tt> 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 :
- <table class="example">
- <tr><th>Example</th><th>Output</th></tr>
- <tr><td>
- \include tut_arithmetic_transpose_aliasing.cpp
- </td>
- <td>
- \verbinclude tut_arithmetic_transpose_aliasing.out
- </td></tr></table>
- 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<Derived>::checkTransposeAliasing(const OtherDerived&) const
- [with OtherDerived = Eigen::Transpose<Eigen::Matrix<int, 2, 2, 0, 2, 2> >, Derived = Eigen::Matrix<int, 2, 2, 0, 2, 2>]:
- Assertion `(!internal::check_transpose_aliasing_selector<Scalar,internal::blas_traits<Derived>::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:
- <table class="example">
- <tr><th>Example</th><th>Output</th></tr>
- <tr><td>
- \include TopicAliasing_block_correct.cpp
- </td>
- <td>
- \verbinclude TopicAliasing_block_correct.out
- </td></tr></table>
- 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
- <tt>a = a.transpose();</tt> with <tt>a = a.transpose().eval();</tt>. 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:
- <table class="example">
- <tr><th>Example</th><th>Output</th></tr>
- <tr><td>
- \include tut_arithmetic_transpose_inplace.cpp
- </td>
- <td>
- \verbinclude tut_arithmetic_transpose_inplace.out
- </td></tr></table>
- 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:
- <table class="manual">
- <tr><th>Original function</th><th>In-place function</th></tr>
- <tr> <td> MatrixBase::adjoint() </td> <td> MatrixBase::adjointInPlace() </td> </tr>
- <tr class="alt"> <td> DenseBase::reverse() </td> <td> DenseBase::reverseInPlace() </td> </tr>
- <tr> <td> LDLT::solve() </td> <td> LDLT::solveInPlace() </td> </tr>
- <tr class="alt"> <td> LLT::solve() </td> <td> LLT::solveInPlace() </td> </tr>
- <tr> <td> TriangularView::solve() </td> <td> TriangularView::solveInPlace() </td> </tr>
- <tr class="alt"> <td> DenseBase::transpose() </td> <td> DenseBase::transposeInPlace() </td> </tr>
- </table>
- In the special case where a matrix or vector is shrunk using an expression like <tt>vec = vec.head(n)</tt>,
- 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.
- <table class="example">
- <tr><th>Example</th><th>Output</th></tr>
- <tr><td>
- \include TopicAliasing_cwise.cpp
- </td>
- <td>
- \verbinclude TopicAliasing_cwise.out
- </td></tr></table>
- 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, <strong>under the
- condition that the destination matrix is not resized</strong>.
- Thus, if \c matA is a \b squared matrix, then the statement <tt>matA = matA * matA;</tt> 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.
- <table class="example">
- <tr><th>Example</th><th>Output</th></tr>
- <tr><td>
- \include TopicAliasing_mult1.cpp
- </td>
- <td>
- \verbinclude TopicAliasing_mult1.out
- </td></tr></table>
- However, this comes at a price. When executing the expression <tt>matA = matA * matA</tt>, %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., <tt>matB = matA * matA</tt>). 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: <tt>matB.noalias() = matA * matA</tt>. This allows %Eigen to evaluate the matrix product
- <tt>matA * matA</tt> directly into \c matB.
- <table class="example">
- <tr><th>Example</th><th>Output</th></tr>
- <tr><td>
- \include TopicAliasing_mult2.cpp
- </td>
- <td>
- \verbinclude TopicAliasing_mult2.out
- </td></tr></table>
- 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:
- <table class="example">
- <tr><th>Example</th><th>Output</th></tr>
- <tr><td>
- \include TopicAliasing_mult3.cpp
- </td>
- <td>
- \verbinclude TopicAliasing_mult3.out
- </td></tr></table>
- 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:
- <table class="example">
- <tr><th>Example</th><th>Output</th></tr>
- <tr><td>
- \include TopicAliasing_mult4.cpp
- </td>
- <td>
- \verbinclude TopicAliasing_mult4.out
- </td></tr></table>
- As for any aliasing issue, you can resolve it by explicitly evaluating the expression prior to assignment:
- <table class="example">
- <tr><th>Example</th><th>Output</th></tr>
- <tr><td>
- \include TopicAliasing_mult5.cpp
- </td>
- <td>
- \verbinclude TopicAliasing_mult5.out
- </td></tr></table>
- \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.
- */
- }
|