求解大型 (n>1000000) 线性方程组答案

【问题标题】：Solving a large (n>1000000) linear system of equations求解大型 (n>1000000) 线性方程组
【发布时间】：2020-05-29 08:21:07
【问题描述】：

我已经达到了以下问题的极限：

作为我的 FEA 代码的一部分（在 MATLAB 中），我需要找到 x，x=A\b。 A 和 b 分别是稀疏、复杂、双精度矩阵和向量。 A 的大小为 (n,n)，b 的大小为 (n,1)，其中 n 为 850000，可以增加到 2000000。此外，A 是对称的，主要是对角线。

我有两个 HPC 集群可供我使用，一个具有 5.7TB 的 RAM，另一个具有 1.5TB（但内核速度更快）。从技术上讲，我可以按原样解决系统，只需等待大约 15 天。但是，我需要对方程组进行求解，每次模拟最多 10 倍。因此，15 天是不可接受的时间量。

我尝试了迭代方法，但这些方法并没有产生与反斜杠运算符相同的结果。在某些情况下也没有获得收敛。

我已将 x=A\b 部分转换为 mexa64 格式，以减少时间。但是，我担心这仍然需要几天时间。

关于如何解决这个问题的任何建议？是否有任何商业代码可以更快/更有效地做到这一点？当模型的节点数超过 1m 时，商业 FEA 包如何解决这个问题？

非常感谢任何帮助！提前致谢。

【问题讨论】：

一旦达到这个级别，您需要特定问题的代码来加速。通常，大规模问题不会将矩阵A 存储在内存中，这可以带来巨大的速度提升。可以这样做是因为操作A*x 和 `A^t*b¬ 可以即时计算。例如我制作了一个断层扫描工具箱，可以动态计算这些操作，这意味着不需要内存，并且从矩阵存储解决方案中加速超过 1000。但它仅适用于特定类型的断层扫描。你需要找到适合你的情况
顺便说一句，matlabs A\b 在后台使用高度加速的代码，我怀疑你可以比使用相同方法提供的 MATLAB 更快
你需要解几个方程组。 A 矩阵对他们中的一些人来说是通用的吗？任何可以使用的结构或限制？否则，我同意 Ander 的观点：在求解线性系统方面似乎很难超越 Matlab 的速度
是的，矩阵 A 是对称的，并且大部分是对角线。我也相应地更新了这个问题。我知道 A\b 是在 MATLAB 中解决问题的一种极其有效的方式。但是，商业 FEA 包在处理此类问题时如何解决它。
如果您需要更快的速度，您可能不得不使用迭代求解器。如果您尝试过的任何迭代求解器收敛太慢或根本不收敛（即使它适合您拥有的矩阵类型），那么您需要找到更好的预处理器。这取决于问题的结构、矩阵的类型以及您可以承担的任何领域知识。

标签： matlab large-data inverse finite-element-analysis

【解决方案1】：

如果您的矩阵是 Hermitian (Aij = Aji*) 并且是正定矩阵，那么使用共轭梯度方法就不会出错。如果你有一个右手边并且保证收敛，那几乎总是最快的方法。此外，这允许您使用可能是有益的猜测解决方案（例如，从前一个时间步长）。虽然解决方案与直接分解不同（它怎么可能没有无数的系数），但它可能并没有那不同。

如果它实际上是一个对称复矩阵 (Aij = Aji)，那么这很奇怪，但可能仍然比一般的复矩阵好。无论哪种情况，请尝试使用只需要您传递一半矩阵的例程（我认为 matlab 不支持此功能）。这会将您的内存使用量减少一半，并且可能会使事情变得更快（移动 16 GB 的数字比 8GB 慢）。

最后，我找不到关于 Matlab 是否支持英特尔 MKL Inspector-executor 稀疏 BLAS 例程的答案。根据我的经验，在我的瓶颈中使用 MKL 总是可以提高性能。它是封闭源代码，你不知道里面发生了什么——可能是魔法。

【讨论】：

【解决方案2】：

稀疏矩阵计算本身就是一个很大的领域。您唯一能做的就是更具体地使用它并使用更可定制的解决方案，因为 Matlab 不允许您这样做并查看它的源代码，看看您是否可以做任何改进。据我所知，主要有两种解决稀疏线性系统的方法，直接和迭代。对于每种类型，都有许多算法。例如，Cholesky 分解具有非常好的性能和低内存使用率，但仅适用于 SPD（对称正定）矩阵。其他方法有LU、LDL、QR等，每一种都满足特定的矩阵类型来求解。

Matlab 也在后台使用其中一种方法（不确定，但我在某处读到它在 A/b 运算符下使用了 SuiteSparse 代码）。根据矩阵的类型（无论是Symmetric Positive Definite 还是Non-symetric 等），您需要选择一个模式（我的意思是迭代或直接）。我认为直接方法更适合这样的大型系统（因为舍入误差）。此外，还有基于 GPU 的求解器可用于直接和迭代方法并且高度并行化，但可能不适合您的大问题，即被分解的矩阵本身超过 8GB。

你应该做的步骤：

确定矩阵的类型（SPD、对称或一般）。
为您的矩阵选择最佳的直接方法。
找到一个可以处理这种分解的免费或商业软件包（例如Eigen、Intel MKL、SuitSpare），并优化代码也支持并行计算。
让你的手有点脏，编写一个非常小的程序，从文件中获取矩阵A 和向量b，并使用其中一个包和例如 C++ 代码找到x（他们已经准备好使用示例代码）。

另外请记住，每个方法和包都可以在开始计算之前改变矩阵的顺序，选择一个好的顺序会对速度甚至准确性产生很大的影响。

这个链接应该也很有趣：

Do not invert that matrix

还有两本关于稀疏矩阵的书：

Sergio Pissanetzky 的稀疏矩阵技术
Timothy A. Davis 的稀疏线性系统的直接方法

【讨论】：