文章目录

混淆example 1 解决混淆问题混淆和component级的操作混淆和矩阵的乘法总结

---

博客转自：fengbingchun Eigen/OpenBLAS

混淆

在Eigen中，当变量同时出现在左值和右值，赋值操作可能会带来混淆问题，此篇博客解释什么是混淆，什么时候是有害的，怎么使用做。

example 1

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

在 mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2); 赋值中展示了混淆。

mat(1,1) 在bottomRightCorner(2,2)和topLeftCorner(2,2)都存在。赋值结果中mat(2,2)本应该赋予操作前mat(1,1)的值=5。但是，最终程序结果mat(2,2)=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);

下面会解释如何通过eval()来解决这个问题。

混淆还会在缩小矩阵时出现，比如 vec = vec.head(n) 和 mat = mat.block(i,j,r,c)。

一般来说，混淆在编译阶段很难被检测到。比如第一个例子，如果mat再大一些可能就不会出现混淆了。但是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_NO_DEBUG宏，在编译时关闭运行时的断言。

解决混淆问题

Eigen需要把右值赋值为一个临时matrix/array，然后再将临时值赋值给左值，便可以解决混淆。eval()函数实现了这个功能。

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.transpose().eval(); ，当然我们最好使用transposeInPlace()。如果存在xxxInPlace函数，推荐使用这类函数，它们更加清晰地标明了你在做什么。提供的这类函数：

OriginIn-placeMatrixBase::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()。

混淆和component级的操作

组件级是指整体的操作，比如matrix加法、scalar乘、array乘等，这类操作是安全的，不会出现混淆。

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;

输出

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

混淆和矩阵的乘法

在Eigen中，矩阵的乘法一般都会出现混淆。除非是方阵（实质是元素级的乘）。

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

其他的操作，Eigen默认都是存在混淆的。所以Eigen对矩阵乘法自动引入了临时变量，对的matA=matA*matA这是必须的，但是对matB=matA*matA这样便是不必要的了。我们可以使用noalias()函数来声明这里没有混淆，matA*matA的结果可以直接赋值为matB。

matB.noalias() = matA * matA;

从Eigen3.3开始，如果目标矩阵resize且结果不直接赋值给目标矩阵，默认不存在混淆。

MatrixXf A(2,2), B(3,2); B << 2, 0, 0, 3, 1, 1; A << 2, 0, 0, -2; A = (B * A).cwiseAbs();//cwiseAbs（）不直接赋给目标 //A = (B * A).eval().cwiseAbs() cout << A;

当然，对于任何混淆问题，都可以通过matA=(matB*matA).eval()来解决。

---

总结

当相同的矩阵或array在等式左右都出现时，很容易出现混淆。

compnent级别的操作不用考虑混淆。矩阵相乘，Eigen默认会解决混淆问题，如果你确定不会出现混淆，可以使用noalias（）来提效。混淆出现时，可以用eval()和xxxInPlace()函数解决。