Eigen 3.2稀疏矩阵入门

原创
2013/10/06 14:59
阅读数 8.4K

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%左右吧。

相关文章

C++矩阵运算库推荐

欢迎访问计算机视觉研究笔记cvnote.info和关注新浪微博@cvnote

展开阅读全文
加载中

作者的其它热门文章

打赏
0
3 收藏
分享
打赏
1 评论
3 收藏
0
分享
返回顶部
顶部