Eigen学习笔记26：求解稀疏线性矩阵

在Eigen中，当系数矩阵是稀疏矩阵的时候，存在这几种方法用于求解线性系统。

Because of the special representation of this class of matrices, special care should be taken in order to get a good performance.

有关Eigen中稀疏矩阵的详细介绍，请参见稀疏矩阵操作。

此页面列出了Eigen中可用的稀疏求解器。还介绍了所有这些线性求解器共有的主要步骤。

Depending on the properties of the matrix, the desired accuracy, the end-user is able to tune those steps in order to improve the performance of its code.

请注意，不需要深入了解这些步骤背后的内容： the last section presents a benchmark routine that can be easily used to get an insight on the performance of all the available solvers.

稀疏求解器列表

Eigen当前提供了应用很广泛的内置求解器，以及外部求解器库的包装器。下表中汇总了它们： They are summarized in the following tables:

内置直接求解器

类	求解器类型	Matrix kind	Features related to performance	License	注意：
SimplicialLLT #include<Eigen/SparseCholesky>	Direct LLt factorization	SPD	Fill-in reducing	LGPL	通常最好使用SimplicialLDLT
SimplicialLDLT #include<Eigen/SparseCholesky>	Direct LDLt factorization	SPD	Fill-in reducing	LGPL	建议用于非常稀疏且不太大的问题（例如2D泊松方程）
SparseLU #include<Eigen/SparseLU>	LU factorization	Square	Fill-in reducing, Leverage fast dense algebra	MPL2	optimized for small and large problems with irregular patterns
SparseQR #include<Eigen/SparseQR>	QR factorization	Any, rectangular	Fill-in reducing	MPL2	r建议用于最小二乘问题, has a basic rank-revealing feature

内置迭代求解器

Class	求解器类型	Matrix kind	Supported preconditioners, [default]	License	注意：
ConjugateGradient #include<Eigen/IterativeLinearSolvers>	Classic iterative CG	SPD	IdentityPreconditioner, [DiagonalPreconditioner], IncompleteCholesky	MPL2	推荐用于大型对称问题（例如3D泊松方程）
LeastSquaresConjugateGradient #include<Eigen/IterativeLinearSolvers>	CG for rectangular least-square problem	Rectangular	IdentityPreconditioner, [LeastSquareDiagonalPreconditioner]	MPL2	Solve for min |A'Ax-b|^2 without forming A'A
BiCGSTAB #include<Eigen/IterativeLinearSolvers>	Iterative stabilized bi-conjugate gradient	Square	IdentityPreconditioner, [DiagonalPreconditioner], IncompleteLUT	MPL2	加速收敛,  try it with the IncompleteLUT preconditioner.

外部求解器的包装器

类	模块	求解器类型	Matrix 类	与性能相关的特征	Dependencies,License	注意
PastixLLT PastixLDLT PastixLU	PaStiXSupport	Direct LLt, LDLt, LU factorizations	SPD SPD Square	Fill-in reducing, Leverage fast dense algebra, Multithreading	Requires the PaStiX package, CeCILL-C	针对棘手的问题和对称模式进行了优化
CholmodSupernodalLLT	CholmodSupport	Direct LLt factorization	SPD	Fill-in reducing, Leverage fast dense algebra	Requires the SuiteSparse package, GPL
UmfPackLU	UmfPackSupport	Direct LU factorization	Square	Fill-in reducing, Leverage fast dense algebra	Requires the SuiteSparse package, GPL
SuperLU	SuperLUSupport	Direct LU factorization	Square	Fill-in reducing, Leverage fast dense algebra	Requires the SuperLU library, (BSD-like)
SPQR	SPQRSupport	QR factorization	Any, rectangular	fill-in reducing, multithreaded, fast dense algebra	requires the SuiteSparse package, GPL	推荐用于线性最小二乘法问题, has a rank-revealing feature
PardisoLLT PardisoLDLT PardisoLU	PardisoSupport	Direct LLt, LDLt, LU factorizations	SPD SPD Square	Fill-in reducing, Leverage fast dense algebra, Multithreading	Requires the Intel MKL package, Proprietary	optimized for tough problems patterns, see also using MKL with Eigen

此处SPD表示对称正定。

稀疏求解器概念

#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求解器，第二个可选模板参数允许指定必须使用哪个三角形部分，例如：

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

 

在以上示例中，仅考虑输入矩阵A的上三角部分进行求解。The opposite triangle might either be empty or contain arbitrary values.

In the case where multiple problems with the same sparsity pattern have to be solved, then the "compute" step can be decomposed as follow:

SolverClassName<SparseMatrix<double> > solver;
solver.analyzePattern(A);   // for this step the numerical values of A are not used
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);
...
A = ...;                    // modify the values of the nonzeros of A, the nonzeros pattern must stay unchanged
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);
...

The compute() method is equivalent to calling both analyzePattern() and factorize().

每个求解器都提供一些特定功能，例如行列式，对因子的访问，迭代控制等。有关更多详细信息，请参见相应类的文档。

Finally, most of the iterative solvers, can also be used in a matrix-free context, see the following example .

计算步骤

在compute()函数中，通常对矩阵进行分解：LLT用于自伴矩阵，LDLT用于hermitian matrices，LU用于non hermitian matrices，QR用于rectangular matrices。这些是使用直接求解器的结果。对于这类求解器，计算步骤可进一步细分为analyzerPattern()和factorize()。

analyticsPattern()的目标是对矩阵的非零元素进行重新排序，以使the factorization step creates less fill-in。此步骤仅利用矩阵的结构。Hence, the results of this step can be used for other linear systems where the matrix has the same structure.但是请注意，有时某些外部求解器（例如SuperLU）要求在此步骤中设置矩阵的值，例如，以平衡矩阵的行和列。在这种情况下，此步骤的结果不应与其他矩阵一起使用。

Eigen提供几种有限的方法来在此步骤中对矩阵进行重新排序，方法是内置的（COLAMD，AMD）或外部的（METIS）。这些方法在求解器的模板参数列表中设置：

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

See the  OrderingMethods module for the list of available methods and the associated options.

In factorize(), the factors of the coefficient matrix are computed. 每当矩阵值更改时，都应调用此步骤。However, the structural pattern of the matrix should not change between multiple calls.

对于迭代求解器，计算步骤用于最终设置预处理器。 例如，使用ILUT预处理器，在此步骤中计算不完全因子L和U。请记住，基本上， the goal of the preconditioner is to speedup the convergence of an iterative method by solving a modified linear system where the coefficient matrix has more clustered eigenvalues. 对于实际问题，迭代求解器应始终与前置条件一起使用。在Eigen中，只需将预处理器作为模板参数添加到迭代求解器对象中，即可选择该预处理器。

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

The member function preconditioner() returns a read-write reference to the preconditioner to directly interact with it.

  有关可用方法的列表，请参见迭代求解器模块和每个类的文档。

求解步骤

solve()函数计算具有一个或多个右侧的线性系统的解。

X = solver.solve(B);

在此，B可以是向量或矩阵，其中各列形成不同的右侧。 The solve() function can be called several times as well, for instance when all the right hand sides are not available at once.

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

对于直接方法，解决方案以机器精度计算。有时，解决方案不必太精确。在这种情况下，迭代方法更合适，并且可以在使用setTolerance（）的求解步骤之前设置所需的精度。有关所有可用功能，请参阅“ 迭代求解器”模块的文档。

基准程序

在大多数情况下，您所需要做的就是了解解决系统所需的时间，并希望有什么是最合适的解决器。在Eigen中，我们提供了可用于此目的的基准例程。It is very easy to use. In the build directory, navigate to bench/spbench and compile the routine by typing make spbenchsolver. 使用–help选项运行它以获取所有可用选项的列表。Basically, the matrices to test should be in MatrixMarket Coordinate format, and the routine returns the statistics from all available solvers in Eigen.

要以矩阵市场格式导出矩阵和右侧矢量，可以使用不受支持的SparseExtra模块：

#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");

下表提供了一些Eigen内置和外部求解器的XML统计信息示例。

 

Matrix	N	NNZ		UMFPACK	SUPERLU	PASTIX LU	BiCGSTAB	BiCGSTAB+ILUT	GMRES+ILUT	LDLT	CHOLMOD LDLT	PASTIX LDLT	LLT	CHOLMOD SP LLT	CHOLMOD LLT	PASTIX LLT	CG
vector_graphics	12855	72069	Compute Time	0.0254549	0.0215677	0.0701827	0.000153388	0.0140107	0.0153709	0.0101601	0.00930502	0.0649689
			Solve Time	0.00337835	0.000951826	0.00484373	0.0374886	0.0046445	0.00847754	0.000541813	0.000293696	0.00485376
			Total Time	0.0288333	0.0225195	0.0750265	0.037642	0.0186552	0.0238484	0.0107019	0.00959871	0.0698227
			Error(Iter)	1.299e-16	2.04207e-16	4.83393e-15	3.94856e-11 (80)	1.03861e-12 (3)	5.81088e-14 (6)	1.97578e-16	1.83927e-16	4.24115e-15
poisson_SPD	19788	308232	Compute Time	0.425026	1.82378	0.617367	0.000478921	1.34001	1.33471	0.796419	0.857573	0.473007	0.814826	0.184719	0.861555	0.470559	0.000458188
			Solve Time	0.0280053	0.0194402	0.0268747	0.249437	0.0548444	0.0926991	0.00850204	0.0053171	0.0258932	0.00874603	0.00578155	0.00530361	0.0248942	0.239093
			Total Time	0.453031	1.84322	0.644241	0.249916	1.39486	1.42741	0.804921	0.862891	0.4989	0.823572	0.190501	0.866859	0.495453	0.239551
			Error(Iter)	4.67146e-16	1.068e-15	1.3397e-15	6.29233e-11 (201)	3.68527e-11 (6)	3.3168e-15 (16)	1.86376e-15	1.31518e-16	1.42593e-15	3.45361e-15	3.14575e-16	2.21723e-15	7.21058e-16	9.06435e-12 (261)
sherman2	1080	23094	Compute Time	0.00631754	0.015052	0.0247514	-	0.0214425	0.0217988
			Solve Time	0.000478424	0.000337998	0.0010291	-	0.00243152	0.00246152
			Total Time	0.00679597	0.01539	0.0257805	-	0.023874	0.0242603
			Error(Iter)	1.83099e-15	8.19351e-15	2.625e-14	1.3678e+69 (1080)	4.1911e-12 (7)	5.0299e-13 (12)
bcsstk01_SPD	48	400	Compute Time	0.000169079	0.00010789	0.000572538	1.425e-06	9.1612e-05	8.3985e-05	5.6489e-05	7.0913e-05	0.000468251	5.7389e-05	8.0212e-05	5.8394e-05	0.000463017	1.333e-06
			Solve Time	1.2288e-05	1.1124e-05	0.000286387	8.5896e-05	1.6381e-05	1.6984e-05	3.095e-06	4.115e-06	0.000325438	3.504e-06	7.369e-06	3.454e-06	0.000294095	6.0516e-05
			Total Time	0.000181367	0.000119014	0.000858925	8.7321e-05	0.000107993	0.000100969	5.9584e-05	7.5028e-05	0.000793689	6.0893e-05	8.7581e-05	6.1848e-05	0.000757112	6.1849e-05
			Error(Iter)	1.03474e-16	2.23046e-16	2.01273e-16	4.87455e-07 (48)	1.03553e-16 (2)	3.55965e-16 (2)	2.48189e-16	1.88808e-16	1.97976e-16	2.37248e-16	1.82701e-16	2.71474e-16	2.11322e-16	3.547e-09 (48)
sherman1	1000	3750	Compute Time	0.00228805	0.00209231	0.00528268	9.846e-06	0.00163522	0.00162155	0.000789259	0.000804495	0.00438269
			Solve Time	0.000213788	9.7983e-05	0.000938831	0.00629835	0.000361764	0.00078794	4.3989e-05	2.5331e-05	0.000917166
			Total Time	0.00250184	0.00219029	0.00622151	0.0063082	0.00199698	0.00240949	0.000833248	0.000829826	0.00529986
			Error(Iter)	1.16839e-16	2.25968e-16	2.59116e-16	3.76779e-11 (248)	4.13343e-11 (4)	2.22347e-14 (10)	2.05861e-16	1.83555e-16	1.02917e-15
young1c	841	4089	Compute Time	0.00235843	0.00217228	0.00568075	1.2735e-05	0.00264866	0.00258236
			Solve Time	0.000329599	0.000168634	0.00080118	0.0534738	0.00187193	0.00450211
			Total Time	0.00268803	0.00234091	0.00648193	0.0534865	0.00452059	0.00708447
			Error(Iter)	1.27029e-16	2.81321e-16	5.0492e-15	8.0507e-11 (706)	3.00447e-12 (8)	1.46532e-12 (16)
mhd1280b	1280	22778	Compute Time	0.00234898	0.00207079	0.00570918	2.5976e-05	0.00302563	0.00298036	0.00144525	0.000919922	0.00426444
			Solve Time	0.00103392	0.000211911	0.00105	0.0110432	0.000628287	0.00392089	0.000138303	6.2446e-05	0.00097564
			Total Time	0.0033829	0.0022827	0.00675918	0.0110692	0.00365392	0.00690124	0.00158355	0.000982368	0.00524008
			Error(Iter)	1.32953e-16	3.08646e-16	6.734e-16	8.83132e-11 (40)	1.51153e-16 (1)	6.08556e-16 (8)	1.89264e-16	1.97477e-16	6.68126e-09
crashbasis	160000	1750416	Compute Time	3.2019	5.7892	15.7573	0.00383515	3.1006	3.09921
			Solve Time	0.261915	0.106225	0.402141	1.49089	0.24888	0.443673
			Total Time	3.46381	5.89542	16.1594	1.49473	3.34948	3.54288
			Error(Iter)	1.76348e-16	4.58395e-16	1.67982e-14	8.64144e-11 (61)	8.5996e-12 (2)	6.04042e-14 (5)