-
-
Notifications
You must be signed in to change notification settings - Fork 12.2k
Merge generalized linear algebra routines into numpy.linalg #3217
Description
The features in gh-2954 should be brought into numpy.linalg.
A suggested plan for action:
- Combine the *_lite libraries.
- Merge in Linear Algebra routines as generalized ufuncs #2954, without exposing the features in a public API.
- Change (Add generalized ufunc linalg functions and make numpy.linalg use them #3220) the implementation of routines in numpy.linalg to use umath_linalg while preserving full backward compatibility, but allow "stacked" input.
- Add the new routines from
_gufunc_linalg.pyinto numpy.linalg (keep eye on Scipy function signatures, if applicable) - Figure out how to deal with error handling in numpy.linalg and implement any changes.
- Figure out how to allow numpy.linalg routines compute single-precision in single precision (rather than using double precision LAPACK)
- Extend umath_linalg to cover all functions (if possible?) in numpy.linalg, and drop the old lapack_lite wrapper.
Open questions:
A. Error handling. Raising LinAlgErrors is not necessarily the right thing to do especially in the stacked case, so this should be made user controllable.
B. The numpy.linalg routines compute results for single-precision by using double precision internally. This should perhaps be made tunable?
C. Which of the generalized ufunc arguments should be exposed in the numpy.linalg API?
D. The umath_linalg error handling goes through by mis-using the invalid FP flag. Should ufunc error handling be extended here? (scipy.special would also have an use case)
Possible solutions to questions:
A. Add keyword argument raise_errors to all functions (inv, solve, eig, ...) that used to raise a LinAlgError on LAPACK errors. Perhaps also add global state, e.g. np.linalg.seterr for tuning the default value --- however, global state here may have unexpected effects in 3rd party libraries.
B. Changing the default can cause changes in results, and may make some algorithms not converge?
C. Dtype, perhaps, and maybe also axis control arguments?