Numpy inaccurate matrix inverse

Question

I have been getting seemingly unacceptably high inaccuracies when computing matrix inverses (solving a linear system) in numpy.

• Is this a normal level of inaccuracy?
• How can I improve the accuracy of this computation?
• Also, is there a way to solve this system more efficiently in numpy or scipy (scipy.linalg.cho_solve seemed promising but does not do what I want)?

In the code below, cholM is a 128 x 128 matrix. The matrix data is too large to include here but is located on pastebin: cholM.txt.

Also, the original vector, ovec, is being randomly selected, so for different ovec's the accuracy varies, but, for most cases, the error still seems unacceptably high.

Edit Solving the system using the singular value decomposition produces significantly lower error than the 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 )

Answer 1

As noted by @Craig J Copi and @pv, the condition number of the cholM matrix is high, around 10^16, indicating that to achieve higher accuracy in the inverse, much greater numerical precision may be required.

Condition number can be determined by the ratio of maximum singular value to minimum singular value. In this instance, this ratio is not the same as the ratio of eigenvalues.