Numpy 不准确的矩阵逆

Question

在 numpy 中计算矩阵逆（求解线性系统）时，我似乎出现了令人无法接受的高错误。

• 这是正常的不准确程度吗？
• 如何提高此计算的准确性？
• 另外，有没有一种方法可以在 numpy 或 scipy 中更有效地解决这个系统（scipy.linalg.cho_solve看起来很有希望，但没有达到我想要的效果）？

在下面的代码中，cholM是一个 128 x 128 矩阵。矩阵数据太大，无法包含在此处，但位于 Pastebin：cholM.txt上。

此外，原始向量ovec是随机选择的，因此对于不同的ovec，精度会有所不同，但在大多数情况下，误差似乎仍然高得令人无法接受。

编辑使用奇异值分解求解系统所产生的误差明显低于其他方法。

import numpy.random as rnd
import numpy.linalg as lin
import numpy as np

cholM=np.loadtxt('cholM.txt')

dims=len(cholM)
print 'Dimensions',dims

ovec=rnd.normal(size=dims)
rvec=np.dot(cholM.T,ovec)
invCholM=lin.inv(cholM.T)

svec=np.dot(invCholM,rvec)
svec1=lin.solve(cholM.T,rvec)

def back_substitute(M,v):
    r=np.zeros(len(v))
    k=len(v)-1
    r[k]=v[k]/M[k,k]
    for k in xrange(len(v)-2,-1,-1):
        r[k]=(v[k]-np.dot(M[k,k+1:],r[k+1:]))/M[k,k]

    return r

svec2=back_substitute(cholM.T,rvec)

u,s,v=lin.svd(cholM)
svec3=np.dot(u,np.dot(np.diag(1./s),np.dot(v,rvec)))

for k in xrange(dims):
    print '%20.3f%20.3f%20.3f%20.3f'%(ovec[k]-svec[k],ovec[k]-svec1[k],ovec[k]-svec2[k],ovec[k]-svec3[k])

assert np.all( np.abs(ovec-svec)<1e-5 ) 
assert np.all( np.abs(ovec-svec1)<1e-5 )

Answer 1

http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html

我们可以使用矩阵逆来找到解向量：...但是，最好使用 linalg.solve 命令，它可以更快且数值更稳定

编辑 - 来自 MATLAB 的 Steve Lord http://www.mathworks.com/matlabcentral/newsreader/view_thread/63130

你为什么要反转？如果你要通过反斜杠来求解一个系统，那么不要——通常你会想使用反斜杠。

但是，对于条件数约为 1e17 的系统（条件数必须大于或等于 1，因此我假设您帖子中的 1e-17 数字是 RCOND 的条件数的倒数），您将不会得到无论如何，结果都非常准确。