大矩阵求逆方法

Question

您好，我做了一些关于矩阵求逆（线性代数）的研究，想使用C++模板编程来实现算法。我发现有许多方法，如高斯-约旦消元或LU分解，并在c++ boost库中找到了函数LU_factorize。

1. 我想知道是否有其他方法，从程序员或数学家的角度来看哪个更好（优缺点）？

2. 如果没有其他更快的方法，在boost库中是否已经有一个（矩阵）求逆函数？因为我搜索了很多，没有找到。

- Ismail Marmoush

---

看一下沃尔夫拉姆所说的：http://mathworld.wolfram.com/MatrixInverse.html。最大值。 - Max

它是一般矩阵还是具有特殊结构（例如三角形、正交、稀疏/密集）？ - Staffan

此外，Boost.ublas 对于除了相当小的矩阵以外的任何内容都非常缓慢。正如Vitor Py所说，如果你想要一个快速的仅限于C++的线性代数库，请看一下Eigen。 - Staffan

@Staffan：不，这是一般的矩阵，为什么您认为uBlas很慢？ - Ismail Marmoush

哪些方面应该被模板化？我已经在这里使用数字类型对代码进行了模板化：http://github.com/victorliu/RNP/blob/master/inc/LinearSolve.h - Victor Liu

@ismail marmoush：据我所知，ublas不进行任何块分解，这会导致大矩阵性能下降，因为缓存使用不佳。有一个旧的Eigen基准测试包括ublas。测试的矩阵大小相当小。如果您需要最大的性能和多核支持，您应该查看更经典的BLAS/LAPACK实现（ATLAS、IMKL、ACML等）。 - Staffan

Answer 1

正如您所提到的，标准方法是执行LU分解，然后解出单位矩阵。例如，可以使用LAPACK库来实现此方法，使用dgetrf（分解）和dgetri（计算逆）。大多数其他线性代数库具有大致等效的功能。

还有一些较慢的方法，当矩阵是奇异或接近奇异时，它们会更加平稳，并因此而被使用。例如，Moore-Penrose伪逆等于逆矩阵，如果矩阵可逆，则即使矩阵不可逆也经常有用；可以使用奇异值分解来计算它。

- Stephen Canon

---

@ismail：当然，但实际上这并不重要。它是一个库；无论它是用C、Fortran、Prolog还是Scheme实现的，你都可以从C或C++代码中很好地调用它。 - Stephen Canon

Answer 2

我建议你看一下Eigen源代码。

- Vitor Py

Answer 3

根据矩阵的大小，你可能需要在任何给定时间仅将一小部分列保留在内存中。这可能需要覆盖矩阵元素的低级写入和读取操作，我不确定Eigen是否允许这样做，尽管它是一个相当不错的库。

对于那些矩阵非常大的情况，有一个名为StlXXL的库，专门设计用于访问大多存储在磁盘上的数组。

编辑更加精确的说法是，如果您有一个矩阵无法存储在可用RAM中，首选方法是进行分块逆转。该矩阵被递归地分割，直到每个矩阵都适合RAM（当然，这是算法的调优参数）。这里的棘手部分是要避免在从磁盘中拉出和插入矩阵时使CPU缺乏需要反转的矩阵。这可能需要调查适当的并行文件系统，因为即使使用了StlXXL，这也很可能是主要的瓶颈。尽管如此，让我重申这句口号：“过早的优化是所有编程恶魔的根源”。只有通过“编码、执行和配置文件”的净化仪式才能驱逐这种邪恶。

- lurscher

Answer 4

请先在谷歌或维基百科上搜索以下专业术语。首先，请确保您真的需要求逆矩阵。解决一个系统并不需要求逆矩阵。可以通过解n个带有单位基向量作为右侧的系统来执行矩阵求逆。因此，我将重点讨论解决系统，因为这通常是您想要的。这取决于“大”是什么意思。基于分解的方法通常必须存储整个矩阵。一旦您分解了矩阵，就可以同时解决多个右侧（从而轻松地求逆矩阵）。我不会在此处讨论分解方法，因为您可能已经知道它们。请注意，当矩阵很大时，其条件数很可能接近零，这意味着该矩阵“在数值上无法求逆”。解决方法：预处理。请查看维基百科上的相关文章。如果矩阵很大，则不要存储它。如果它有很多零，则是稀疏矩阵。它可能具有结构（例如，带状对角线，块状矩阵等），并且您可以使用专门用于解决涉及这些矩阵的系统的方法，或者它没有。当您面对一个没有明显结构的稀疏矩阵，或者您不想存储的矩阵时，必须使用迭代方法。它们仅涉及矩阵向量乘法，这不需要特定形式的存储：您可以在需要时计算系数，或以所需的方式存储非零系数等。这些方法包括：对于对称正定矩阵：共轭梯度法。简而言之，解Ax = b相当于最小化1/2 x ^ T Ax-x ^ T b。用于一般矩阵的双共轭梯度法。但不稳定。最小残差方法或最佳GMRES。有关详细信息，请查看维基百科文章。您可能需要尝试多次迭代才能重新启动算法。最后，您可以使用专门设计的算法对稀疏矩阵进行某种因式分解，以最小化要存储的非零元素数量。

- Alexandre C.

Answer 5

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

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

- John D. Cook

---

IMKL是LAPACK的实现（以及更多），而不是包装器。 - Staffan

我并不是想贬低LAPACK。我的意思是，它是一个C++接口，用于连接预计是用Fortran编写的代码，除非英特尔重写了所有内容。 - John D. Cook

是的，LAPACK是用Fortran编写的，但IMKL的价格是399美元!! :D - Ismail Marmoush

ATLAS是一款开源且非常优秀的软件。ACML据我所知是免费的，但如果您想要进行分发，则需要签署额外的再分发协议。 - Staffan

感谢Staffan的帮助，我现在正在查看所有的库并学习。 - Ismail Marmoush