簡介

別名混亂Aliasing是指在賦值表達式中，一個Eigen對象(矩陣、數組、向量)同時出現在左值和右值表達式中，比如v = v*2; m = m.transpose();;

別名混亂會引起錯誤，從而產生問題，比如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 problemmat.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 7After 'mat = 2 * mat', mat = 2 48 14After the subtraction, it becomes1 48 13After squaring, it becomes1 1664 169Doing everything at once yields1 1664 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 44 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()函數。

總結

以上是生活随笔為你收集整理的Eigen入门之密集矩阵 9 - 别名混乱Aliasing的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。