Eigen稀疏线性求解

在 Eigen 中，当系数矩阵稀疏时，有多种方法可用于求解线性系统。 由于此类矩阵的特殊表示，应特别注意以获得良好的性能。 有关 Eigen 中稀疏矩阵的详细介绍，请参阅稀疏矩阵操作。 此页面列出了 Eigen 中可用的稀疏求解器。 还介绍了所有这些线性求解器共有的主要步骤。 根据矩阵的属性、所需的精度，最终用户能够调整这些步骤以提高其代码的性能。 请注意，不需要深入了解这些步骤背后隐藏的内容：最后一部分介绍了一个基准程序，可以轻松地使用它来深入了解所有可用求解器的性能。

稀疏求解器列表

直接求解

迭代求解

外部求解器的包装器

稀疏求解器概念

所有这些求解器都遵循相同的一般概念。 这是一个典型和一般的例子：

#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 的上三角部分进行求解。 相反的三角形可能是空的或包含任意值。

在需要解决多个具有相同稀疏模式的问题的情况下，“计算”步骤可以分解如下：

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

...

compute() 方法等效于调用analyzePattern() 和factorize()。

每个求解器都提供了一些特定的功能，例如行列式、对因子的访问、迭代控制等。 更多详细信息可在相应类的文档中找到。

计算步骤

在compute()函数中，矩阵通常被分解：LLT用于自伴随矩阵，LDLT用于一般埃尔米特矩阵，LU用于非埃尔米特矩阵，QR用于矩形矩阵。 这些是使用直接求解器的结果。 对于这类求解器，计算步骤进一步细分为analyzePattern() 和factorize()。

analysisPattern() 的目标是对矩阵的非零元素重新排序，以便分解步骤创建较少的填充。 此步骤仅利用矩阵的结构。 因此，此步骤的结果可用于矩阵具有相同结构的其他线性系统。 但请注意，有时，某些外部求解器（如 SuperLU）要求在此步骤中设置矩阵的值，例如平衡矩阵的行和列。 在这种情况下，此步骤的结果不应与其他矩阵一起使用。

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

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

有关可用方法和相关选项的列表，请参阅 OrderingMethods 模块。

在factorize() 中，计算系数矩阵的因子。 每次矩阵的值发生变化时都应调用此步骤。 但是，矩阵的结构模式不应在多次调用之间发生变化。

对于迭代求解器，计算步骤用于最终设置预处理器。 例如，使用 ILUT 预处理器，在此步骤中计算不完整因子 L 和 U。 请记住，基本上，预处理器的目标是通过求解系数矩阵具有更多聚类特征值的修正线性系统来加速迭代方法的收敛。 对于实际问题，迭代求解器应始终与预处理器一起使用。 在 Eigen 中，通过简单地将预处理器作为模板参数添加到迭代求解器对象中来选择预处理器。

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

成员函数 preconditioner() 返回对预处理器的读写引用以直接与其交互。 有关可用方法的列表，请参阅迭代求解器模块和每个类的文档。

求解步骤

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

X = solver.solve(B); 

这里，B 可以是向量或矩阵，其中列形成不同的右侧。 也可以多次调用 solve() 函数，例如，当所有右侧都无法同时使用时。

x1 = solver.solve(b1);

// Get the second right hand side b2

x2 = solver.solve(b2);

// ...

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