CODEWITHSUNDEEP
numpyscipylinear-algebra
user1149913
user1149913

Reputation: 4523

Numpy inaccurate matrix inverse

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

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 )

Upvotes: 5

Views: 3652

Answers (2)

hpaulj
hpaulj

Reputation: 231385

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

We could find the solution vector using a matrix inverse: ... However, it is better to use the linalg.solve command which can be faster and more numerically stable

edit - from Steve Lord at MATLAB http://www.mathworks.com/matlabcentral/newsreader/view_thread/63130

Why are you inverting? If you're inverting to solve a system, don't -- generally you would want to use backslash instead.

However, for a system with a condition number around 1e17 (condition numbers must be greater than or equal to 1, so I assume that the 1e-17 figure in your post is the reciprocal condition number from RCOND) you're not going to get a very accurate result in any case.

Upvotes: -1

user1149913
user1149913

Reputation: 4523

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.

Upvotes: 2

Related Questions