Handling ndarrays with singular matrices in linalg functions #6027

richardotis · 2015-06-30T00:06:00Z

I am trying to solve several independent systems of equations at the same time using numpy.linalg.solve, e.g., a has shape (N, M, M). The challenge I'm running into is how to deal with the case when, for some values of N, the last two dimensions comprise a singular matrix.

Consider this example

import numpy as np
arr = [[[1, 0], [0, 1]], [[1, 0], [0, 1]], [[0, 0], [0, 0]]]
np.linalg.solve(arr, 1)

On numpy 1.9.2, this will raise a LinAlgError, but I still would like to know the solutions for all the full-rank matrices. (I would be okay with getting back NaNs for the singular cases.) I'm aware I can compute the singular values and use fancy indexing to slice the array, but in my algorithm I have to do several slicing/filtering steps and I would prefer not to lose the alignment of a with b -- my real-world a typically has 6-7 dimensions and all the indexing arrays make the code hard to follow.

My ideal solution seems to involve masked arrays since I end up doing several other filtering steps, but it appears that most (all?) functions in the linalg family ignore array masks.

Is there an alternative solution already available, or perhaps a suggestion on how I could contribute a solution that I could package as a PR?

The text was updated successfully, but these errors were encountered:

njsmith · 2015-06-30T00:44:09Z

I don't know of any really good workarounds, no.

There was some discussion of this issue when we first merged the ability to solve stacks of equations in a single call (which happened in #3220, a follow-up to #2954, and see also #3217). Basically, we all agreed that the LinAlgError behaviour is suboptimal for the reasons you suggest (inability to get out partial results), but that getting this right required some proper thought, and we didn't want to try and make up some hack in a few minutes while in the middle of getting the basic infrastructure in place. So we decided to stick with the crude LinAlgError approach at first, defer a proper fix until someone found time and attention to do it properly.

I posted some general thoughts on how this might work, and the questions we'd need to answer, here.

If I were going to spitball a solution off the top of my head, it'd be something like:

We add a new thread-local storage that stores "error flags" (in the sense of np.seterr); each error type gets assigned a single flag bit, and the ufunc code is adapted to add a way to set these flag bits and also check for them at the end of a computation. The idea here is that, right now, the only way for a ufunc loop to communicate an error condition to the outside world is through the floating point error flags that are built into the CPU. (This is how np.seterr etc. work right now.) This is awkward because it means that we are restricted to the IEEE754 error flags, and also because explicitly setting/clearing these flags requires cumbersome (and sometimes slow) FPU code. Of course we would still have to check the FPU flags in general, but having a second set of flags that we control makes things much more flexible.
Take advantage of this flexibility by adding a new linalg error type to np.seterr
Teach the functions in np.linalg to set this flag instead of directly raising an exception.

So the idea would be that by default we would continue to error (or at least warn) when singular matrices were passed to solve (i.e., the default setting for the linalg error handler would be error or warn), but for your code where you are expecting singular matrices and are prepared for them, you could write:

with np.errstate(linalg="ignore"):
    np.linalg.solve(stack_with_singular_matrices)

and silently get back matrices full of NaNs when appropriate.

That's just one idea, though; any real proposal would need a bit of thought into the questions posted in that link above, and to then get posted on the mailing list for wider review to make sure we didn't accidentally miss some important issue.

ovillellas · 2015-06-30T08:49:30Z

As far as I remember, the underlying gufunc used to implement solve does report any error in the underlying call to LAPACK by setting the results to NaNs. The linalg module function provides an "extobj" so that if that condition happens the exception is raised.

Maybe using seterr to 'ignore' for invalid will prevent the exception from being raised, with invalid entries having its results set to NaN. It could also be possible to call the underlying gufunc without extobj (and using the appropriate seterr as well).

Maybe adding a keyword to the linalg wrappers so that they don't raise exceptions and just return the results with NaNs for invalid entries would be useful (this could be implemented in all the linalg functions that are 'gufunc-like' and that can raise exceptions).

Another option would be to have the results available even if an exception is raised (maybe in the exception object). This would allow the user the following 3 options:

Ignore the exception by recovering the results with NaNs in invalid entries, maybe applying corrective code to the results.
Handle like an error in a way that is compatible with current code.
Use different codepaths for "all-correct" and "some-failed" results.

pv · 2015-06-30T09:22:56Z

IIRC, what was not decided is the precise API to be used for the error handling. . Raising exceptions is IMHO a bad idea in the first place. Small fp perturbations can make the LAPACK algorithms return with no error, but with garbage results --- the exceptions don't guard against this and are just a nuisance. Error tolerances / condition number limits could be used, but their values would need to be tunable by the user. Many other compsci software just emits warnings on potential numerical instability, returning best-efforts results. . Implementation-wise, controlling error state with keyword arguments in the linalg functions should be fairly simple to implement just by changing the code in linalg.py that raises the errors. Extending np.seterr with an additional linalg flag requires more extensive changes. . Internally, numpy.linalg abuses the 'invalid' flag for reporting linalg errors. seterr has currently no effect on raising errors; an explicit error object is always used.

ovillellas · 2015-06-30T10:33:41Z

In fact, when implementing the gufuncs, error handling was a bit of a caveat, but
I went the NaN-way as it was enough for my purposes (another caveat was not being
able to 'compute' output dimensions).

IMHO, a plain exception raised when any of the elements has an error condition is
not appropriate without extra information. At least it should be possible to have access
to the elements that didn't have problems when computing them, as well as knowing
which elements had errors -and ideally which errors-. That would require some tweaks
in the ufunc machinery to support error conditions. Another option would be to have an
extra output argument with per-element error codes that could be handled by the wrapper
function.

Note that implementing a "noexcept" keyword would just be a matter of not providing and extobj to
the underlying ufunc. That will provide the "errors as NaNs" semantics, respecting the seterr settings.
As long as the default is False, the old behaviour would be met.

pv · 2015-06-30T11:03:56Z

The linalg behavior probably shouldn't be controlled by the value of the `invalid` seterr flag, I think, because linalg errors aren't directly connected to the FP flags. noexpect=True then should pass in an extobj that specifies 'ignore'. . Also, I think it is not a good approach to go with an approach where "error" directly means "LAPACK routine returned error code" in all cases. What 'error' means may be fairly different for eigenvalue computation and matrix inversion. Although for the latter LAPACK returns an error if it encounters pivoting failure, the amount of precision loss is a continuum and a 'successful' return may well contain garbage.

glyg mentioned this issue Dec 16, 2015

Numerical instability in solving A\dot{X} + B = 0 bnoi/kt_simul#6

Closed

pv mentioned this issue Jan 27, 2018

Inverse of singular matrix. Should be error but in fact not. #10471

Closed

pmli mentioned this issue Feb 23, 2019

np.linalg.inv doesn't throw error on singular matrix inversion #13011

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Uh oh!

Handling ndarrays with singular matrices in linalg functions #6027

Handling ndarrays with singular matrices in linalg functions #6027

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Handling ndarrays with singular matrices in linalg functions #6027

Handling ndarrays with singular matrices in linalg functions #6027

Comments

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!