| namespace Eigen { |
| |
| /** \page 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. |
| |
| <b>Table of contents</b> |
| - \ref TopicAliasingExamples |
| - \ref TopicAliasingSolution |
| - \ref TopicAliasingCwise |
| - \ref TopicAliasingMatrixMult |
| - \ref TopicAliasingSummary |
| |
| |
| \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. |
| |
| Note that if \c mat were a bigger, then the blocks would not overlap, and there would be no aliasing |
| problem. This means that in general aliasing cannot be detected at compile time. 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 tranposition, 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> |
| |
| |
| \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 .eval() 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. Thus, if \c matA is a |
| 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> |
| |
| |
| \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. |
| |
| */ |
| } |