numpy中矩阵求逆的复杂性

计算科学 麻木的 密集矩阵
2021-12-21 00:19:32

我正在求解需要反转密集方阵的微分方程。这个矩阵求逆消耗了我大部分的计算时间,所以我想知道我是否使用了最快的算法。

我目前的选择是numpy.linalg.inv从我的数字中,我看到它按比例缩放O(n3)其中n是行数,所以该方法似乎是高斯消元法。

根据维基百科,有更快的算法可用。有谁知道是否有实现这些的库?

我想知道,为什么 numpy 不使用这些更快的算法?

2个回答

(评论太长了……)

我假设您实际上需要在算法中计算逆。1首先,重要的是要注意这些替代算法实际上并没有声称更快,只是它们具有更好的渐近复杂度(意味着所需的基本运算数量增长更慢)。事实上,在实践中,这些实际上比标准方法慢得多(对于给定的n),原因如下:

  1. O-notation 将常数隐藏在幂的前面n, 它可以是天文数字——大到C1n3可以远小于C2n2.x对于任何n在可预见的将来,任何计算机都可以处理。(例如,Coppersmith–Winograd 算法就是这种情况。)

  2. 复杂性假设每个(算术)运算都花费相同的时间 - 但在实际实践中远非如此:将一堆数字与相同的数字相乘比将相同数量的不同数字相乘要快得多这是因为当前计算的主要瓶颈是将数据放入缓存,而不是对该数据的实际算术运算。因此,可以重新排列以具有第一种情况(称为缓存感知)的算法将比不可能的算法快得多。(例如,Strassen 算法就是这种情况。)

此外,数值稳定性至少与性能一样重要。在这里,标准方法通常会获胜。

出于这个原因,标准的高性能库(BLAS/LAPACK,当你要求 Numpy 计算逆时调用它)通常只实现这种方法。当然,有 Numpy 实现,例如 Strassen 算法,但是O(n3)在组装级别手动调整的算法将完全击败O(n2.x)用高级语言编写的算法,适用于任何合理的矩阵大小。


1但如果我没有指出这很少是真正必要的,那我就错了:任何时候你需要计算一个产品A1b, 你应该解决线性系统Ax=b(例如,使用numpy.linalg.solve)和使用x相反 - 这更加稳定,并且可以完成(取决于矩阵的结构A)快得多如果您需要使用A1多次,您可以预先计算一个因式分解A(这通常是解决方案中最昂贵的部分)并在以后重用它。

您可能应该注意到,深埋在 numpy 源代码中(请参阅https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src) inv 例程尝试调用 dgetrf 函数从您的系统 LAPACK 包中,然后执行原始矩阵的 LU 分解。这在道德上等同于高斯消元法,但可以通过在高性能 BLAS 中使用更快的矩阵乘法算法将复杂度调整为略低。

如果你遵循这条路线,你应该被警告强制整个库链使用新库而不是你的发行版附带的系统是相当复杂的。现代计算机系统的一种替代方法是使用诸如 scaLAPACK 或(在 Python 世界中)petsc4py 之类的包来查看并行化方法。然而,这些通常更适合用作线性代数系统的迭代求解器,而不是应用于直接方法和 PETSc,尤其是针对稀疏系统而不是密集系统。