[Numpy-discussion] Inversion of near singular matrices.

Bruce Southey bsouthey at gmail.com
Sun Jan 30 10:40:10 EST 2011


On Sun, Jan 30, 2011 at 9:15 AM, Sturla Molden <sturla at molden.no> wrote:
> 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
>
>
If the matrix is not full rank then it is not invertible
(http://en.wikipedia.org/wiki/Matrix_inverse) and (matrix) rank
(http://en.wikipedia.org/wiki/Matrix_rank) can be computed from the
above code. You do have to beware that you can get a generalized
inverse (numpy.linalg provides pinv) when your system has an infinite
number of solutions. (A generalized inverse is not the problem, the
problem is when you expect a unique solution and do not get it.)

Bruce



More information about the NumPy-Discussion mailing list