python - Efficient & pythonic check for singular matrix

31
votes

Working on some matrix algebra here. Sometimes I need to invert a matrix that may be singular or ill-conditioned. I understand it is pythonic to simply do this:

try:
    i = linalg.inv(x)
except LinAlgErr as err:
    #handle it

but am not sure how efficient that is. Wouldn't this be better?

if linalg.cond(x) < 1/sys.float_info.epsilon:
    i = linalg.inv(x)
else:
    #handle it

Does numpy.linalg simply perform up front the test I proscribed?

4
pythonnumpylinear-algebra
Trying something and reacting to error conditions after they occur is usually the Pythonic way to do things, but the Pythonic way is not always the most efficient way. It's presumed that if you care about doing things "Pythonically," efficiency is a secondary priority. - David Z
Say not that efficiency is a secondary priority; say instead that I want to perform bivariate optimization: pythonicity + efficiency (hence the post title). - Dr. Andrew
In numerical computing it is usually considered bad practice to explicitly calculate the inverse. In most cases it is much better to calculate the LU decomposition with scipy.linalg.lu_factor, then later you can solve it quickly for many vectors using scipy.linalg.lu_solve. - DaveP
Correct, but in some cases, the inverse is actually required. I'm not inverting simply to use the result to solve a system of equations. I need the result. - Dr. Andrew
For others with this issue, as myself, the LLinAlgError need to be loaded from numpy, as numpy.linalg.linalg.LinAlgError. - Austin Downey

4 Answers

21
votes

So based on the inputs here, I'm marking my original code block with the explicit test as the solution:

if linalg.cond(x) < 1/sys.float_info.epsilon:
    i = linalg.inv(x)
else:
    #handle it

Surprisingly, the numpy.linalg.inv function doesn't perform this test. I checked the code and found it goes through all it's machinations, then just calls the lapack routine - seems quite inefficient. Also, I would 2nd a point made by DaveP: that the inverse of a matrix should not be computed unless it's explicitly needed.

15
votes

Your first solution catches the case where your matrix is so singular that numpy cannot cope at all - potentially quite an extreme case. Your second solution is better, because it catches the case where numpy gives an answer, but that answer is potentially corrupted by rounding error - this seems much more sensible.

If you are trying to invert ill-conditioned matrices, then you should consider using singular value decomposition. If used carefully, it can give you a sensible answer where other routines fail.

If you don't want SVD, then see also my comment about using lu_factor instead of inv.

5
votes

You should compute the condition number of the matrix to see if it is invertible.

import numpy.linalg

if numpy.isfinite(numpy.linalg.cond(A)):
    B = numpy.linalg.inv(A)
else:
    # handle it
0
votes

Why not just check if the determinant is nonzero?

det = numpy.linalg.det(A)
if det != 0:
   #proceed