首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >mpmath矩阵求逆的替代方案或加速比

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

Stack Overflow用户
提问于 2013-03-10 21:21:45
回答 3查看 2K关注 0票数 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机器上使用python3mpmath 0.17 (不确定gmpy版本)。我不知道为什么mpmath会慢这么多,但是有没有开源库可以达到mathematica的速度(即使是1/2的速度也不错)?

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

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

EN

回答 3

Stack Overflow用户

发布于 2013-03-11 14:42:58

不幸的是,mpmath中的Linera代数相当慢。有很多库可以更好地解决这个问题(例如Sage)。这就是说,作为Stuart建议的后续,在Python中使用定点算法在不安装任何库的情况下进行相当快速的高精度矩阵乘法是相当容易的。下面是一个使用mpmath矩阵作为输入和输出的版本:

代码语言:javascript
复制
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中的浮点类型编写您自己的矩阵函数,这在本质上应该是同样快的。

票数 2
EN

Stack Overflow用户

发布于 2013-03-11 00:21:26

我假设双精度对于最终结果的精度不是问题,但对于某些矩阵,它会导致逆矩阵的中间结果出现问题。在这种情况下,让我们将正常numpy (双精度)逆的结果仅视为良好的近似值,然后将其作为牛顿方法的几次迭代的起点来求解逆。

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

代码语言:javascript
复制
X = X*(2I - AX)

对于大型矩阵,与求逆相比,计算上述几次迭代的努力几乎微不足道,而且它可以极大地提高最终结果的准确性。试一试。

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

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

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

代码语言:javascript
复制
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
票数 1
EN

Stack Overflow用户

发布于 2013-07-08 15:10:09

Multiprecision toolbox for MATLAB使用128位精度(核心i7 930)提供以下计时:

20x20 - 0.007秒

30x30 - 0.019秒

60x60 - 0.117秒

200x200 - 3.2秒

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

票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/15322686

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档