Eigen 3.2稀疏矩陣入門
https://my.oschina.net/cvnote/blog/166980
Eigen自帶的稀疏矩陣分解功能包括LDLt、LLt分解（即Cholesky分解，這個功能是LGPL許可，不是Eigen的MPL許可）、LU分解、QR分解（這個是3.2版本之後正式Release的）、共軛梯度解矩陣等。另外還提供了到第三方稀疏矩陣庫的C++介面，包括著名的SuiteSparse系列（這個系列非常強大，有機會要好好研究一下）的SparseQR、UmfPack等。(歡迎訪問計算機視覺研究筆記cvnote.info和關注新浪微博@cvnote )
基本稀疏矩陣操作
Eigen中使用 Eigen::Triplet<Scalar>來記錄一個非零元素的行、列、值，填充一個稀疏矩陣，只需要將所有表示非零元素的Triplet放在一個 std::vector中即可傳入即可。除了求逆等功能外，Eigen::SparseMatrix 有和 Eigen::Matrix幾乎一樣的各種成員操作函式，並且可以方便混用。
比如這樣：
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
#include <iostream> #include "eigen/Eigen/Eigen"
int main ( ) { // 矩陣： // 0 6.1 0 0 // 0 0 7.2 0 // 0 0 8.3 0 // 0 0 0 0 using namespace Eigen ; SparseMatrix < double > A ( 4 , 4 ) ; std :: vector < Triplet < double > > triplets ;
int r [ 3 ] = { 0 , 1 , 2 } ; // 非零元素的行號 int c [ 3 ] = { 1 , 2 , 2 } ; // 非零元素的列號 double val [ 3 ] = { 6.1 , 7.2 , 8.3 } ; // 非零元素的值 for ( int i = 0 ; i < 3 ; ++ i ) triplets . emplace_back ( r [ i ] , c [ i ] , val [ i ] ) ; // 填充Triplet
A . setFromTriplets ( triplets . begin ( ) , triplets . end ( ) ) ; // 初始化係數矩陣 std :: cout << "A = " << A << std :: endl ;
MatrixXd B = A ; // 可以和普通稠密矩陣方便轉換 std :: cout << "B = \n" << B << std :: endl ; std :: cout << "A * B = \n" << A * B << std :: endl ; // 可以各種運算 std :: cout << "A * A = \n" << A * A << std :: endl ; return 0 ; } |
用Eigen進行矩陣分解
首先複習一下Cholesky（LLt）、QR和LU分解，說的不對的地方歡迎數學大牛和數值計算大牛來指教和補充。一般來講LLt分解可以理解成給矩陣開平方，類比於開平方一般針對正數而言，LLt分解則限定矩陣需為正定的 Hermitian矩陣（自共軛矩陣，即對稱的實數矩陣或對稱元素共軛的複數矩陣）。LU分解則稍微放鬆一點，可用於一般的方陣（順便提一句LU分解是圖靈發明的）。QR則可用於一般矩陣，結果也是最穩定的。分解演算法的效率，三者都是O(n^3)的，具體系數三者大概是Cholesky:LU:QR=1:2:4。Google可以找到很多相關資料，比如我看了這個。
下面試一下用Eigen自帶的Eigen::SparseQR 進行我最喜歡的QR分解（其實我更喜歡SVD）。
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 |
#include <iostream> #include "eigen/Eigen/Eigen" #include "eigen/Eigen/SparseQR"
int main ( ) { using namespace Eigen ; SparseMatrix < double > A ( 4 , 4 ) ; std :: vector < Triplet < double > > triplets ;
// 初始化非零元素 int r [ 3 ] = { 0 , 1 , 2 } ; int c [ 3 ] = { 1 , 2 , 2 } ; double val [ 3 ] = { 6.1 , 7.2 , 8.3 } ; for ( int i = 0 ; i < 3 ; ++ i ) triplets . emplace_back ( r [ i ] , c [ i ] , val [ i ] ) ;
// 初始化稀疏矩陣 A . setFromTriplets ( triplets . begin ( ) , triplets . end ( ) ) ; std :: cout << "A = \n" << A << std :: endl ;
// 一個QR分解的例項 SparseQR < SparseMatrix < double > , AMDOrdering < int > > qr ; // 計算分解 qr . compute ( A ) ; // 求一個A x = b Vector4d b ( 1 , 2 , 3 , 4 ) ; Vector4d x = qr . solve ( b ) ; std :: cout << "x = \n" << x ; std :: cout << "A x = \n" << A * x ;
return 0 ; } |
這樣就用QR分解解了一個係數矩陣，A和上面的例子是一樣的，注意到這個Ax=b其實是沒有確定解的，A(1:3, 2:3)是over determined的，剩下的部分又是非滿秩under determined的，這個QR分解對於A(1:3, 2:3)給了最小二乘解，其他位返回了0。
另一個注意的地方就是 SparseQR<SparseMatrix<double>, AMDOrdering<int> >的第二個模板引數，是一個矩陣重排列（ordering）的方法，為什麼要重排列呢，wikipedia的LU分解詞條給了一個例子可以大概解釋一下，某些矩陣沒有重排直接分解可能會失敗。Eigen提供了三種重排列方法，參見OrderingMethods Module。關於矩陣重排列的細節求數值計算牛人指點！我一般就隨便選一個填進去了>_<。
除了解方程，這個QR例項也可以用下面程式碼返回Q和R矩陣：
| 1 2 3 |
SparseMatrix < double > Q , R ; qr . matrixQ ( ) . evalTo ( Q ) ; R = qr . matrixR ( ) ; |
注意到Q和R的返回方法不一樣，猜測是因為 matrixQ() 成員好像是沒有完整儲存Q矩陣（元素太多？）。
用 Eigen::SPQR 呼叫SuiteSparseQR進行QR分解
SuiteSparseQR效率很高，但是C風格介面比較不好用，Eigen提供了 Eigen::SPQR 的介面封裝比如和上面同樣的程式可以這樣寫：
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 |
#include <iostream> #include "eigen/Eigen/Eigen" #include "eigen/Eigen/SPQRSupport"
int main ( ) { using namespace Eigen ; SparseMatrix < double > A ( 4 , 4 ) ; std :: vector < Triplet < double > > triplets ;
// 初始化非零元素 int r [ 3 ] = { 0 , 1 , 2 } ; int c [ 3 ] = { 1 , 2 , 2 } ; double val [ 3 ] = { 6.1 , 7.2 , 8.3 } ; for ( int i = 0 ; i < 3 ; ++ i ) triplets . emplace_back ( r [ i ] , c [ i ] , val [ i ] ) ;
// 初始化稀疏矩陣 A . setFromTriplets ( triplets . begin ( ) , triplets . end ( ) ) ; std :: cout << "A = \n" << A << std :: endl ;
// 一個QR分解的例項 SPQR < SparseMatrix < double > > qr ; // 計算分解 qr . compute ( A ) ; // 求一個A x = b Vector4d b ( 1 , 2 , 3 , 4 ) ; Vector4d x = qr . solve ( b ) ; std :: cout << "x = \n" << x ; std :: cout << "A x = \n" << A * x ;
return 0 ; } |
如果你有用過SuiteSparseQR的話，會覺得這個介面真心好用多了。編譯這個程式除了spqr庫還需要連結blas庫、lapack庫、cholmod庫（SuiteSparse的另一元件），有一點麻煩。比如我在ubuntu，使用 apt-get install libsuitesparse-* 安裝了suitesparse標頭檔案到/usr/include/suitesparse 目錄，使用如下命令編譯。
| 1 |
g ++ spqr . cpp - std = c ++ 11 - I / usr / include / suitesparse - lcholmod - lspqr - llapack - lblas |
注意lapack要在blas前面，spqr要在lapack前面。用了c++11是因為上面程式碼偷懶用了emplace_back ，和矩陣庫沒關係。
一些效率的經驗
SuiteSparseQR畢竟實現更好一些，我的一些經驗是比自帶Eigen::SparseQR快50%左右吧。
相關文章
歡迎訪問計算機視覺研究筆記cvnote.info和關注新浪微博@cvnote