Eigen的速度為什麽這麽快？

2015-03-06財經

歡迎來到C++數值最佳化的世界！Eigen能寫的快和它的設計思路有關，涵蓋了演算法加速的幾個方法。Eigen的設計思路，是把所有能最佳化的步驟放在 編譯時 去最佳化。要想執行時變快，寫Eigen程式碼的工程師需要顯式地告訴Eigen矩陣的特性，告訴的越多，越能提供空間讓編譯器在編譯時加速演算法

比如一個矩陣乘法

#include <Eigen/Dense> Eigen :: MatrixXd Jacobian_i = Eigen :: MatrixXd :: Random ( 10 , 10 ); Eigen :: MatrixXd Jacobian_j = Eigen :: MatrixXd :: Random ( 10 , 10 ); Eigen :: MatrixXd Hessian = Eigen :: MatrixXd :: Zero ( 10 , 10 ); Hessian += Jacobian_i . transpose () * Jacobian_j ;

這個是一個稠密矩陣相乘。首先需要知道的是Eigen並不是一步一步地先做轉置，再去乘，而是使用lazy evaluation的方法。 Lazy evalution ，是把計算本身盡可能放在最後做，減少記憶體存取。

在程式碼上來看，就是 transpose()函式並沒有做轉置，和operator*() 本身也沒有做乘法，只是更新了一下flag，真正執行計算的是 operator+=()。從Eigen/src/Core/EigenBase.h 可以看到，是在operator+=() 這邊才真正去做記憶體讀取和計算。

第二個可以用的最佳化技巧是改變記憶體分配的方式，使用Eigen時應該盡可能 用靜態記憶體代替動態記憶體 。在上述矩陣計算中，三個矩陣都是動態記憶體分配，因為Eigen::MatrixXd 是如下的縮寫：

typedef MatrixXd Matrix < double , Dynamic , Dynamic , ColMajor >

MatrixBase第二和第三個選項是行列的長度，有一項是Dynamic就會用動態記憶體分配。所以在已知矩陣大小時應盡可能聲明大小，比如Matrix<double, 10, 10>。如果記憶體在整個程式中大小會變，但知道最大可能的大小，都可以告知Eigen，Eigen同樣會選擇用靜態記憶體。

#include <Eigen/Dense> Eigen :: Matrix < double , Eigen :: Dynamic , Eigen :: Dynamic , Eigen :: ColMajor , 10 , 10 > Jacobian_i ; Eigen :: Matrix < double , Eigen :: Dynamic , Eigen :: Dynamic , Eigen :: ColMajor , 10 , 10 > Jacobian_j ; Eigen :: Matrix < double , Eigen :: Dynamic , Eigen :: Dynamic , Eigen :: ColMajor , 10 , 10 > Hessian = Eigen :: Matrix < double , Eigen :: Dynamic , Eigen :: Dynamic , Eigen :: ColMajor , 10 , 10 >:: Zero (); Hessian += Jacobian_i . transpose () * Jacobian_j ;

靜態記憶體分配不但讓我們節省了new/delete的開銷，還給Eigen內部繼續最佳化提供了可能。Eigen內建Single-Instruction-Multiple-Data( SIMD )指令集，對稠密矩陣有很好的最佳化，如果能觸發CPU SIMD的指令，能收獲成倍的計算效率。

最後是記憶體的讀寫。Lazy evaluation的另一個特點是不生成中間變量，減少記憶體搬運次數，而Eigen為了防止矩陣覆蓋自己，對矩陣-矩陣乘法會生成一個中間變量。如果我們知道等式左右兩邊沒有相同的項，則可以通知Eigen去取消中間變量。

Hessian . noalias () += Jacobian_i . transpose () * Jacobian_j ;

以上是普通稠密矩陣計算。如果矩陣本身有自身的性質，都可以通知Eigen，讓Eigen用對應的加速方式。比如正定矩陣可以只用上三角進行計算，並且在求解時使用Eigen::LLT這樣又快又數值穩定的解法，詳情見Linear algebra and decompositions。

總之Eigen的高效既來自於演算法最佳化技巧，又來自於工程師本身對公式和記憶體計算的理解。

參考資料：

Eigen原始碼：Files · master · libeigen / eigen

Eigen官方教程：The Matrix class