混叠问题(Aliasing)

混叠问题指的是在进行赋值操作的时候，赋值表达式的两边存在重叠的矩阵区域，例如 mat = 2*mat; 或 mat = mat.transpose(); 。前一个表达式中的混叠问题不会对结果产生影响，但后一个表达式中的混叠问题会导致不可预料的结果。本节将解释什么是混叠问题，何时有害，如何处理。

举例

MatrixXi mat(3,3);
mat << 1,2,3,4,5,6,7,8,9;
cout << "mat:\n" << mat << endl;

mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);
cout << "After the assignment:\n" << mat << endl;

输出：

mat:
1 2 3
4 5 6
7 8 9
After the assignment: 
1 2 3
4 1 2
7 4 1

可见，输出并不像我们所期望的一样将左上角的 $2\times2$ 子矩阵赋值给右下角的 $2\times2$ 子矩阵，这种情况正是由混叠问题引起的，赋值表达式的两端存在重叠的矩阵区域，即 $mat(1,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(1,1) 已经不再是旧的值，而已经被赋上了新的值。

尝试进行矩阵压缩的时候很容易出现混叠问题，例如 vec = vec.head(n); 或 mat = mat.block(i,j,row,col); 等。

通常情况下，混叠问题在编译过程中是不会被识别出来的，但运行时会以异常信息 (run-time assertion) 做出提示，可以使用 EIGEN_NO_DEBUG 退出debug模式关闭异常提示，但最好不要。

如何解决混叠问题

显然，最简单的解决方法就是把右侧的计算结果赋给一个临时变量，计算结束后再将这个临时变量赋值给左侧的变量，Eigen中使用 eval() 函数实现。用 eval() 纠正上例中出现的问题：

MatrixXi mat(3,3);
mat << 1,2,3,4,5,6,7,8,9;
cout << "mat:\n" << mat << endl;

mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2).eval();
cout << "After the assignment:\n" << mat << endl;

输出：

mat:
1 2 3
4 5 6
7 8 9
After the assignment: 
1 2 3
4 1 2
7 4 5

又例如对矩阵作转置时，直接使用 a = a.transpose(); 会导致混叠问题：

Matrix2i a; a << 1, 2, 3, 4;
cout << "Matrix a:\n" << a << endl;
 
a = a.transpose(); // !!! do NOT do this !!!
cout << "Transposed a:\n" << a << endl;

输出：

Matrix a:
1 2
3 4
Transposed a:
1 2
2 4

与上面相同的解决方法，即使用 a = a.transpose().eval() 代替。此外，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

类似的原地处理函数还有：

MatrixBase::adjoint() <-> MatrixBase::adjointInPlace()
DenseBase::reverse() <-> DenseBase::reverseInPlace()
LDLT::solve() <-> LDLT::solveInPlace()
LLT::solve() <-> LLT::solveInPlace()
TriangularView::solve() <-> TriangularView::solveInPlace()
DenseBase::transpose() <-> DenseBase::transposeInPlace()

针对矩阵/向量压缩处理（取几行几列/取前几个元素），例如 vec = vec.head(n) ，可以使用 conservativeResize() 实现原地处理。

混叠问题与元素级操作

在元素级操作中，尽管赋值表达式两边也有可能出现同一个矩阵（例如矩阵加法，数组乘法，标量乘法等），但不会对结果产生影响，即元素级操作中的混叠是安全的，因此也就不需要 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

总而言之，当表达式右侧的 (i,j) 元素仅仅依赖于表达式左侧对应的 (i,j) 元素时，赋值操作是安全的，也就不再需要多余的 eval() 操作。

混叠问题与矩阵乘法

在目标矩阵未调整大小的情况下，矩阵乘法是Eigen中唯一默认进行混叠处理的操作。换言之，若矩阵 matA 为平方型矩阵（行数 = 列数），则赋值语句 matA = matA*matA 是安全的。

除此之外的所有操作，Eigen都默认不存在由混叠导致的错误，要么是因为结果被赋值给了一个不同的矩阵，要么是因为当前进行的是元素级操作。

不过Eigen对矩阵乘法进行的这种处理是需要付出一定计算资源的代价的。在进行上述乘法的过程中，Eigen首先将乘积计算入一个临时矩阵，结束计算后再将这个临时矩阵赋值给 matA ，在当前情况下这是合理的。但相同的操作也会对 matB = matA * matA 进行，这就会导致一定程度上的计算资源浪费。在这种情况下，Eigen提供了 noalias() 函数进行直接赋值而不用经过中间的临时矩阵：matB.noalias() = matA * matA 。这种操作一定要慎用！！！

总结

当表达式两端同时出现相同的矩阵/数组时会出现混叠

元素级操作时无害，包括标量/数组乘法和矩阵/数组加法
两个行列相等的矩阵相乘，Eigen默认进行混叠处理；若明确不会出现混叠，使用 noalias()
其他情况下Eigen默认不进行混叠处理，因此容易在某些情况下产生错误答案。使用 eval() 或 xxxInPlace() 函数来避免此类情况