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

Question

我正在用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的速度（即使是它的一半速度也不错）？

我并不需要任意精度，128位的精度应该就足够了。我也不明白为什么mpmath会慢得如此离谱。它一定是使用了非常不同的矩阵求逆算法。具体来说，我使用的是M**-1。

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

Answer 1

MATLAB的多精度工具箱提供了以下使用128位精度的时间数据（在Core i7 930处理器上测试）：

20x20 - 0.007秒

30x30 - 0.019秒

60x60 - 0.117秒

200x200 - 3.2秒

需要注意的是，这些数字在现代CPU上会更低。

Answer 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

Answer 3

在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中的浮点类型来自己编写矩阵函数，这样的速度应该也差不多。