Eigen入门教程（5）稀疏矩阵

myleader

Eigen提供了稀疏矩阵的相关功能，包括创建、赋值、子矩阵、分解等一系列功能。

稀疏矩阵部分有3个模块：

SparseCore提供了基本的稀疏矩阵功能，比如基本的线性代数运算等，包括SparseMatrix类和SparseVector类，使用时#include <Eigen/SparseCore>

SparseCholesky提供了直接LLT和LLDT Cholesky分解，用于解自伴随正定稀疏矩阵问题，使用时#include <Eigen/SparseCholesky>

IterativeLinearSolvers提供了大型通用线性方阵的迭代求解器，也适用于自伴随正定稀疏矩阵问题，使用时#include <Eigen/IterativeLinearSolvers>

当然你也可以直接#include <Eigen/Sparse>同时使用这3个模块。

稀疏矩阵的定义

稀疏矩阵的定义方法和满矩阵其实差不多

1. Eigen::SparseMatrix<double> A(m,m);
   
2. SparseMatrix<std::complex<float> > mat(1000,2000)；
   
3. SparseMatrix<double,RowMajor> mat(1000,2000);
   
4. SparseVector<std::complex<float> > vec(1000);
   
5. SparseVector<double,RowMajor> vec(1000);
   

复制代码

稀疏矩阵的赋值

常规矩阵只要存储值就可以了，它在内存中的位置就是它在矩阵中的位置。稀疏矩阵内部数据的存储和常规矩阵是不一样的，由于0元素不会存储在内存中，所以数据在内存中的位置不代表其在矩阵中的位置，通常不同的稀疏矩阵解决方案会使用不同的格式来存储稀疏矩阵的内容，要存储的内容包括：值、所在位置、矩阵的最大尺寸等信息。Eigen的稀疏矩阵类内部使用一种独特的数据结构来存储，而且每个单元存储的内容、含义和尺寸都不一样，所以没有办法使用通用的类来充当其单元。

为了能够输入数据，Eigen引入了Eigen::Triplet<>类，使用这个类再配合专用的读取该类的函数就可以把稀疏矩阵压缩到SparseMatrix类中存储起来。

首先定义一个Eigen::Triplet<>类

1. Eigen::Triplet<double> b(0,0,2);     //Eigen::Triplet<double> b=Eigen::Triplet<double>(0,0,2);

复制代码

其中构造函数第1个参数是行位置，第2个参数是列位置，第3个参数是值，第1和第2个参数必须是正的int型，而且还必须小于稀疏矩阵的尺寸，第3个参数必须和定义类时尖括号里面生命的类型一致。

然后你必须定义一系列Eigen::Triplet<>类并把它们按顺序存放在内存中，最常用的方法当然是std::vector，

1. std::vector<Eigen::Triplet<double> > input;
   
2. for(int i=0;i<somevalue;++i)
   
3. {
   
4. input.pushback(Eigen::Triplet<double>(row,col,value));
   
5. }

复制代码

SparseMatrix类提供了一个setFromTriplets()函数用以从Eigen::Triplet<>的数组构造一个稀疏矩阵，如果这个数列中的元素存在对系数矩阵中同一位置重复赋值的情况，则以第一个值为准

1. A.setFromTriplets(input.begin(), input.end());

复制代码

input当中的Triplet数组就被压缩进稀疏矩阵A中，注意Eigen并不是把Triplet类按顺序排列在内存中，否则这和代成员函数的std::vector还有什么区别呢，实际上Eigen对数据做了更进一步的压缩，内存中的数据排列也不再是简单地顺序排列。具体可以参看
http://eigen.tuxfamily.org/dox/TutorialSparse.html

这种方法不要求Triplet按照矩阵中的位置在内存中顺序存放，而且允许多个元素指向矩阵中同一个位置，这些元素的值会被加在一起，存在矩阵中那个位置。另外A矩阵中原来存有的所有数据都会被销毁，这个稀疏矩阵的名义尺寸也会改变。

这种方法的缺点是显而易见的，input仍然占用着内存，就是说为了利用Eigen::SparseMatrix类的高级计算功能，我们必须要先浪费一点内存，在矩阵构造完成后再使用clear()函数清空这个vector容器。

另外一种方法就是使用insert()函数，注意这个函数和stl容器中的insert含义不同，用法也不一样

1. A.insert(i,j)=val;

复制代码

配合上足够的循环，就可以在节省内存的情况下完成稀疏矩阵的构建了。

但是又有一个问题，如果(i,j)这个位置已经有值了，insert()函数会将其覆盖，除非你在编程时能够确保绝对不会出现这种情况，否则最后的结果很可能不是你所期待的。而且Eigen的稀疏矩阵在内存中并不是简单地把值和位置构成的类按顺序排列，insert()函数每次执行都要把压缩信息调整一次，就是说相比setFromTriplets()增加了额外的计算量。如果(i,j)超过了稀疏矩阵的名义尺寸，会出错。

Eigen提供了另外一个函数coeffRef()来解决这个问题，它会在插入时检查对应的位置是否已经有值了，如果没有才会插入，如果有了则会跳过。这个函数的计算量比insert()还要多，因为多了一步判断那个位置是否已经有值的运算。

对比这3种方法就可以发现setFromTriplets()浪费内存，不过速度却是最快的；insert()和coeffRef()内存占用倒是小了，可是计算量大了；而在insert()和coeffRef()两者之间insert()不够安全，coeffRef()虽然是安全了，但却要付出更多的计算为代价。实际使用中到底该如何取舍取决于具体的问题，由程序员自己把握。

内存预留与压缩

有一种办法可以稍微减少insert()和coeffRef()的计算消耗，就是内存预留。举例来讲，对于流体计算，通常矩阵都是对角占优的，实际上每列往往是3、5或7个有效元素，在这种情况下，我们可以指定该稀疏矩阵每列只占用对应数量的元素，然后去预留内存

1. A.reserve(VectorXd::Constant(n,3));

复制代码

reserve()函数支持一个Eigen::Vector作为参数，该Vector的尺寸必须和稀疏矩阵的列数相等，如果该稀疏矩阵是按行排列的，则必须和稀疏矩阵的行数相等。然后它会在内存中逐行的去预留内存。这与Eigen稀疏矩阵的数据结构有关，但数据结构不是本文主要讲解的内容，所以各位只要知道这么回事就行了（当然，其实楼主对这个东西也不是很懂）。例子中给的是每列都预留3个非0元素的空间，实际使用中，每列可以预留不同的空间。

在完成了内存预留之后，再用insert()或coeffRef()赋值时，稀疏矩阵的解压缩计算中有一部分会跳过，这样可以节省一些计算量，当然，这种节省是以一定的内存占用为代价的，因为我们必须要为可能的数据占用足够的内存，保险起见，有时必须要多占用点。

由于多预留了内存空间，在完成矩阵赋值之后，就必须要把空闲的内存释放掉

1. A.makeCompressed();

复制代码

如果程序员对问题的把握比较准，对稀疏矩阵的元素量把握比较准，那么内存预留不会多太多，这项压缩操作可节省的内存恐怕不会很大，但计算量可能却不小。Eigen的官方手册说不做内存压缩也可以，不过楼主的习惯是每次都做，否则那些没有被赋值的内存就像野指针一样不知道会带来什么问题。

Eigen::SparseMatrix还提供了一个isCompressed()函数来检查是否完成了压缩，它会返回一个bool值

到底是在内存预留时超级精准的预测内存占用还是每次在使用之后都来一次压缩，取决于程序员。

稀疏矩阵的操作

稀疏矩阵也支持cols()、rows()之类的成员函数，也支持resize()，不过resize()只是声明一下这个稀疏矩阵的尺寸，以免在赋值时出现行、列超标的情况，实际上不会在内存中多开辟空间。类似的，即使构造函数给定了这个稀疏矩阵的尺寸，它也不会在内存中开辟那么大的空间，仅仅是给出一个名义尺寸罢了。同样的cols()和rows()也只是返回一个名义值。

稀疏矩阵也支持block函数读取子矩阵，而且返回的还是稀疏矩阵

要注意的是，如果希望提取某一个位置的值，不能直接用A(i,j)的方法，而是要

1. A.coeff(i,j);

复制代码

另外Eigen::SparseMatrix类还提供了大量有用的成员函数，这些成员函数很多都和常规矩阵相同，如conjugate()、cwiseAbs()、cwiseAbs2()、cwiseSqrt()、cwiseProduct()、cwiseQuotient()、cwiseInverse()、cwiseMax(const Scalar &other)、cwiseMax(const Eigen::SparseMatrixBase< OtherDerived > &other)、cwiseMin(const Scalar &other)、cwiseMin(const Eigen::SparseMatrixBase< OtherDerived > &other)、cwiseEqual(const Scalar &s)、cwiseEqual(const Eigen::SparseMatrixBase< OtherDerived > &other)、cwiseNotEqual(const Eigen::SparseMatrixBase< OtherDerived > &other)、real()、imag()、swap(SparseMatrix &other)

要注意稀疏矩阵的成员函数支持常规矩阵作为参数，也就是说你想计算一个常规矩阵和稀疏矩阵的积，直接乘就可以。

稀疏矩阵的求解

既然是稀疏矩阵，那肯定不小，所以就不用考虑矩阵求逆的方法了，直接来看分解法

和常规矩阵求解有多种方法类似，Eigen的稀疏矩阵求解也支持支持多种方法。而且Eigen不仅自己提供了一些求解器，还支持外部求解器。

类	模块	求解器类型	适用矩阵	性能有关的特性	许可证和外部依赖	注意事项
SimplicialLLT	#include <Eigen/SparseCholesky>	直接LLT分解	SPD	Fill-in reducing	Eigen内建LGPL	推荐求解器，解决一般问题
SimplicialLDLT	#include <Eigen/SparseCholesky>	直接LLT分解	SPD	Fill-in reducing	Eigen内建LGPL	适用于非常稀疏，而且又不太大的问题
ConjugateGradient	#include <Eigen/IterativeLinearSolvers>	古典迭代法	SPD	Preconditionning	Eigen内建LGPL	推荐用于较大的对称问题
BiCGSTAB	#include <Eigen/IterativeLinearSolvers>	共轭梯度迭代法	方阵	Preconditionning	Eigen内建LGPL	不一定总收敛
PastixLLT PastixLDLT PastixLU	#include <Eigen/PaStiXSupport>	直接LLT、LLDT、LU分解	SPD SPD 方阵	Fill-in reducing, Leverage fast dense algebra, Multithreading	PaStiX软件包CeCILL-C许可协议	适用于较难的和对称的问题
CholmodSupernodalLLT	#include <Eigen/CholmodSupport>	直接LLT分解	SPD	Fill-in reducing, Leverage fast dense algebra	SuiteSparse软件包GPL
UmfPackLU	#include <Eigen/UmfPackSupport>	直接LU分解	方阵	Fill-in reducing, Leverage fast dense algebra	SuiteSparse软件包GPL
SuperLU	#include <Eigen/SuperLUSupport>	直接LU分解	方阵	Fill-in reducing, Leverage fast dense algebra	SuperLU软件包，类BSD协议

SPD的全文是symmetry positive definite就是对称正定矩阵，不过即使你的矩阵不是严格的“正定”或者严格的“对称”其实前3种方法都可以。

要注意的是，Eigen除了提供了自带的4个求解器以外，还可以通过链接外部库来调用其他求解器，这就要求你的编译器在编译时能找到对应的头文件，在链接是能找到对应的库文件，如果是动态链接，还要求在运行时能找到对应的动态库。UmfPack要怎么安装不是本文讨论的内容，以后有时间楼主再开帖讨论。

有经验的研究者都非常了解UmfPack目前仍然是稀疏矩阵求解的翘楚，不过它是严格的GPL许可的，所以说如果你使用UmfPack，那你就必须要开放你自己的源代码，这可能会让很多人难以接受，那就选择其他许可证不严格的求解器吧，你可能要付出求解速度和收敛稳定性的代价，不过总还算是保住了自己的代码。

在Eigen入门教程（4）中，常规矩阵的分解可以通过Matrix类的成员函数调用分解器，生成分解矩阵，然后调用solve()函数来求解，稀疏矩阵求解器没有这项功能。因为稀疏矩阵的求解比常规矩阵的求解更加的不确定，所以你只能通过构造求解器来求解，这样有助于维持程序的稳定性。

而稀疏矩阵构造求解器方法来解方程的流程和常规矩阵没有太大差别：

在包含了对应的头文件，并声明了命名空间以后。首先你需要构造一个求解器

1. SimplicialLDLT<SparseMatrix<double> > solver;

复制代码

然后读入要求解的矩阵

1. solver.compute(A);

复制代码

接下来的这步和常规矩阵的不太一样，你得判断一下矩阵的分解有没有成功

1. if(solver.info()!=Success)
   
2. {
   
3.   // decomposition failed
   
4.     return(false);
   
5. }

复制代码

注意Success是Eigen内建的

然后求解

1. c=solver.solve(d);

复制代码

要命的是，你还得再判断一下矩阵的求解有没有成功

1. if(solver.info()!=Success)
   
2. {
   
3.   // decomposition failed
   
4.     return(false);
   
5. }

复制代码

如果没问题，那c就是矩阵的解。

这段代码是以SimplicialLDLT为例讲解的，其他的求解器类似。

要注意的是，和常规矩阵类似，方程组的求解总会有一些不确定性，比如不能收敛、无解之类的，这就要求程序员在编制程序进行运算时多考虑一些可能性，有些时候可能要试试其他的方法，说不定就搞定了，当然有时候即使看起来搞定了，实际上也可能根本没搞定，所以这里就要程序员更有耐心和细心的来处理这些问题，而不是把抱怨停留在口头上。

对于某些特殊形式的稀疏矩阵，如主对角占优的矩阵，有一些已经被历史证明的非常快速有效的专用算法如TDMA算法，那些算法的速度往往比这些通用算法的速度要快，只不过可能需要你自己多写一些代码。如果是MATLAB，多出来的for循环可能把算法带来的速度优势完全消耗掉了，而这里是C++编程，循环不会占用太多的计算，不过你可能无法调用矢量指令集来加速。到底应该选择哪种取决于程序员自身，以楼主自己的经验来看，还是那些久经考验的专用算法更快，所以建议大家多敲几行代码，不吃亏的。

稀疏矩阵这个东西楼主本人也不是十分懂，也只能写这么多了，如果有谁可以补充，欢迎在下面跟帖

下一篇计划是几何运算，这部分内容不太多，新年要到了，楼主也要歇一歇了

hoop247

版主好惨啊，都没人回复你！我来回复！Eigen确实是个有意思的数学库！