[136] | 1 | namespace Eigen {
|
---|
| 2 |
|
---|
| 3 | /** \eigenManualPage TopicAliasing Aliasing
|
---|
| 4 |
|
---|
| 5 | In %Eigen, aliasing refers to assignment statement in which the same matrix (or array or vector) appears on the
|
---|
| 6 | left and on the right of the assignment operators. Statements like <tt>mat = 2 * mat;</tt> or <tt>mat =
|
---|
| 7 | mat.transpose();</tt> exhibit aliasing. The aliasing in the first example is harmless, but the aliasing in the
|
---|
| 8 | second example leads to unexpected results. This page explains what aliasing is, when it is harmful, and what
|
---|
| 9 | to do about it.
|
---|
| 10 |
|
---|
| 11 | \eigenAutoToc
|
---|
| 12 |
|
---|
| 13 |
|
---|
| 14 | \section TopicAliasingExamples Examples
|
---|
| 15 |
|
---|
| 16 | Here 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 |
|
---|
| 27 | The output is not what one would expect. The problem is the assignment
|
---|
| 28 | \code
|
---|
| 29 | mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);
|
---|
| 30 | \endcode
|
---|
| 31 | This 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
|
---|
| 34 | right corner should have the value of \c mat(1,1) before the assignment, which is 5. However, the output shows
|
---|
| 35 | that \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
|
---|
| 38 | mat(1,1) = mat(0,0);
|
---|
| 39 | mat(1,2) = mat(0,1);
|
---|
| 40 | mat(2,1) = mat(1,0);
|
---|
| 41 | mat(2,2) = mat(1,1);
|
---|
| 42 | \endcode
|
---|
| 43 | Thus, \c mat(2,2) is assigned the \e new value of \c mat(1,1) instead of the old value. The next section
|
---|
| 44 | explains how to solve this problem by calling \link DenseBase::eval() eval()\endlink.
|
---|
| 45 |
|
---|
| 46 | Aliasing occurs more naturally when trying to shrink a matrix. For example, the expressions <tt>vec =
|
---|
| 47 | vec.head(n)</tt> and <tt>mat = mat.block(i,j,r,c)</tt> exhibit aliasing.
|
---|
| 48 |
|
---|
| 49 | In general, aliasing cannot be detected at compile time: if \c mat in the first example were a bit bigger,
|
---|
| 50 | then the blocks would not overlap, and there would be no aliasing problem. However, %Eigen does detect some
|
---|
| 51 | instances of aliasing, albeit at run time. The following example exhibiting aliasing was mentioned in \ref
|
---|
| 52 | TutorialMatrixArithmetic :
|
---|
| 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 |
|
---|
| 63 | Again, the output shows the aliasing issue. However, by default %Eigen uses a run-time assertion to detect this
|
---|
| 64 | and exits with a message like
|
---|
| 65 |
|
---|
| 66 | \verbatim
|
---|
| 67 | void 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>]:
|
---|
| 69 | Assertion `(!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 |
|
---|
| 73 | The user can turn %Eigen's run-time assertions like the one to detect this aliasing problem off by defining the
|
---|
| 74 | EIGEN_NO_DEBUG macro, and the above program was compiled with this macro turned off in order to illustrate the
|
---|
| 75 | aliasing problem. See \ref TopicAssertions for more information about %Eigen's run-time assertions.
|
---|
| 76 |
|
---|
| 77 |
|
---|
| 78 | \section TopicAliasingSolution Resolving aliasing issues
|
---|
| 79 |
|
---|
| 80 | If you understand the cause of the aliasing issue, then it is obvious what must happen to solve it: %Eigen has
|
---|
| 81 | to evaluate the right-hand side fully into a temporary matrix/array and then assign it to the left-hand
|
---|
| 82 | side. The function \link DenseBase::eval() eval() \endlink does precisely that.
|
---|
| 83 |
|
---|
| 84 | For 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 |
|
---|
| 95 | Now, \c mat(2,2) equals 5 after the assignment, as it should be.
|
---|
| 96 |
|
---|
| 97 | The 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
|
---|
| 99 | better solution. %Eigen provides the special-purpose function
|
---|
| 100 | \link DenseBase::transposeInPlace() transposeInPlace() \endlink which replaces a matrix by its transpose.
|
---|
| 101 | This 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 |
|
---|
| 112 | If an xxxInPlace() function is available, then it is best to use it, because it indicates more clearly what you
|
---|
| 113 | are doing. This may also allow %Eigen to optimize more aggressively. These are some of the xxxInPlace()
|
---|
| 114 | functions 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 |
|
---|
| 126 | In the special case where a matrix or vector is shrunk using an expression like <tt>vec = vec.head(n)</tt>,
|
---|
| 127 | you can use \link PlainObjectBase::conservativeResize() conservativeResize() \endlink.
|
---|
| 128 |
|
---|
| 129 |
|
---|
| 130 | \section TopicAliasingCwise Aliasing and component-wise operations
|
---|
| 131 |
|
---|
| 132 | As explained above, it may be dangerous if the same matrix or array occurs on both the left-hand side and the
|
---|
| 133 | right-hand side of an assignment operator, and it is then often necessary to evaluate the right-hand side
|
---|
| 134 | explicitly. However, applying component-wise operations (such as matrix addition, scalar multiplication and
|
---|
| 135 | array multiplication) is safe.
|
---|
| 136 |
|
---|
| 137 | The following example has only component-wise operations. Thus, there is no need for \link DenseBase::eval()
|
---|
| 138 | eval() \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 |
|
---|
| 149 | In general, an assignment is safe if the (i,j) entry of the expression on the right-hand side depends only on
|
---|
| 150 | 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
|
---|
| 151 | not necessary to evaluate the right-hand side explicitly.
|
---|
| 152 |
|
---|
| 153 |
|
---|
| 154 | \section TopicAliasingMatrixMult Aliasing and matrix multiplication
|
---|
| 155 |
|
---|
| 156 | Matrix multiplication is the only operation in %Eigen that assumes aliasing by default. Thus, if \c matA is a
|
---|
| 157 | matrix, then the statement <tt>matA = matA * matA;</tt> is safe. All other operations in %Eigen assume that
|
---|
| 158 | there are no aliasing problems, either because the result is assigned to a different matrix or because it is a
|
---|
| 159 | component-wise operation.
|
---|
| 160 |
|
---|
| 161 | <table class="example">
|
---|
| 162 | <tr><th>Example</th><th>Output</th></tr>
|
---|
| 163 | <tr><td>
|
---|
| 164 | \include TopicAliasing_mult1.cpp
|
---|
| 165 | </td>
|
---|
| 166 | <td>
|
---|
| 167 | \verbinclude TopicAliasing_mult1.out
|
---|
| 168 | </td></tr></table>
|
---|
| 169 |
|
---|
| 170 | However, this comes at a price. When executing the expression <tt>matA = matA * matA</tt>, %Eigen evaluates the
|
---|
| 171 | product in a temporary matrix which is assigned to \c matA after the computation. This is fine. But %Eigen does
|
---|
| 172 | the same when the product is assigned to a different matrix (e.g., <tt>matB = matA * matA</tt>). In that case,
|
---|
| 173 | it is more efficient to evaluate the product directly into \c matB instead of evaluating it first into a
|
---|
| 174 | temporary matrix and copying that matrix to \c matB.
|
---|
| 175 |
|
---|
| 176 | The user can indicate with the \link MatrixBase::noalias() noalias()\endlink function that there is no
|
---|
| 177 | aliasing, as follows: <tt>matB.noalias() = matA * matA</tt>. This allows %Eigen to evaluate the matrix product
|
---|
| 178 | <tt>matA * matA</tt> directly into \c matB.
|
---|
| 179 |
|
---|
| 180 | <table class="example">
|
---|
| 181 | <tr><th>Example</th><th>Output</th></tr>
|
---|
| 182 | <tr><td>
|
---|
| 183 | \include TopicAliasing_mult2.cpp
|
---|
| 184 | </td>
|
---|
| 185 | <td>
|
---|
| 186 | \verbinclude TopicAliasing_mult2.out
|
---|
| 187 | </td></tr></table>
|
---|
| 188 |
|
---|
| 189 | Of course, you should not use \c noalias() when there is in fact aliasing taking place. If you do, then you
|
---|
| 190 | may get wrong results:
|
---|
| 191 |
|
---|
| 192 | <table class="example">
|
---|
| 193 | <tr><th>Example</th><th>Output</th></tr>
|
---|
| 194 | <tr><td>
|
---|
| 195 | \include TopicAliasing_mult3.cpp
|
---|
| 196 | </td>
|
---|
| 197 | <td>
|
---|
| 198 | \verbinclude TopicAliasing_mult3.out
|
---|
| 199 | </td></tr></table>
|
---|
| 200 |
|
---|
| 201 |
|
---|
| 202 | \section TopicAliasingSummary Summary
|
---|
| 203 |
|
---|
| 204 | Aliasing occurs when the same matrix or array coefficients appear both on the left- and the right-hand side of
|
---|
| 205 | an assignment operator.
|
---|
| 206 | - Aliasing is harmless with coefficient-wise computations; this includes scalar multiplication and matrix or
|
---|
| 207 | array addition.
|
---|
| 208 | - When you multiply two matrices, %Eigen assumes that aliasing occurs. If you know that there is no aliasing,
|
---|
| 209 | then you can use \link MatrixBase::noalias() noalias()\endlink.
|
---|
| 210 | - In all other situations, %Eigen assumes that there is no aliasing issue and thus gives the wrong result if
|
---|
| 211 | aliasing does in fact occur. To prevent this, you have to use \link DenseBase::eval() eval() \endlink or
|
---|
| 212 | one of the xxxInPlace() functions.
|
---|
| 213 |
|
---|
| 214 | */
|
---|
| 215 | }
|
---|