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

Question

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

pythonpython-3.xmatrixinversionmpmath

4

我正在使用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为什么如此慢，但是否有任何开源库可以接近mathematica为此管理的速度（即使1/2的速度也很好）？

我不需要任意精度——128位可能已经足够了。我也不明白mpmath为什么会慢那么多。它一定是使用非常不同的矩阵反演算法。具体来说，我正在使用M ** -1。

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

- user2153813

你是在使用矩阵求逆来解决一组方程吗？如果是的话，那么有更有效率的方法，不需要显式地求逆。我相信LU分解会更加高效。 - Stuart

不，我正在使用它来解决线性规划问题的变体，因此我需要逆矩阵来明确确定成本函数。实际上，问题在于随着成本变得非常小，不精确的逆矩阵可能会引起各种问题。但是我认为提高到128位精度就足够了（至少对于我目前的目的而言）。 - user2153813

当然，我从来不需要实际的逆矩阵，而是需要将其乘以其他一些矩阵。因此，它与解A.x=b并不完全类似，因为我需要A^-1 * b，其中b是一个矩阵而不是向量。但也许有一种方法可以推广找到这样的矩阵解吗？另一方面，我需要多次执行此操作，因此找到逆矩阵可能真的更好。 - user2153813

对于许多b的情况下评估A^-1 * b，LU分解非常适用。同样适用于评估A^-1 * B，其中B矩阵，这与为矩阵B的每个列向量b评估A^-1 * b相同... - Fredrik Johansson

我需要更深入地研究这个可能性。你有什么好的参考资料吗？ - user2153813

3个回答

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

- Stuart

嗯，我不太确定。当与原矩阵相乘时，numpy生成的逆矩阵具有1e-15或1e-16的零（如双精度所预期）。我认为这实际上还不够好。线性规划涉及交换矩阵的一行并重新评估成本，我发现在某些点上，变化变得非常小，成本不再改变（尽管应该改变）。我认为这是低精度的结果。但也许我可以使用您上面的建议来加速mp.math计算逆矩阵。 - user2153813

1

因此，我将这个想法与numpy中内置的float128支持相结合，得到了一个中间解决方案。在我的机器上，float128似乎支持约1e-19或1e-20的精度（比double多约1e3或1e4）。不幸的是，numpy.linalg不支持float128，因此我无法直接进行128位反演。但是，在64位中反演，然后转换为128位并使用如上所述的牛顿方法似乎可以得到一个很好的折衷方案。对于大矩阵（~100行），当计算norm(A*x - I)时，我得到1e-16（128位）而不是1e-13（64位答案）。这也比mpmath快得多。 - user2153813

好的，这听起来像是一个不错的妥协。:) 是的，128位FP在大多数硬件上仍然是模拟的（没有直接的硬件支持），但它比更通用的任意精度例程要快。 - Stuart

顺便说一句，我刚刚检查了一下，在Windows下，Numpy不支持float128。在Linux下，它是被支持的，但并不是真正的四倍精度，实际上只有80位。不过好消息是，这是一个FPU支持的格式，所以不需要模拟。从你观察到的情况来看，从64位浮点数到80位理论上将分辨率从约1e-16增加到约1e-19。 - Stuart

哦，好的，谢谢你告诉我。你能给我这个的参考资料吗？在numpy中，Float128似乎没有得到很好的记录。你有什么想法是否有更接近真正的128位支持的东西呢？ - user2153813

显示剩余4条评论

1

MATLAB多精度工具箱使用128位精度（Core i7 930）提供以下计时：

20x20 - 0.007秒

30x30 - 0.019秒

60x60 - 0.117秒

200x200 - 3.2秒

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

- Pavel Holoborodko

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Fredrik Johansson · Accepted Answer

遗憾的是，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矩阵相乘的速度提高16倍，比mpmath快。这种方法也不难直接编写矩阵求逆程序。当然，如果矩阵条目非常大或非常小，您需要先重新缩放它们。更可靠的解决办法是使用gmpy中的浮点类型编写自己的矩阵函数，这应该与现有方法基本相同。