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