3

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

我正在达到机器精度的极限,所以已经开始尝试使用mpmath任意精度的矩阵求逆,但它非常慢,即使使用gmpy.

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

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

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

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

4

3 回答 3

2

不幸的是,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中的浮点类型编写自己的矩阵函数,这应该基本上一样快。

于 2013-03-11T06:42:58.490 回答
1

我假设双精度对于最终结果的精度来说不是问题,但对于某些矩阵,它会导致逆的中间结果出现问题。在这种情况下,让我们将普通 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
于 2013-03-10T16:21:26.417 回答
1

MATLAB 多精度工具箱使用 128 位精度(Core i7 930)提供以下时序:

20x20 - 0.007 秒

30x30 - 0.019 秒

60x60 - 0.117 秒

200x200 - 3.2 秒

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

于 2013-07-08T07:10:09.300 回答