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 )

- user1149913

---

通常这表明您的矩阵病态。对于您提供的下三角矩阵，最大奇异值与最小奇异值之比（条件数）约为10^16。这绝对是一个问题。 - Craig J Copi

你可以尝试将cholM和invCholM相乘，然后减去单位矩阵，以说服自己条件数是否重要。 - aka.nice

这本质上是我在示例代码中所做的相同的事情。 - user1149913

条件数 根据定义 衡量了最大误差。因此，它无疑会给出一个关于朴素回代精度的可靠估计。你提供的矩阵的条件数约为1e16，而不是44.6。除了正常矩阵（cholM 明显不是）之外，特征值通常不能给出条件数。 - pv.

这里似乎存在一些误解。在数学世界中，条件数不会影响解决方案；但在计算机科学领域中，它们确实会影响解决方案。条件数为1e16绝对是个坏消息。首先确保您使用的是双精度。如果这样做没有帮助，那么解决相同系统的不同方法可能对数值误差的敏感性有很大的差异。我看到的那个回代也完全属于“不要在糟糕条件下尝试这样做”的范畴。 - Eelco Hoogendoorn

Answer 1

如@Craig J Copi和@pv所指出的，cholM矩阵的条件数高达10^16，这意味着为了获得更高的逆矩阵精度，需要使用更高的数值精度。

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

- user1149913

Answer 2

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得到的倒数条件数），您无论如何都不能得到非常准确的结果。

- hpaulj