1. 程式人生 > >Eigen中的noalias(): 解決矩陣運算的混淆問題

Eigen中的noalias(): 解決矩陣運算的混淆問題

需要 右值 什麽 原因 lan sin 一個 eba ner

作者:@houkai
本文為作者原創,轉載請註明出處:http://www.cnblogs.com/houkai/p/6349990.html


目錄

混淆
例子
解決混淆問題
混淆和component級的操作。
混淆和矩陣的乘法
總結

整理下Eigen庫的教程,參考:http://eigen.tuxfamily.org/dox/index.html

混淆

在Eigen中,當變量同時出現在左值和右值,賦值操作可能會帶來混淆問題。這一篇將解釋什麽是混淆,什麽時候是有害的,怎麽使用做。

例子

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 = a.transpose().eval(); ,當然我們最好使用 transposeInPlace()。如果存在xxxInPlace函數,推薦使用這類函數,它們更加清晰地標明了你在做什麽。提供的這類函數:

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

而針對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在等式左右都出現時,很容易出現混淆。

  1. compnent級別的操作不用考慮混淆。
  2. 矩陣相乘,Eigen默認會解決混淆問題,如果你確定不會出現混淆,可以使用noalias()來提效。
  3. 混淆出現時,可以用eval()和xxxInPlace()函數解決。

Eigen中的noalias(): 解決矩陣運算的混淆問題