I get, apparently, unacceptably high inaccuracies when calculating matrix inversions (linear system solution) in numpy.
- Is this a normal level of inaccuracy?
- How can I improve the accuracy of these calculations?
- Also, is there a way to more efficiently solve this system in numpy or scipy (
scipy.linalg.cho_solveseems promising, but doesn't do what I want)
The code below cholMcontains a 128 x 128 matrix. The matrix data is too large to include here, but it is located on pastebin: cholM.txt .
Also, the original vector ovecis randomly selected, so the ovecaccuracy varies for different ones , but in most cases the error still seems unacceptably high.
Change . Solving a system using a decomposition of singular values gives a significantly lower error than other methods.
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 )
source
share