大矩阵求逆方法答案

【问题标题】：Large matrix inversion methods大矩阵求逆方法
【发布时间】：2010-07-28 20:24:13
【问题描述】：

您好，我一直在做一些关于矩阵求逆（线性代数）的研究，我想使用 C++ 模板编程算法，我发现有很多方法，例如：Gauss-Jordan Elimination 或 LU Decomposition我找到了函数 LU_factorize (c++ boost library)

从程序员或数学家的角度，我想知道是否还有其他方法，哪个更好（优点/缺点）？
如果没有其他更快的方法，boost 库中是否已经存在（矩阵）反转函数？，因为我搜索了很多，没有找到。

【问题讨论】：

看看 wolfram 不得不说：mathworld.wolfram.com/MatrixInverse.html Max.
它是通用矩阵还是有特殊结构（例如三角形、正交、稀疏/密集）？
此外，Boost.ublas 对于除相当小的矩阵之外的任何东西都是 SLOW。正如 Vitor Py 所说，如果你想要一个快速的仅限 C++ 的线性代数库，请查看 Eigen。
@Staffan：不是一般矩阵，为什么你认为uBlas很慢？
应该对哪些方面进行模板化？我在这里的数字类型上有代码模板：github.com/victorliu/RNP/blob/master/inc/LinearSolve.h

标签： c++ algorithm math linear-algebra matrix-inverse

【解决方案1】：

正如您提到的，标准方法是执行 LU 分解，然后求解恒等式。这可以使用 LAPACK 库来实现，例如，使用dgetrf（因子）和dgetri（计算逆）。大多数其他线性代数库具有大致相同的功能。

当矩阵是奇异的或接近奇异的时，有一些较慢的方法可以更优雅地降级，并且出于这个原因使用。例如，如果矩阵可逆，Moore-Penrose pseudoinverse 等于逆矩阵，即使矩阵不可逆也经常有用；它可以使用奇异值分解来计算。

【讨论】：

@ismail：当然，但这并不重要。这是一个库；不管它实际上是用 C 或 Fortran 或 Prolog 或 Scheme 实现的，你都可以从 C 或 C++ 代码中调用它。

【解决方案2】：

我建议你看看Eigen 源代码。

【讨论】：

【解决方案3】：

请在 Google 或 Wikipedia 上查找以下流行语。

首先，确保你真的想要逆向。求解系统确实不需要反转矩阵。矩阵求逆可以通过求解 n 个系统来执行，右手边是单位基向量。所以我将专注于解决系统问题，因为它通常是你想要的。

这取决于“大”的含义。基于分解的方法通常必须存储整个矩阵。分解矩阵后，您可以一次求解多个右手边（从而轻松反转矩阵）。我不会在这里讨论分解方法，因为您可能已经知道它们。

请注意，当一个矩阵很大时，它的条件数很可能接近于零，这意味着该矩阵是“数值不可逆的”。补救措施：预处理。检查维基百科。文章写的很好。

如果矩阵很大，你不想存储它。如果它有很多零，则它是一个稀疏矩阵。它要么具有结构（例如，带对角线，块矩阵，...），并且您具有解决涉及此类矩阵的系统的专门方法，或者它没有。

当你面对一个没有明显结构的稀疏矩阵，或者你不想存储的矩阵时，你必须使用迭代方法。它们只涉及矩阵向量乘法，不需要特定形式的存储：您可以在需要时计算系数，或者以您想要的方式存储非零系数，等等。

方法有：

对于对称定正矩阵：共轭梯度法。简而言之，求解 Ax = b 等于最小化 1/2 x^T A x - x^T b。
双共轭梯度法用于一般矩阵。但不稳定。
最小残差方法，或最佳 GMRES。有关详细信息，请查看维基百科文章。您可能需要在重新启动算法之前试验迭代次数。

最后，您可以对稀疏矩阵执行某种分解，使用专门设计的算法来最小化要存储的非零元素的数量。

【讨论】：

【解决方案4】：

根据矩阵的实际大小，您可能只需要在任何给定时间将一小部分列保留在内存中。这可能需要覆盖对矩阵元素的低级写入和读取操作，我不确定 Eigen 是否会允许您这样做。

对于这些矩阵非常大的非常狭窄的情况，有StlXXL 库设计用于对主要存储在磁盘中的数组进行内存访问

编辑更准确地说，如果您有一个未在可用 RAM 中修复的矩阵，首选方法是执行 blockwise inversion。矩阵被递归分割，直到每个矩阵都适合 RAM（这当然是算法的调整参数）。这里棘手的部分是避免在将矩阵拉入和拉出磁盘时使矩阵的 CPU 饿死以进行反转。这可能需要在适当的并行文件系统中进行调查，因为即使使用 StlXXL，这也可能是主要瓶颈。虽然，让我重复一遍咒语； 过早的优化是所有编程罪恶的根源。这种邪恶只能通过编码、执行和配置文件

的净化仪式来消除

【讨论】：

【解决方案5】：

您可能希望在 LAPACK 周围使用 C++ 包装器。 LAPACK 是非常成熟的代码：经过充分测试、优化等。

Intel Math Kernel Library 就是一个这样的包装器。

【讨论】：

IMKL 是 LAPACK 实现（以及更多），而不是包装器。
我并没有贬低 LAPACK 的意思。我的意思是它是一个 C++ 接口，我认为它是用 Fortran 编写的，除非 Intel 重写了所有代码。
是的 LAPACK 是用 Fortran 编写的，但是 man IMKL 是 399 美元！ :D
ATLAS 是开源的，非常好。据我所知，ACML 是免费的，但如果您想分发它，则需要额外的再分发协议。
感谢 Staffan 的帮助我现在正在查看所有的库和学习