【C++】Eigen入门之密集矩阵 8 - 别名混乱Aliasing

参考：https://blog.csdn.net/whereismatrix/article/details/104569080

简介

别名混乱Aliasing是指在赋值表达式中，一个Eigen对象(矩阵、数组、向量)同时出现在左值和右值表达式中，比如v = v*2; m = m.transpose();；

别名混乱会引起错误，从而产生问题；

这里将介绍什么是别名混乱，如何来避免发生错误的情况。

aliasing 示例

先给一个示例，从例子中来获得直观印象。本示例期望将左上角的2X2的块覆盖右下角的块。

//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;
}

执行结果：

$ 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

从结果可以看到，实际的结果并不如我们锁期望。很明显，mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);的写法有问题。

系数mat(1,1)在bottomRightCorner(2,2) , topLeftCorner(2,2)都存在，在赋值后，右下角的块的系数(2,2)将保存为原mat(1,1)的值5，但实际其值为1。

问题是 Eigen使用的是lazy evaluation延迟求值。所以，操作的过程类似：

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

于是，mat(2,2) 被赋予的新值最终来自于mat(0,0), mat(1,1) <-- mat(0,0)，所以结果不是预期的样子。

在收缩(shrink)一个矩阵时，更容易有代码产生别名混乱。就比如：表达式语句 vec = vec.head(n)； mat = mat.block(i,j,r,c)；都产生错误。

比较麻烦的是，别名混乱 aliasing 在编译时，编译器并不能识别处理。然而，Eigen可以在执行时侦测到某些别名混乱。

下面的示例：

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

当然，我们知道了这有别名混乱的问题，结果也确实如此。 但Eigen在执行时，使用了运行时断言，给出了提升消息。如此这样：

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.

提示： 要关闭这样的断言，请定义一个宏：EIGEN_NO_DEBUG，并编译程序。然后你就可以编译运行，得到这样的错误结果。

解决别名混乱

只有知道了别名混乱的根源，那么就可以解决问题了。方法就是，对右值就行计算求值，并赋予临时变量，再为左值赋值，就可以解决问题了。

对右值的求值，使用函数.eval()完成。

Eigen提供了示例：

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

上面问题示例中的a = a.transpose();会产生问题，简单地修改为a = a.transpose().eval();就可以正确了。

但更通用，更好的方式是：使用Eigen提供的专用函数来操作 – Eigen提供了transposeInPlace()函数可以完成这个任务。

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

只有一个操作有对应的Eigen函数xxxInPlace()，那就放心地使用，其更高效，更正确。

下面列表列出了Eigen提供的对应函数：

原函数	In-place函数
MatrixBase::adjoint()	MatrixBase::adjointInPlace()
DenseBase::reverse()	DenseBase::reverseInPlace()
LDLT::solve()	LDLT::solveInPlace()
LLT::solve()	LLT::solveInPlace()
TriangularView::solve()	TriangularView::solveInPlace()
DenseBase::transpose()	DenseBase::transposeInPlace()

其他的情况：
在上面提到的a = a.head(n)，可以使用函数conservativeResize()来执行操作。

别名混乱与面向组件的操作

如上所述，如果矩阵或者数组对象同时在赋值表达式的左右两边，这时就可能发生错误。但可以通过在右边显式调用求值函数eval()来改进。

然而，采用面向组件的操作，比如矩阵相加，乘以标量，或者数组乘法，这些都是安全的。

下面示例只有面向组件的操作，它们是安全的，所有无需eval()求值的显式调用。

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

别名混乱与矩阵乘法

缺省状态下，在Eigen内，在目标矩阵不调整尺寸大小时，矩阵乘法是唯一一个假定有别名混乱的操作。所以，如果matA是一个方阵，那么matA = matA * matA;计算就是安全的。所有其它的操作都认为是安全的，没有别名混乱的问题，要么是因为计算结果赋予了一个不同的矩阵对象，要么使用了面向组件的操作函数。

如下有简单示例：

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

matA = matA * matA;

cout << matA;

执行结果：

4 0
0 4

然而，这也带来了一定的代价。比如：matA = matA * matA;运算，这里的矩阵乘法表达式，会赋予一个临时变量对象，然后计算结果再赋予matA对象。而赋值表达式matB = matA * matA;也会一样多使用一个临时变量。这种情况下，更有效率的方式是直接计算乘积，然后赋值给matB对象，而不多用一个临时变量。

这就需要使用noalias()函数来完成工作，其告知这里没有别名混乱。使用方式如此：matB.noalias() = matA * matA;。

示例:

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

// 这是使用了临时变量周转的计算赋值模式。
matB = matA * matA;
cout << matB << endl << endl;

// 这是比较高效的模式。
matB.noalias() = matA * matA;
cout << matB;

执行结果：

4 0
0 4

4 0
0 4

注意：从Eigen3.3开始，如果乘法计算时，目标矩阵的尺寸大小会发生变化时，也认为不会有别名混乱的错误。这时，你需要自己处理。

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;

这里执行（3X2）矩阵被（2X2）矩阵乘，执行结果得到（3X2）矩阵。矩阵A从（2X2）–> （3X2）：

4 0
0 6
2 2

但这里依然有错。

正确的实现方法应该如下：

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;

注意：这里的A = (B * A).eval().cwiseAbs();，显式调用了eval()。

总结

别名混乱发生于矩阵或数组对象在赋值操作时同时位于左右两边，对同样的多个系数进行操作了。

• 但在面向系数计算的别名混乱是无害的，这包括：标量乘于矩阵或数组，矩阵或数组的加法。
• 而在2个矩阵相乘时，Eigen假定会发生别名混乱。如果我们确定这里并没有发生别名混乱的错误，我们可以使用noalias()。
• 剩下的其他情况下，Eigen假定没有别名混乱的问题，但如果，此时真实情况是有别名混乱的问题发生的，则会导致错误发生。为防止此情况发生，必须调用eval()或者xxxInPlace()函数。