Eigen3.2稀疏矩阵入门

Eigen3.2稀疏矩阵入门
转自 /cvnote/blog/166980
Eigen 自带的稀疏矩阵分解功能包括LDLt 、LLt 分解（即Cholesky 分解，这个功能是LGPL 许可，不是Eigen 的MPL 许可）、LU 分解、QR 分解（这个是3.2版本之后正式Release 的）、共轭梯度解矩阵等。

另外还提供了到第三方稀疏矩阵库的C++接口，包括著名的SuiteSparse 系列（这个系列非常强大，有机会要好好研究一下）的SparseQR 、UmfPack 等。

(欢迎访问计算机视觉研究笔记 和关注新浪微博@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 #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 ;
19 20
21 22
23
24
25
26
27
28
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 #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 ;
25
26
27
28
29
30
31
32
} 这样就用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
#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 ] ) ;
13 14
15 16 17 18 19 20 21 22 23 24 25 26
27
28
29
30
31
32
// 初始化稀疏矩阵
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。