Numpy不准确的矩阵逆答案

【问题标题】：Numpy inaccurate matrix inverseNumpy不准确的矩阵逆
【发布时间】：2014-03-17 09:09:17
【问题描述】：

在 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 )

【问题讨论】：

这通常表明您的矩阵是病态的。对于您提供的下三角矩阵，最大奇异值与最小奇异值（条件数）的比率似乎约为 10 ^ 16。这绝对是个问题。
为什么条件数在这里很重要？我正在取逆，所以通常情况下，如果矩阵接近奇异，我会担心。 cholM 矩阵是三角形的，所以对角线上非常小的值肯定会影响解，但这种情况不会发生。
您可以尝试将 cholM 和 invCholM 相乘，然后减去单位矩阵以说服自己条件数是否重要
条件号by definition 衡量您可以获得的最大误差。因此，它无疑将忠实地估计朴素反向替换的准确性。您提供的矩阵的条件数是 ~1e16，而不是 44.6。特征值通常不会给你条件数，除了正态矩阵（cholM 显然不是）。
这里似乎有误会。条件数不会影响数学世界中的解决方案；它们确实影响了计算机科学领域的解决方案。条件编号 1e16 绝对是个坏消息。确保您对初学者使用双精度。如果这没有帮助；求解同一系统的不同方法对数值误差的敏感性可能大不相同。我看到的这种反向替代也完全属于“不要在条件不好的情况下尝试这个”阵营。

标签： numpy scipy linear-algebra

【解决方案1】：

如@craig j copi和@pv所指出的，condition number cholM矩阵很高，大约为10 ^ 16，表明为了实现更高的准确性，可能需要更大的数值精度。

条件数可以通过最大奇异值与最小奇异值的比率来确定。在这种情况下，该比率与特征值的比例不同。

【讨论】：

【解决方案2】：

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

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

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

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

但是，对于条件编号在 1e17 左右的系统（条件编号必须大于或等于 1，所以我假设 1e-17 中的数字你的帖子是来自 RCOND 的倒数条件数）你不会在任何情况下都能得到非常准确的结果。

【讨论】：