Eigen学习笔记(10)-混淆

原文:Eigen官网-Aliasing

在Eigen中,当变量同时出现在左值和右值,赋值操作可能会带来混淆问题。比如:mat = 2 * matmat = mat.transpose() ,对于第一个表达式,是无害的,但是第二个表达式则会造成意想不到的后果。这一篇将解释什么是混淆,什么时候是有害的,怎么使用做。

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

出现上述混淆的原因是:Eigen使用了lazy evaluation

一般来说,混淆在编译阶段很难被检测到。比如第一个例子,如果mat再大一些可能就不会出现混淆了。但是Eigen可以在运行时检测某些混淆,如前面讲的例子。

2. 解决混淆问题

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 = a.transpose(), 可以使用a = a.transpose().eval()进行替换,另外也可以使用transposeInPlace()来替换transpose()函数,以解决混淆问题。

如果存在xxxInPlace()函数,推荐使用这类函数,它们更加清晰地标明了你在做什么。如下是提供的这类函数:

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

而针对vec = vec.head(n)这种情况,推荐使用conservativeResize()

3. 混淆和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;
// 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

4. 混淆和矩阵的乘法

在假设目标矩阵的尺寸不变的情况下,矩阵乘法是Eigen中唯一默认会出现混淆的操作。因此,对于方阵,matA = matA * matA是安全的。因为Eigen默认进行了混淆问题的解决。

对于Eigen中的其他所有操作,Eigen默认是不存在混淆问题的,或者是因为目标矩阵变成了具有不同尺寸大小的矩阵,或者是元素级的操作。因此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

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

当然,对于任何混淆问题,都可以通过matA=(matB*matA).eval()在赋值之前显示求表达式的值来解决。

5. 总结

当同一个matrix或array在等号左右都出现时,很容易出现混淆。

  • compnent级别的操作不存在混淆问题,比如标量乘法、matrix加法、array加法。
  • 矩阵相乘,Eigen默认会解决混淆问题,如果你确定不会出现混淆,可以使用noalias()来提高效率。
  • 在其他情况中,Eigen假设不存在混淆问题,因此当混淆出现时,Eigen会给出错误的结果,但是可以用eval()或者xxxInPlace()函数解决。

参考:

“Eigen教程(10)”

© 版权声明
THE END
喜欢就支持一下吧
点赞510 分享
评论 抢沙发
头像
欢迎您留下宝贵的见解!
提交
头像

昵称

取消
昵称表情代码图片

    暂无评论内容