lesson18 Aliasing

Aliasing

Dense matrix and array manipulation

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.

Examples

Here is a simple example exhibiting aliasing:

Example	Output
MatrixXi mat(3,3); mat << 1, 2, 3, 4, 5, 6, 7, 8, 9; cout << "Here is the matrix mat:\n" << mat << endl; // This assignment shows the aliasing problem mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2); cout << "After the assignment, mat = \n" << mat << endl;	Here is the matrix mat: 1 2 3 4 5 6 7 8 9 After the assignment, mat = 1 2 3 4 1 2 7 4 1

The output is not what one would expect. The problem is the assignment

mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);

This assignment exhibits aliasing: the coefficient 

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 

mat(1,1)

 before the assignment, which is 5. However, the output shows that 

mat(2,2)

 is actually 1. The problem is that
Eigen uses lazy evaluation (see Expression templates in Eigen)
for 

mat.topLeftCorner(2,2)

. The result is similar to

mat(1,1) = mat(0,0);

mat(1,2) = mat(0,1);

mat(2,1) = mat(1,0);

mat(2,2) = mat(1,1);

Thus, 

mat(2,2)

 is assigned the new value of 

mat(1,1)

 instead of the old value. The next section explains how to solve this problem by calling eval().
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 

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 Matrix
and vector arithmetic :

Example	Output
Matrix2i a; a << 1, 2, 3, 4; cout << "Here is the matrix a:\n" << a << endl; a = a.transpose(); // !!! do NOT do this !!! cout << "and the result of the aliasing effect:\n" << a << endl;	Here is the matrix a: 1 2 3 4 and the result of the aliasing effect: 1 2 2 4

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

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.

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 Assertions for more information about Eigen's run-time
assertions.

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 functioneval() does precisely
that.
For example, here is the corrected version of the first example above:

Example	Output
MatrixXi mat(3,3); mat << 1, 2, 3, 4, 5, 6, 7, 8, 9; cout << "Here is the matrix mat:\n" << mat << endl; // The eval() solves the aliasing problem mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2).eval(); cout << "After the assignment, mat = \n" << mat << endl;	Here is the matrix mat: 1 2 3 4 5 6 7 8 9 After the assignment, mat = 1 2 3 4 1 2 7 4 5

Now, 

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 transposeInPlace() which
replaces a matrix by its transpose. This is shown below:

Example	Output
MatrixXf a(2,3); a << 1, 2, 3, 4, 5, 6; cout << "Here is the initial matrix a:\n" << a << endl; a.transposeInPlace(); cout << "and after being transposed:\n" << a << endl;	Here is the initial matrix a: 1 2 3 4 5 6 and after being transposed: 1 4 2 5 3 6

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 function	In-place function
MatrixBase::adjoint()	MatrixBase::adjointInPlace()
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 useconservativeResize() .

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.

Example

Output

MatrixXf mat(2,2); mat << 1, 2, 4, 7; cout << "Here is the matrix mat:\n" << mat << endl << endl; mat = 2 * mat; cout << "After 'mat = 2 * mat', mat = \n" << mat << endl << endl; mat = mat - MatrixXf::Identity(2,2); cout << "After the subtraction, it becomes\n" << mat << endl << endl; ArrayXXf arr = mat; arr = arr.square(); cout << "After squaring, it becomes\n" << arr << endl << endl; // Combining all operations in one statement: mat << 1, 2, 4, 7; mat = (2 * mat - MatrixXf::Identity(2,2)).array().square(); cout << "Doing everything at once yields\n" << mat << endl << endl;

Here is the matrix mat: 1 2 4 7 After 'mat = 2 * mat', mat = 2 4 8 14 After the subtraction, it becomes 1 4 8 13 After squaring, it becomes 1 16 64 169 Doing everything at once yields 1 16 64 169

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.

Aliasing and matrix multiplication

Matrix multiplication
is the only operation in Eigen that assumes aliasing by default. Thus, if 

matA

 is a 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.

Example	Output
MatrixXf matA(2,2); matA << 2, 0, 0, 2; matA = matA * matA; cout << matA;	4 0 0 4

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 

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 

matB

 instead of evaluating it first into a temporary
matrix and copying that matrix to 

matB

.
The user can indicate with the noalias() function
that there is no aliasing, as follows: 

matB.noalias() = matA * matA

. This allows Eigen to evaluate the matrix product 

matA * matA

 directly into 

matB

.

Example	Output
MatrixXf matA(2,2), matB(2,2); matA << 2, 0, 0, 2; // Simple but not quite as efficient matB = matA * matA; cout << matB << endl << endl; // More complicated but also more efficient matB.noalias() = matA * matA; cout << matB;	4 0 0 4 4 0 0 4

Of course, you should not use 

noalias()

 when there is in fact aliasing taking place. If you do, then you may get wrong results:

Example	Output
MatrixXf matA(2,2); matA << 2, 0, 0, 2; matA.noalias() = matA * matA; cout << matA;	0 0 0 0

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 noalias().
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 eval() or
one of the xxxInPlace() functions.