mpmath 矩阵求逆的替代方案或加速答案

【问题标题】：alternatives or speedups for mpmath matrix inversionmpmath 矩阵求逆的替代方案或加速
【发布时间】：2013-03-10 13:21:45
【问题描述】：

我正在用 python 编写一些代码，这些代码需要经常反转大型方阵（100-200 行/列）。

我正在达到机器精度的极限，因此开始尝试使用mpmath 进行任意精度矩阵求逆，但速度非常慢，即使使用gmpy。

在精度为 30（十进制）的情况下反转大小为 20、30、60 的随机矩阵需要约 0.19、0.60 和 4.61 秒，而在 mathematica 中的相同操作需要 0.0084、0.015 和 0.055 秒。

这是在arch linux机器上使用python3和mpmath 0.17（不确定gmpy版本）。我不确定为什么 mpmath 慢得多，但是是否有任何开源库可以接近数学为此管理的速度（即使快 1/2 也不错）？

我不需要任意精度——128 位可能就足够了。我也只是不明白 mpmath 怎么会这么慢。它必须使用非常不同的矩阵求逆算法。具体来说，我使用的是M**-1。

有没有办法让它使用更快的算法或加速它。

【问题讨论】：

你只是用矩阵求逆来解一组方程吗？如果是这样，那么有更有效的方法不需要明确要求逆。 LU 分解我相信效率更高。
不，我在线性规划问题的变体中使用它，所以我需要逆来明确确定成本函数。事实上，问题在于，随着成本变得非常小，不精确的逆可能会导致各种问题。但我认为达到 128 位精度就足够了（至少对于我目前的目的而言）。
当然我从不需要实际的逆，而是需要它乘以其他一些矩阵。所以它与求解 A.x=b 不太相似，因为我需要 A^-1 * b 和 b 是矩阵而不是向量。但也许有一种方法可以概括为矩阵找到这样的解决方案？ OTOH，我需要这样做很多次，所以找到逆可能真的会更好。
为许多 b 评估 A^-1 * b 正是 LU 分解的优点。这同样适用于评估 A^-1 * B 其中 B 矩阵，这与评估矩阵 B 的每一列 b 的 A^-1 * b 相同...
我需要进一步研究这种可能性。你有好的参考吗？

标签： python python-3.x matrix inversion mpmath

【解决方案1】：

不幸的是，mpmath 中的线性代数相当慢。有许多库可以更好地解决这个问题（例如Sage）。也就是说，作为 Stuart 建议的后续，在 Python 中进行相当快速的高精度矩阵乘法相当容易，无需安装任何库，使用定点算术。这是一个使用 mpmath 矩阵进行输入和输出的版本：

def fixmul(A, B, prec):
    m = A.rows; p = B.rows; n = B.cols;
    A = [[A[i,j].to_fixed(prec) for j in range(p)] for i in range(m)]
    B = [[B[i,j].to_fixed(prec) for j in range(n)] for i in range(p)]
    C = [([0] * n) for r in range(m)]
    for i in range(m):
        for j in range(n):
            s = 0
            for k in range(p):
                s += A[i][k] * B[k][j]
            C[i][j] = s
    return mp.matrix(C) * mpf(2)**(-2*prec)

在 256 位精度下，这将两个 200x200 矩阵相乘比我的 mpmath 快 16 倍。这样直接写一个矩阵求逆程序也不难。当然，如果矩阵条目非常大或非常小，您需要先重新调整它们。一个更可靠的解决方案是使用gmpy 中的浮点类型编写自己的矩阵函数，这应该基本上一样快。

【讨论】：

谢谢，但我需要的加速比要多得多。即使对于大小为 100-200 的矩阵，我也需要将反演步骤保持在几分之一秒内。我在 cmets 中发布给 Stuart 的混合 128 位解决方案在 N=150 的情况下进行了约 1/100 秒的反转（而 mpmath 将花费更多更长的时间），增益为小数点后 3-4地方精度。尽管我希望 float128 真的是 128 位，但这可能就足够了。我想我可以在 gmp 中实现上述内容，或者只使用 C 并将其插入 python。
好吧，我在 Mathematica 中以 128 位精度将 200x200 矩阵乘法计时为 4.4 秒，在 Sage 中计时为 3.3 秒（全矩阵求逆的成本大致相同），而仅使用上述代码进行一次矩阵乘法需要 2.7 秒。我估计在 C 中使用 GMP 或 MPFR 将把它推低到大约 0.5 到 1 秒。如果您需要它在几分之一秒内发生，最好的办法可能是研究双双或四双算术。

【解决方案2】：

我假设双精度对最终结果的精度来说不是问题，但对于某些矩阵，它会导致逆的中间结果出现问题。在这种情况下，让我们将普通 numpy（双精度）逆的结果视为一个很好的近似值，然后以此为起点，进行几次牛顿法迭代来求解逆。

让 A 是我们尝试求逆的矩阵，X 是我们对逆矩阵的估计。牛顿方法的一个迭代简单地包括：

X = X*(2I - AX)

对于大型矩阵，计算上述几次迭代的工作与找到逆矩阵的工作相比几乎微不足道，并且可以大大提高最终结果的准确性。试试这个。

顺便说一句，I 是上式中的单位矩阵。

编辑以添加代码以测试浮动类型的精度。

使用此代码测试浮点类型的精度。

x = float128('1.0')
wun = x
two = wun + wun
cnt = 1
while True:
   x = x/two
   y = wun + x
   if y<=wun: break
   cnt +=1

print 'The effective number of mantissa bits is', cnt

【讨论】：

好吧，我不太确定。 numpy 产生的逆矩阵乘以原始矩阵，具有 1e-15 或 1e-16 阶的零（正如预期的那样）。我认为这实际上还不够好。线性程序涉及换出矩阵的一行并重新评估成本，我发现在某些时候变化变得如此之小，成本不再变化（即使它应该）。我认为这是精度低的结果。但也许我可以使用您上面的建议来加快 mp.math 的逆计算。
它可能比 mp.A**(-1) 快一点。我现在刚刚在一个随机的 200x200 矩阵上对其进行了测试。正常的浮点 numpy inv(A) 很快，但正如您所指出的，mpmath 等价物非常慢。然后我使用 mpmath 测试了上述算法的一次迭代（从正常的 inv(A) 结果开始）并将 mp.norm(A*X - mp.eye(200)) 从大约 1E-12 提高到 1E- 18.所以看起来改进是可能的。唯一的问题是，对于这么大的矩阵，即使只是上述迭代的 mpmath 矩阵乘法也很慢（但仍比 mpmath 逆矩阵快）
谢谢，这实际上看起来很有可能。精度似乎提高得非常快......即使一次迭代似乎也能让我达到相当高的精度。不幸的是，正如您所说，mpmath 乘法仍然很慢，但也许我可以找到其他一些可以进行快速高精度矩阵乘法的库。感谢您的建议！
所以我将这个想法与numpy内置的float128支持结合起来，得到了一个中间解决方案。在我的机器上，float128 似乎支持大约 1e-19 或 1e-20 的精度（大约是 1e3 或 1e4 的两倍多）。不幸的是 numpy.linalg 不支持 float128 所以我不能直接反转 128 位。但是在 64 位中反转然后转换为 128 位并使用上面建议的牛顿方法似乎给出了一个很好的折衷方案。对于计算 norm(A*x - I) 时的大型矩阵（约 100 行），我得到 1e-16（128 位）而不是 1e-13（64 位答案）。这也比 mpmath 快很多。
我应该添加以上所有内容都是牛顿方法的一次迭代。更多的迭代似乎并没有显着提高精度。

【解决方案3】：

Multiprecision toolbox for MATLAB 使用 128 位精度（Core i7 930）提供以下时序：

20x20 - 0.007 秒

30x30 - 0.019 秒

60x60 - 0.117 秒

200x200 - 3.2 秒

请注意，对于现代 CPU，这些数字要低得多。

【讨论】：