[Numpy-discussion] Inversion of near singular matrices.

[Sturla Molden]

Den 30.01.2011 07:28, skrev Algis Kabaila:
> Why not simply numply.linalg.cond? This gives the condition
> number directly (and presumably performs the inspection of
> sv's). Or do you think that sv's give more useful information?

You can use the singular value decomposition to invert the matrix, solve 
linear systems and solve least squares problems. Looking at the topic 
you don't just want to compute condition numbers, but invert the 
ill-conditioned (nearly singular) matrix.

Say you want to do:

    invA =  np.linalg.inv(a)

With matrix a nearly singular, you can proceed like this:

    U, s, Vh = np.linalg.svd(a, full_matrices=False)

Edit small singular values. This will add a small bias but reduce 
rounding error:

    singular = s < threshold
    invS = 1/s
    invS[singular] = 0 # truncate inf to 0 actually helps...

Et voilá:

    invA = np.dot(np.dot(U,np.diag(invS)),Vh)

I hope this helps :)

There is a chapter on SVD in "Numerical Receipes" and almost any linear 
algebra textbook.

Just to verify:

 >>> a = np.diag([1,2,3])
 >>> np.linalg.inv(a)
array([[ 1.        ,  0.        ,  0.        ],
        [ 0.        ,  0.5       ,  0.        ],
        [ 0.        ,  0.        ,  0.33333333]])
 >>> U, s, Vh = np.linalg.svd(a, full_matrices=False)
 >>> np.dot(np.dot(U,np.diag(1/s)),Vh)
array([[ 1.        ,  0.        ,  0.        ],
        [ 0.        ,  0.5       ,  0.        ],
        [ 0.        ,  0.        ,  0.33333333]])


Sturla