Eigen教程5 - 求解稀疏线性方程组

Eigen中有一些求解稀疏系数矩阵的线性方程组。由于稀疏矩阵的特殊的表示方式，因此获得较好的性能需要格外注意。查看《Eigen教程3 - 稀疏矩阵操作》，了解更多有关稀疏矩阵的内容。

本文列出了Eigen中的稀疏求解器。同时也介绍了所有线性求解器的共同步骤。

用户可以根据矩阵的性质，准确度的要求，来调整求解步骤，以提升代码的性能。

注意：没有必要深入了解这些步骤后面隐藏的内容。

本文最后一部分列出了一个基准例程，可以用于窥探所有的求解器的性能。

稀疏求解器的列表

Eigen提供了一系列的内建求解器和封装了一些外部求解器库。

内建的直接求解器

类	头文件	求解器类型	矩阵类型	关于性能的特性	许可证	备注
SimplicialLLT	include	直接LLT分解	SPD	Fill-in reducing	LGPL	一般使用它
SimplicialLDLT	include	直接LDLT分解	SPD	Fill-in reducing	LGPL	当矩阵非常稀疏，且问题不太大时使用它
SparseLU	include	LU分解	方阵	Fill-in reducing，Leverage fast dense algebra	MPL2	对于不规则小或大问题进行了优化
SparseQR	include	QR分解	任意，矩形	Fill-in reducing	MPL2	推荐最小方差问题使用，

内建的迭代求解器

类	头文件	求解器类型	矩阵类型	Supported preconditioners, [default]	许可证	备注
ConjugateGradient	include	经典的迭代共轭梯度	SPD	IdentityPreconditioner, [DiagonalPreconditioner], IncompleteCholesky	MPL2	适用于大对称矩阵问题（比如，3D泊松）
LeastSquaresConjugateGradient	include	矩形最小方差问题的共轭梯度	矩形	IdentityPreconditioner, [LeastSquareDiagonalPreconditioner]	MPL2	在不计算A’A的情况下，求解 `$min
BiCGSTAB	include	Iterative stabilized bi-conjugate gradient	方阵	IdentityPreconditioner, [DiagonalPreconditioner], IncompleteLUT	MPL2	为了加速收敛，尝试使用 IncompleteLUT preconditioner

对第三方求解器的封装

类	模块	求解器类型	矩阵类型	有关性能的特性	依赖，许可证	备注
PastixLLT, PastixLDLT, PastixLU	PaStiXSupport	直接LLT，LDLT，LU分解	SPD，SPD，Square	Fill-in reducing, Leverage fast dense algebra, 多线程	需要PaStiX，CeCILL-C	对于难题和对称模式进行了优化
CholmodSupernodalLLT	CholmodSupport	直接LLT分解	SPD	Fill-in reducing, Leverage fast dense algebra	需要SuiteSparse，GPL
UmfPackLU	UmfPackSupport	直接LU分解	Square	Fill-in reducing, Leverage fast dense algebra	需要SuiteSparse，GPL
SuperLU	SuperLUSupport	直接LU分解	方阵	Fill-in reducing, Leverage fast dense algebra	需要SuperLU，(BSD-like)
SPQR	SPQRSupport	QR分解	任意，矩形	fill-in reducing, 多线程，fast dense algebra	需要SuiteSparse，GPL	适用于线性最小方差问题，has a rank-revealing feature
PardisoLLT，PardisoLDLT，PardisoLU	PardisoSupport	直接LLT，LDLT，LU分解	SPD，SPD，方阵	Fill-in reducing, Leverage fast dense algebra, 多线程	需要 Intel MKL，Proprietary	对于难题模式进行了优化，详情查阅 using MKL with Eigen

* 说明：此处的SPD (symmetric positive definite) 表示对称正定。

稀疏求解器的概念

所有的求解器遵循相同的概念，下面是一个典型的例子。

#include <Eigen/RequiredModuleName>
// ...
SparseMatrix<double> A;
// fill A
VectorXd b, x;
// fill b
// solve Ax = b
SolverClassName<SparseMatrix<double> > solver;
solver.compute(A);
if(solver.info()!=Success) {
// decomposition failed
return;
}
x = solver.solve(b);
if(solver.info()!=Success) {
// solving failed
return;
}
// solve for another right hand side:
x1 = solver.solve(b1);

对于SPD求解器，允许指定第二个模板参数来决定使用哪一个三角部分（上三角或下三角）。比如：

ConjugateGradient<SparseMatrix<double>, Eigen::Upper> solver;
x = solver.compute(A).solve(b);

上面的代码仅使用了A的上三角部分进行求解。相对的三角部分可能为空，也可能包含任意值。

当需要求解具有相同稀疏模式的多个问题时，可以将compute()函数拆解成2步：

SolverClassName<SparseMatrix<double> > solver;
solver.analyzePattern(A);   // 在这一步中，A的数值并未使用
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);
...
A = ...;                    // 修改A中非零元素的值，但并未改变非零模式
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);

compute()方法等价于analyzePattern() 加 factorize()。

每一个求解器提供了一些特性，比如行列式，访问因子和控制迭代等。

绝大多数的迭代求解器，也可以用于matrix-free context。

计算步骤

在compute()函数中，将对矩阵进行分解。LLT（三角矩阵分解）用于自伴随矩阵，LDLT用于Hermitian矩阵，LU用于非Hermitian矩阵，QR用于矩形矩阵。为了更加精确的求解，compute ()函数计算步骤又进一步细分为2步：analyzePattern() 和 factorize()。

analyzePattern()的目的是为了记录矩阵中的非零元素，以便分解步骤中创建更少的fill-in。这一步仅仅使用了矩阵的结构。因此，这一步的结果也可用于具有相同结构的其它线性方程组。注意，一些第三方求解器（比如，SuperLU）在本步中需要矩阵的元素值，比如为了平衡矩阵的行和列。这种情形下，这一步的结果不能用于其它线性方程组（矩阵）。

Eigen提供了有限的几种方法在本步来记录矩阵。内建方法（COLAMD，AMD），第三方方法（METIS）。这些方法可以通过设置求解器的模板参数进行设置：

DirectSolverClassName<SparseMatrix<double>, OrderingMethod<IndexType> > solver;

在factorize()中，计算系数矩阵的因子。只要矩阵的值发生变化，就需要调用该函数。然而在多次调用之间，应该保证矩阵的结构不变化。

对于迭代求解器，compute()步骤是用于设置preconditioner。比如使用ILUT preconditioner的求解器，就是在计算步骤中计算非完备因子L和U。记住一点：一般preconditioner都是为了加速迭代方法收敛的速度。Eigen可以通过给迭代求解器对象添加一个模板参数，来选择一个preconditioner：

IterativeSolverClassName<SparseMatrix<double>, PreconditionerName<SparseMatrix<double> > solver;

成员函数preconditioner()返回一个preconditioner 的可读写的引用。

求解步骤

solve()函数计算线性方程组的解。等号右边可以是一个向量，也可以是多个向量（矩阵）。

X = solver.solve(B);

此处B是一个向量或矩阵。

可以多次调用solve()函数。

x1 = solver.solve(b1);
// Get the second right hand side b2
x2 = solver.solve(b2);

直接求解方法的计算结果是机器最小误差。有时结果没必要如此精确，这时，迭代方法更加适用。可以通过setTolerance()来设置希望的误差。

基准例程

绝大多数时候，你只需知道求解你的线性方程组需要的时间和最适合的求解器是什么。Eigen提供了一个基准例程来实现这个功能。

位置：bench/spbench/

编译命令：

make spbenchsolver

矩阵格式： MatrixMarket Coordinate format

该程序将返回所有求解器的花费时间统计信息

通过使用SparseExtra模块（非官方支持），可以将矩阵和等号右边向量转换成matrix-market格式。

#include <unsupported/Eigen/SparseExtra>
...
Eigen::saveMarket(A, "filename.mtx");
Eigen::saveMarket(A, "filename_SPD.mtx", Eigen::Symmetric); // if A is symmetric-positive-definite
Eigen::saveMarketVector(B, "filename_b.mtx");

参考

http://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html