Efficient & pythonic check for singular matrix

2019-01-17 11:28发布

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条回答
老娘就宠你
2楼-- · 2019-01-17 11:53

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
查看更多
爷的心禁止访问
3楼-- · 2019-01-17 12:07

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.

查看更多
时光不老,我们不散
4楼-- · 2019-01-17 12:07

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.

查看更多
我欲成王,谁敢阻挡
5楼-- · 2019-01-17 12:12

Why not just check if the determinant is nonzero?

det = numpy.linalg.det(A)
if det != 0:
   #proceed
查看更多
登录 后发表回答