C++ 求解线性方程组的最有效方法

Most efficient way to solve a system of linear equations

提问by aesir

I have an (n x n) symmetric matrix A and a (n x 1) vector B. Basically, I just need to solve Ax = b for x. The issue is that A is going to potentially be massive. So I'm looking for the most efficient algorithm for solving linear equations in C++. I looked over the Eigen library. Apparently it has an SVD method, but I've been told it's slow. Solving x=inverse(A)*b also seems like it would be suboptimal. Is uBLAS faster? Are there any more efficient methods? Thanks.

我有一个 (nxn) 对称矩阵 A 和一个 (nx 1) 向量 B。基本上，我只需要为 x 求解 Ax = b。问题是 A 可能会很大。所以我正在寻找在 C++ 中求解线性方程的最有效算法。我查看了 Eigen 图书馆。显然它有一个 SVD 方法，但我被告知它很慢。求解 x=inverse(A)*b 似乎也不是最理想的。uBLAS 更快吗？有没有更有效的方法？谢谢。

Edit: matrix A is positive definite and not sparse.

编辑：矩阵 A 是正定的而不是稀疏的。

回答by Pavan Yalamanchili

The best way to solve a system of linear equations of the form Ax = bis to do the following.

求解该形式的线性方程组的最佳方法Ax = b是执行以下操作。

• decompose Ainto the format A = M1 * M2(where M1and M2are triangular)
• Solve M1 * y = bfor yusing back substitution
• Solve M2 * x = yfor xusing back substitution

• 分解A成格式A = M1 * M2（其中M1和M2是三角形）
• 解决M1 * y = b了y使用回代
• 解决M2 * x = y了x使用回代

For square matrices, step 1 would use LU Decomposition.

对于方阵，步骤 1 将使用LU 分解。

For non square matrices, step 1 would use QR Decomposition.

对于非方阵，步骤 1 将使用QR 分解。

If matrix A is positive definite and not sparseyou'd use Cholesky Decompositionfor the first step.

如果矩阵 A 是正定的且不是稀疏的，则您可以在第一步中使用Cholesky 分解。

If you want to use eigen, you will have to first decompose itand then triangular solveit.

如果要使用本征，则必须先分解它，然后对其进行三角求解。

If this is still slow, thankfully, there are numerous linear algebra libraries available that can help improve the performance. The routine you should be looking for is dpotrs. Some libraries that have this implemented are as follows:

如果这仍然很慢，幸运的是，有许多线性代数库可以帮助提高性能。您应该寻找的例程是dpotrs. 一些实现了这个的库如下：

• Netlib's lapack: Documentationand Download(free)
• Intel's MKL: Documentationand Download(Free for non commercial use).
• AMD's ACML: Download(free)
• PLASMA: Download(free, multi core optimized)
• MAGMA : Download(free, Implemented in CUDA, OpenCL)
• CULA: Download(freemium, Implemented in CUDA).

• Netlib 的 lapack：文档和下载（免费）
• Intel 的 MKL：文档和下载（非商业用途免费）。
• AMD 的 ACML：下载（免费）
• PLASMA：下载（免费，多核优化）
• MAGMA：下载（免费，在 CUDA、OpenCL 中实现）
• CULA：下载（免费增值，在 CUDA 中实现）。

If you are using eigen in the overall project, you can interface the LAPACK routine you need as described here.

如果您在整个项目中使用 eigen，您可以按照此处所述连接您需要的 LAPACK 例程。