[Numpy-discussion] Inversion of near singular matrices.
Sturla Molden
sturla at molden.no
Sun Jan 30 10:15:34 EST 2011
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
More information about the NumPy-Discussion
mailing list