reference resources: https://blog.csdn.net/whereismatrix/article/details/104569080

brief introduction

Aliasing means that in the assignment expression, an Eigen object (matrix, array, vector) appears in both lvalue and rvalue expressions, such as v = v*2; m = m.transpose();；

Alias confusion can cause errors and problems;

Here we will introduce what alias confusion is and how to avoid errors.

Aligning example

First give an example to get an intuitive impression from the example. This example expects the block of 2X2 in the upper left corner to cover the block in the lower right corner.

//matrix_aliasing1.cpp
#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;


int main()
{
    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;
}

Execution results:

$ g++   -I /usr/local/include/eigen3 matrix_aliasing1.cpp -o matrix_aliasing1
promote:eigen david$ 
promote:eigen david$ ./matrix_aliasing1 
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

From the results, we can see that the actual results are not as good as we expected. Obviously, mat bottomRightCorner(2,2) = mat. topLeftCorner(2,2); There's something wrong with your writing.

The coefficient mat(1,1) exists in both bottomrightcorner (2,2) and topleftcorner (2,2). After assignment, the coefficient (2,2) of the block in the lower right corner will be saved as the value 5 of the original mat(1,1), but its actual value is 1.

The problem is that Eigen uses lazy evaluation delay evaluation. Therefore, the operation process is similar:

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

Therefore, the new value assigned to mat(2,2) finally comes from mat (0,0), mat (1,1) < -- mat (0,0), so the result is not as expected.

When shrinking a matrix, it is easier for code to generate alias confusion. For example, the expression VEC = VEC head(n)； mat = mat.block(i,j,r,c)； All produce errors.

The trouble is that aliasing is confusing. At compile time, the compiler does not recognize the processing. However, Eigen can detect some alias confusion at execution time.

The following example:

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;

If compiling and running, the execution effect should be as follows (but more operations are required, see the following):

Here is the matrix a:
1 2
3 4
and the result of the aliasing effect:
1 2
2 4

Of course, we know that there is a problem of alias confusion, and the result is true. However, when Eigen executes, it uses runtime assertions and gives a promotion message. So:

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.

Tip: to turn off such assertions, define a macro: EIGEN_NO_DEBUG and compile the program. Then you can compile and run to get such an error result.

Resolve alias confusion

Only knowing the root cause of alias confusion can solve the problem. The method is to calculate and evaluate the right value, assign a temporary variable, and then assign a value to the left value to solve the problem.

For the evaluation of the right value, use the function eval() complete.

Eigen provides examples:

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;

Execution results:

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

A = A. transfer () in the above problem example; It will cause problems. Simply modify it to a = a. transfer() eval(); You can be right.

But a more general and better way is to use the special function provided by Eigen to operate - Eigen provides the transportinplace () function to complete this task.

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;

Execution results:

Here is the initial matrix a:
1 2 3
4 5 6
and after being transposed:
1 4
2 5
3 6

If only one operation has the corresponding Eigen function xxinplace (), it can be used safely, which is more efficient and correct.

The following list lists the corresponding functions provided by Eigen:

Primitive 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()

Other situations:
In a = a.head(n) mentioned above, you can use the function conserveresize () to perform the operation.

Alias confusion and component oriented operation

As mentioned above, if the matrix or array object is on the left and right sides of the assignment expression at the same time, an error may occur. But it can be improved by explicitly calling the evaluation function eval() on the right.

However, it is safe to use component-oriented operations, such as matrix addition, scalar multiplication, or array multiplication.

The following example has only component oriented operations, which are secure, and all explicit calls that do not require eval() evaluation.

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;

Execution results:

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

Alias confusion and matrix multiplication

By default, in Eigen, matrix multiplication is the only operation that assumes aliasing confusion when the target matrix is not resized. So, if matA is a square matrix, then matA = matA * matA; Computing is safe. All other operations are considered safe without alias confusion, either because the calculation result is given a different matrix object, or component-oriented operation functions are used.

Here is a simple example:

MatrixXf matA(2,2); 
matA << 2, 0,  0, 2;

matA = matA * matA;

cout << matA;

Execution results:

4 0
0 4

However, this also brings a certain price. For example: matA = matA * matA; Operation, the matrix multiplication expression here will be given a temporary variable object, and then the calculation result will be given to the matA object. The assignment expression matB = matA * matA; Just as many temporary variables will be used. In this case, a more efficient way is to calculate the product directly and assign it to the matB object instead of using a temporary variable.

This requires the use of the noalias () function, which tells you that there is no alias confusion. Use this way: matb noalias() = matA * matA;.

Example:

MatrixXf matA(2,2), matB(2,2); 
matA << 2, 0,  0, 2;

// This is a calculation assignment mode using temporary variable turnover.
matB = matA * matA;
cout << matB << endl << endl;

// This is a more efficient model.
matB.noalias() = matA * matA;
cout << matB;

Execution results:

4 0
0 4

4 0
0 4

Note: from eigen3 3 at the beginning, if the size of the target matrix will change during multiplication calculation, it is also considered that there will be no alias confusion error. At this time, you need to deal with it yourself.

MatrixXf A(2,2), B(3,2);
B << 2, 0,  0, 3, 1, 1;
A << 2, 0, 0, -2;

A = (B * A).cwiseAbs();

cout << A;

Here, the execution (3X2) matrix is multiplied by the (2X2) matrix to obtain the (3X2) matrix. Matrix A starts from (2X2) – > (3X2):

4 0
0 6
2 2

But there are still mistakes.

The correct implementation method should be as follows:

MatrixXf A(2,2), B(3,2);
B << 2, 0,  0, 3, 1, 1;
A << 2, 0, 0, -2;

A = (B * A).eval().cwiseAbs();
cout << A;

Note: here a = (b * a) eval(). cwiseAbs();， Eval() was explicitly called.

summary

Alias confusion occurs when matrix or array objects are located on both left and right sides during assignment operation, and the same multiple coefficients are operated.

• However, alias confusion in coefficient oriented calculation is harmless, including scalar multiplication of matrix or array and addition of matrix or array.
• When two matrices are multiplied, Eigen assumes alias confusion. If we are sure that there is no alias confusion error, we can use noalias().
• In other cases, Eigen assumes that there is no alias confusion, but if the real situation is that there is an alias confusion, an error will occur. To prevent this, you must call the eval() or xxinplace() function.