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

Question

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

我已经达到了机器精度的极限，所以已经开始尝试使用mpmath来进行任意精度的矩阵求逆，但速度非常慢，甚至使用gmpy。在

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

这是在archlinux机器上使用python3和{}（不确定gmpy版本）。我不知道为什么mpmath要慢得多，但是有没有开源库可以达到mathematica为此所管理的速度（即使是速度的1/2也不错）？在

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

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

Answer 1

不幸的是，mpmath中的Linera代数相当慢。有许多库可以更好地解决这个问题（例如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中的浮点类型编写自己的矩阵函数，这基本上应该同样快。在

Answer 2

我假设双精度不是最终结果精度的问题，而是某些矩阵的问题，它在求逆的中间结果中引起了问题。在这种情况下，让我们把一个普通的numpy（双精度）逆的结果仅仅当作一个好的近似值，然后用它作为牛顿方法的几个迭代的起点来求解逆。在

设A是我们要求逆的矩阵，X是我们对逆矩阵的估计。牛顿法的迭代简单地包括：

X = X*(2I - AX)

对于大型矩阵，与求逆相比，计算上述迭代次数的工作量几乎是微不足道的，它可以大大提高最终结果的精度。试试这个。在

顺便说一下，I是上述等式中的单位矩阵。在

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

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

^{pr2}$

Answer 3

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

20x20-0.007秒

30x30-0.019秒

60x60-0.117秒

200x200-3.2秒

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