解决线性方程组的最有效方法

Question

我有一个(n x n)的对称矩阵A和一个(n x 1)的向量B。基本上，我只需要求解Ax = b以得到x。问题是A可能很大。因此，我正在寻找在C++中解线性方程组最有效的算法。我查看了Eigen库。显然它有一个SVD方法，但我被告知它很慢。解决x=inverse(A)*b也似乎不是最优解。uBLAS更快吗？还有其他更有效的方法吗？谢谢。

编辑：矩阵A是正定且非稀疏的。

- aesir

---

嗯...高斯消元？ - Eutherpy

你的矩阵是否正定？那么，(P)CG可能是一个选项。除此之外，如果你的系统太大，可以考虑迭代方法（如高斯-塞德尔、雅可比等），而不是直接求解器。 - Zeta

你的矩阵平均大小是多少？你的矩阵是稀疏的吗？ - ggael

根据个人经验，稳定性比速度更重要。如果一个快速的求解器产生错误结果，而现实世界中的问题通常不稳定，那么它就没有什么意义了。高斯方法较差，会积累舍入误差，而LU分解却能得到优秀的结果。现在自己去开发殆已无意义，可以尝试所有方法，并寻找正确的结果。 - Hans Passant

这可能是适合您的正确选择：Eigen Cholesky（也许是具有稳定性的带有枢轴变换的变体）。 - Tobias

Answer 1

解决形如 Ax = b 的线性方程组的最佳方法是按照以下步骤进行：

• 将 A 分解为格式为 A = M1 * M2（其中 M1 和 M2 是三角形的）
• 使用反向代换求解 M1 * y = b，得到 y
• 使用反向代换求解 M2 * x = y，得到 x

对于 方阵 ，第一步将使用LU分解。

对于非方阵，第一步将使用QR分解。

如果矩阵A是正定且不稀疏，则第一步将使用Cholesky分解。

如果您想使用Eigen，则必须首先将其分解，然后三角形求解。

如果仍然很慢，幸运的是，有许多线性代数库可用于提高性能。您应该寻找的程序是dpotrs。以下是一些实现了此功能的库：

• Netlib's LAPACK: 文档 和下载 (免费)
• 英特尔的 MKL: 文档 和下载 (非商业用途免费)
• AMD 的 ACML：下载 (免费)
• PLASMA：下载 (免费，多核优化)
• MAGMA：下载 (免费，实现在CUDA、OpenCL中)
• CULA：下载 (freemium，实现在CUDA中)

如果在整个项目中使用 eigen，您可以根据此处的介绍，将需要的 LAPACK 例程与之接口。

- Pavan Yalamanchili

---

很好的答案。此外还有一个C++头文件库Eigen： Eigen::Matrix<double, N, N> A; Eigen::Matrix<double, N, 1> b; // ... 载入A，b Eigen::Matrix<double, N, 1> x = A.fullPivLu().solve(b); fullPivLU()函数执行A = L*U分解，但你也可以选择其他约10种不同的分解方式。 - wcochran