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