37

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?

Praveen
  • 6,872
  • 3
  • 43
  • 62
Dr. Andrew
  • 2,532
  • 3
  • 26
  • 42
  • 1
    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 Nov 06 '12 at 11:03
  • 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 Nov 06 '12 at 11:30
  • 3
    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 Nov 07 '12 at 07:03
  • 1
    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 Nov 07 '12 at 13:19
  • 2
    For others with this issue, as myself, the LLinAlgError need to be loaded from numpy, as `numpy.linalg.linalg.LinAlgError`. – Austin Downey May 26 '17 at 14:10
  • whats wrong with your suggested solution? – Charlie Parker Oct 21 '17 at 01:48
  • Also see https://github.com/numpy/numpy/issues/2074 for cases where the first solution will fail due to floating point inaccuracy. A very simple example: `vec = atleast_2d([4,9]); linalg.inv(vec.T @ vec)` does not throw an exception but obviously yields a completely nonsensical result. – Eike P. Sep 24 '18 at 15:04

4 Answers4

27

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.

Dr. Andrew
  • 2,532
  • 3
  • 26
  • 42
  • Instead of `sys.float_info.epsilon` I prefer `numpy.spacing` since it can be given 1.0 as float16, float32, float64, ..., whatever type `x` is in. – Ahmed Fasih Jan 10 '14 at 02:04
  • 4
    You should use the epsilon for the dtype of the array itself, i.e. `np.finfo(x.dtype).eps`, which may differ from `sys.float_info.epsilon` – ali_m Jul 05 '15 at 11:54
  • From performance point of view first way (catching exception) is much faster for big matrices: `In [56]: paA = np.identity(10000) In [57]: %timeit np.linalg.cond(paA) 7min 20s ± 11.8 s per loop (mean ± std. dev. of 7 runs, 1 loop each) In [58]: %timeit np.linalg.det(paA) 11.2 s ± 817 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) In [59]: %timeit np.linalg.inv(paA) 40.4 s ± 1.14 s per loop (mean ± std. dev. of 7 runs, 1 loop each)` – Prokhozhii Dec 04 '18 at 10:05
20

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.

DaveP
  • 6,952
  • 1
  • 24
  • 37
  • The 1st is what I guess (and seems confirmed) is most pythonic, but I wouldn't call it "my solution" :-). After catching the ill-conditionedness, I already know multiple regularization methods I can use to fix it. – Dr. Andrew Nov 07 '12 at 13:22
  • I guess my main message is that before you worry about being "most pythonic", you first make sure that you are being most numerically correct. The second solution you propose is better from that point of view. – DaveP Nov 07 '12 at 21:50
  • Completely agree with you there DaveP. – Dr. Andrew Nov 08 '12 at 09:55
6

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
Nicolas Barbey
  • 6,639
  • 4
  • 28
  • 34
  • 3
    The matrix will be ill-conditioned when the condition number approaches 1/epsilon. So the second solution which was proposed in the question is better than your solution, which only catches an extremely large condition number. – DaveP Nov 07 '12 at 07:06
0

Check if the determinant is nonzero (but be careful, see comments below):

det = numpy.linalg.det(A)
if not numpy.isclose(det, 0):
   #proceed
Bjørn Kjos-Hanssen
  • 764
  • 7
  • 21
  • 7
    This would work if computers could represent floats *exactly*. However they can't, so you will most of the time end up with an answer like `det = 3.123e-14` and you never can be sure whether it's 0 or just a small number. But of course, if you could calculate exact numbers, checking for a nonzero determinant is the way to go. – tamasgal Jan 09 '17 at 08:05
  • @tamasgal: ... as some people have noticed empirically. ;) https://github.com/numpy/numpy/issues/2074 – Eike P. Sep 24 '18 at 14:58
  • numpy computes the determinant through approximation. for matrix `A`, numpy uses LU decomposition and then takes the product of the diagonals. in other words, numpy doesnt compute the determinant analytically but approximates it. – mynameisvinn Jun 04 '20 at 00:46
  • 1
    [`np.isclose()`](https://numpy.org/doc/stable/reference/generated/numpy.isclose.html) might be a potential workaround for using the determinant – ijuneja Jan 13 '22 at 21:35