Eigen教程-Solving Sparse Linear Systems

转载http://blog.csdn.net/xuezhisdc/article/details/54634080

• 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 AVectorXd b, x;// fill b// solve Ax = bSolverClassName<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 b2x2 = 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-definiteEigen::saveMarketVector(B, "filename_b.mtx");

参考

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