[Eigen中文文档] 求解稀疏线性系统

news/2025/2/19 4:08:35/

文档总目录

本文目录

      • 稀疏求解器列表
        • 内置直接求解器
        • 内置迭代求解器
        • 外部求解器的包装器
      • 稀疏求解器概念
      • 计算步骤
      • 基准测试例程

英文原文(Solving Sparse Linear Systems)

在Eigen中,有多种方法可用于求解稀疏系数矩阵的线性系统。由于此类矩阵的特殊表示,必须特别小心以获得良好的性能。有关Eigen中稀疏矩阵的详细介绍,请参见“[稀疏矩阵操作](#5.1 稀疏矩阵操作)”。本页列出了Eigen中可用的稀疏求解器。还介绍了所有这些线性求解器共同的主要步骤。根据矩阵的属性、所需的准确度,最终用户可以调整这些步骤以提高其代码的性能。请注意,并不需要深入了解这些步骤背后的内容:最后一节介绍了一个基础例程,可轻松使用以获取所有可用求解器的性能洞察。

稀疏求解器列表

Eigen 目前提供了一组广泛的内置求解器,以及外部求解器库的包装器。它们总结在下表中:

内置直接求解器

类及头文件求解器种类矩阵种类与性能相关的功能备注
SimplicialLLT
#include<Eigen/SparseCholesky>
直接LLt分解对称正定(SPD)通过对矩阵的重排列、消元等操作,来减少生成的新元素,以减少计算时间和空间。(Fill-in reducing)SimplicLDLT 通常更可取
SimplicialLDLT
#include<Eigen/SparseCholesky>
直接 LDLt 分解对称正定(SPD)Fill-in reducing推荐用于非常稀疏且不太大的问题(例如,2D 泊松方程)
SparseLU
#include<Eigen/SparseLU>
LU 分解方阵快速的密集代数运算,包括基于GPU的加速、多线程优化、并行计算等。(Leverage fast dense algebra)、Fill-in reducing针对不规则模式的大小问题进行了优化
SparseQR
#include<Eigen/SparseQR>
QR 分解矩形阵Fill-in reducing推荐用于最小二乘问题,具有基本的秩显示特性。

内置迭代求解器

类及头文件求解器种类矩阵种类默认支持的预处理器备注
ConjugateGradient
#include<Eigen/IterativeLinearSolvers>
经典共轭梯度法(CG)SPD IdentityPreconditioner, [DiagonalPreconditioner], IncompleteCholesky推荐用于大型对称问题(例如 3D 泊松方程)
LeastSquaresConjugateGradient
#include<Eigen/IterativeLinearSolvers>
矩形阵最小二乘问题的 CG矩形阵 IdentityPreconditioner, [LeastSquareDiagonalPreconditioner]求解 $min
BiCGSTAB
#include<Eigen/IterativeLinearSolvers>
迭代稳定双共轭梯度方阵IdentityPreconditioner, [DiagonalPreconditioner], IncompleteLUT要加速收敛,请尝试使用 IncompleteLUT 预调节器。
  • 共轭梯度法:(Conjugate Gradient method,简称CG),它是一种迭代求解线性方程组的方法。CG方法通常用于解决大规模稀疏对称正定线性方程组,具有高效、收敛快等优点。在Classic iterative CG中,每次迭代会得到一个新的解向量,同时利用前一次的残差和搜索方向来求解下一次的迭代解。
  • 稳定双共轭梯度法:是求解非对称矩阵线性方程组(Ax = b)的迭代算法。与传统的共轭梯度法不同,稳定双共轭梯度法可以处理非对称矩阵的情况,并且在求解过程中,该方法可以保持数值稳定性,避免出现数值震荡等问题。该方法的主要思路是采用两条共轭的方向(与传统共轭梯度法相同)来求解线性方程组,但在每次迭代过程中,引入了一个稳定因子,该稳定因子的作用是保持迭代过程的数值稳定性。该稳定因子可以通过计算残差向量和搜索方向向量之间的内积来确定。

外部求解器的包装器

模块求解器类型矩阵类型与性能相关的功能依赖/证书备注
PastixLLT
PastixLDLT
PastixLU
PaStiXSupportDirect LLt, LDLt, LU factorizationsSPD
SPD
方阵
Fill-in reducing, Leverage fast dense algebra, MultithreadingRequires the PaStiX package, CeCILL-C针对棘手问题和对称模式进行了优化
CholmodSupernodalLLTCholmodSupportDirect LLt factorizationSPDFill-in reducing, Leverage fast dense algebraRequires the SuiteSparse package, GPL
UmfPackLU UmfPackSupportDirect LU factorization方阵Fill-in reducing, Leverage fast dense algebraRequires the SuiteSparse package, GPL
KLU KLUSupportDirect LU factorization方阵Fill-in reducing, suitted for circuit simulationRequires the SuiteSparse package, GPL
SuperLUSuperLUSupportDirect LU factorization方阵Fill-in reducing, Leverage fast dense algebraRequires the SuperLU library, (BSD-like)
SPQR SPQRSupportQR factorization矩形阵fill-in reducing, multithreaded, fast dense algebraRequires the SuiteSparse package, GPL推荐用于最小二乘问题,具有基本的秩显示特性。
PardisoLLT
PardisoLDLT
PardisoLU
PardisoSupportDirect LLt, LDLt, LU factorizationsSPD
SPD
方阵
Fill-in reducing, Leverage fast dense algebra, MultithreadingRequires the Intel MKL package, Proprietary针对棘手问题模式进行了优化,另请参阅将 MKL 与 Eigen 结合使用
AccelerateLLT
AccelerateLDLT
AccelerateQR
AccelerateSupportDirect LLt, LDLt, QR factorizationsSPD
SPD
方阵
Fill-in reducing, Leverage fast dense algebra, MultithreadingRequires the Apple Accelerate package, Proprietary

这里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 failedreturn;
}
x = solver.solve(b);
if(solver.info()!=Success) {// solving failedreturn;
}
// 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);
...
// modify the values of the nonzeros of A, the nonzeros pattern must stay unchanged
A = ...; 
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);
...

compute() 方法相当于同时调用 analyzePattern()factorize()

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

最后,大多数迭代求解器也可以在无矩阵上下文中使用,请参阅以下示例。

计算步骤

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 来编译该例程。使用 --help 选项运行它以获取所有可用选项的列表。基本上,要测试的矩阵应该遵循 MatrixMarket 坐标格式,并且该例程返回Eigen中所有可用求解器的统计信息。

要以 MatrixMarket 格式导出矩阵和右侧向量,可以使用不受支持的 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 统计信息的示例:

在这里插入图片描述


http://www.ppmy.cn/news/531180.html

相关文章

【谈古论今】观HP Origins道科技企业文化

目录 1 The purpose of a company&#xff08;公司的宗旨&#xff09; 2 Contribution&#xff08;贡献&#xff09; 3 Good friends&#xff08;革命战友情谊&#xff09; 4 Management by walking around&#xff08;亲临现场管理工作&#xff09; 5 Death of the slide…

惠普渠道重新回归

新生力量加入&#xff0c;老兵谱写新传。面对IT新型态带来的机遇与挑战&#xff0c;惠普渠道在聚齐精兵强将后&#xff0c; 凭借充满吸引力的渠道计划和奖励政策&#xff0c;以及全新的面向IT新型态的销售工具&#xff0c;宣告那个久违了的强大的惠普渠道正重新回归。 “IT新型…

企业租赁idc服务器的服务商选择

企业租赁idc服务器的服务商选择 服务器租赁主机托管广州市主机托管中国香港服务器租赁服务器机柜大网络带宽器租赁Https证书云服务器 1.产品介绍 携手腾佑在全国各地布署经营多个营运商等级大数据中心&#xff0c;遮盖华北地区、华东地区、华南地区、西北等地区等地区。 2.产品…

[转]惠普知识管理的前世今生

这是从Google资讯订阅上获取的&#xff0c;已精读一遍 惠普知识管理的前世今生 2008-8-18 8:48:25 佚名/KMCenter 一个新人进来&#xff0c;没有人告诉他该去见什么样的客户&#xff0c;见客户的时候应该做什么&#xff0c;应该说什么&#xff0c;很多人都是凭自己的感觉&…

惠普公司将拆成两家公司

惠普公司的邮件称&#xff1a; 2015 年 11 月 1 日&#xff0c;中国惠普公司( Hewlett-Packard Company )将合法拆分成两家上市公司&#xff1a;其中一家公司将由惠普具备市场领先地位的企业技术基础设施、软件和服务业务构成&#xff0c;并以惠普企业公司( Hewlett Packard …

企业应用软件商的平台怎么挣钱?

我昨天讲了一句话&#xff1a; 所谓平台型公司&#xff0c;一定要下沉做。就跟俄罗斯方块一样&#xff0c;你越在底层&#xff0c;你越是平台。你越垒到最高处了&#xff0c;上面已经没空间了&#xff0c;你就不可能成为平台。 所以大家会看到微软Office全球销量最高客户群最大…

HP 性能测试工具LoadRunner 12.00 中的新增功能详解:附下载地址

HP 性能测试工具LoadRunner 12.00 中的新增功能详解&#xff1a;附下载地址 LoadRunner 12.00 中的新增功能 支持云上的 Load Generator 可直接从 Controller 配置云 Load Generator 可在 Amazon EC2 Cloud 中的 Load Generator 上运行测试 改进了 Controller 和 Load Generato…

微软反盗版瞄准16经销商:不排除法律诉讼可能

“91%被篡改默认安全设置。”“59%在销售之前已经感染了恶意软件。”“72%已被更改浏览器设置。” 称16家电脑经销商预装盗版Windows&#xff0c;不排除法律诉讼可能 新京报讯 (记者 林其玲) 从2008年“黑屏”行动以来&#xff0c;微软在华的反盗版行动不断升级。昨天&#x…