中国在线观看免费_ Eigen 3.2稀疏矩陣入門(mén)

轉自 http://my.oschina.net/cvnote/blog/166980

Eigen自帶的稀疏矩陣分解功能包括LDLt、LLt分解（即Cholesky分解，這個(gè)功能是LGPL許可，不是Eigen的MPL許可）、LU分解、QR分解（這個(gè)是3.2版本之后正式Release的）、共軛梯度解矩陣等。另外還提供了到第三方稀疏矩陣庫的C++接口，包括著(zhù)名的SuiteSparse系列（這個(gè)系列非常強大，有機會(huì )要好好研究一下）的SparseQR、UmfPack等。(歡迎訪(fǎng)問(wèn)計算機視覺(jué)研究筆記cvnote.info和關(guān)注新浪微博 @cvnote )

基本稀疏矩陣操作

Eigen中使用 Eigen::Triplet<Scalar>來(lái)記錄一個(gè)非零元素的行、列、值，填充一個(gè)稀疏矩陣，只需要將所有表示非零元素的Triplet放在一個(gè) std::vector中即可傳入即可。除了求逆等功能外，Eigen::SparseMatrix 有和 Eigen::Matrix幾乎一樣的各種成員操作函數，并且可以方便混用。

比如這樣：

 #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進(jìn)行矩陣分解

首先復習一下Cholesky（LLt）、QR和LU分解，說(shuō)的不對的地方歡迎數學(xué)大牛和數值計算大牛來(lái)指教和補充。一般來(lái)講LLt分解可以理解成給矩陣開(kāi)平方，類(lèi)比于開(kāi)平方一般針對正數而言，LLt分解則限定矩陣需為正定的 Hermitian矩陣（自共軛矩陣，即對稱(chēng)的實(shí)數矩陣或對稱(chēng)元素共軛的復數矩陣）。LU分解則稍微放松一點(diǎn)，可用于一般的方陣（順便提一句LU分解是圖靈發(fā)明的）。QR則可用于一般矩陣，結果也是最穩定的。分解算法的效率，三者都是O(n^3)的，具體系數三者大概是Cholesky:LU:QR=1:2:4。Google可以找到很多相關(guān)資料，比如我看了這個(gè)。

下面試一下用Eigen自帶的Eigen::SparseQR 進(jìn)行我最喜歡的QR分解（其實(shí)我更喜歡SVD）。

 #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 ; 
      // 一個(gè)QR分解的實(shí)例 
      SparseQR < SparseMatrix < double > ,   AMDOrdering < int >   >   qr ; 
      // 計算分解 
      qr . compute ( A ) ; 
      // 求一個(gè)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分解解了一個(gè)系數矩陣，A和上面的例子是一樣的，注意到這個(gè)Ax=b其實(shí)是沒(méi)有確定解的，A(1:3, 2:3)是over determined的，剩下的部分又是非滿(mǎn)秩under determined的，這個(gè)QR分解對于A(yíng)(1:3, 2:3)給了最小二乘解，其他位返回了0。

另一個(gè)注意的地方就是 SparseQR<SparseMatrix<double>, AMDOrdering<int> >的第二個(gè)模板參數，是一個(gè)矩陣重排列（ordering）的方法，為什么要重排列呢，wikipedia的LU分解詞條給了一個(gè)例子可以大概解釋一下，某些矩陣沒(méi)有重排直接分解可能會(huì )失敗。Eigen提供了三種重排列方法，參見(jiàn)OrderingMethods Module。關(guān)于矩陣重排列的細節求數值計算牛人指點(diǎn)！我一般就隨便選一個(gè)填進(jìn)去了>_<。

除了解方程，這個(gè)QR實(shí)例也可以用下面代碼返回Q和R矩陣：

 SparseMatrix < double >   Q ,   R ; 
 qr . matrixQ ( ) . evalTo ( Q ) ; 
 R   =   qr . matrixR ( ) ; 

注意到Q和R的返回方法不一樣，猜測是因為 matrixQ() 成員好像是沒(méi)有完整保存Q矩陣（元素太多？）。

用 Eigen::SPQR 調用SuiteSparseQR進(jìn)行QR分解

SuiteSparseQR效率很高，但是C風(fēng)格接口比較不好用，Eigen提供了 Eigen::SPQR 的接口封裝比如和上面同樣的程序可以這樣寫(xiě)：

 #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 ; 
      // 一個(gè)QR分解的實(shí)例 
      SPQR < SparseMatrix < double >   >   qr ; 
      // 計算分解 
      qr . compute ( A ) ; 
      // 求一個(gè)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 ; 
 } 

如果你有用過(guò)SuiteSparseQR的話(huà)，會(huì )覺(jué)得這個(gè)接口真心好用多了。編譯這個(gè)程序除了spqr庫還需要鏈接blas庫、lapack庫、cholmod庫（SuiteSparse的另一組件），有一點(diǎn)麻煩。比如我在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 ，和矩陣庫沒(méi)關(guān)系。

一些效率的經(jīng)驗

SuiteSparseQR畢竟實(shí)現更好一些，我的一些經(jīng)驗是比自帶Eigen::SparseQR快50%左右吧。

欧美性猛交XXXX免费看蜜桃,成人网18免费韩国,亚洲国产成人精品区综合,欧美日韩一区二区三区高清不卡,亚洲综合一区二区精品久久

基本稀疏矩陣操作

用Eigen進(jìn)行矩陣分解

用 Eigen::SPQR 調用SuiteSparseQR進(jìn)行QR分解

一些效率的經(jīng)驗

相關(guān)文章