Add specification for computing the pseudo-inverse (linalg: pinv)

[kgryte]

This PR

• specifies the interface for computing the (Moore-Penrose) pseudo-inverse.
• is derived from comparing signatures across array libraries.

Notes

• Following Torch, MXNet, TF, NumPy, and JAX, this proposal allows for providing a stack of square matrices. CuPy does not currently support providing stacks.

• Dask does not provide an API for computing the pseudo-inverse.

• TF supports a validate_args argument for embedding additional validations within its computational graph.

• NumPy, MXNet, and Torch (latest master) supporting providing a hermitian keyword argument to indicate that more efficient computation methods be used. This PR omits this keyword, as more of an implementation detail, than a generalizable API.

• NumPy, MXNet, CuPy, and Torch set the default rcond value to 1e-15, while JAX and TF compute a default value based on the machine epsilon associated with the input array data type and the number of rows/cols. This PR follows JAX and TF in computing the default value (as 1e-15 does not make sense for non-float64 input, such as float32 or bfloat16) and requiring that rcond be a broadcast compatible array (or a float).

• This proposal renames the rcond keyword argument to rtol in order to unify keyword arguments for pinv, lstsq, and matrix_rank which all support specifying relative tolerances. The default value is also the same across these APIs.

• Question: should this return a namedtuple to allow for a variable number of returns (see API for variable number of returns in linalg #95)? SciPy, e.g., does support returning multiple values (the matrix B along with the effective rank of the result). In theory, other info could be returned, such as error info, but not clear whether this is enough of a forward-looking concern.

  • Answer: no.

[rgommers]

Question: should this return a tuple to allow for a variable number of returns (see #95)? SciPy, e.g., does support returning multiple values (the matrix B along with the effective rank of the result). In theory, other info could be returned, such as error info, but not clear whether this is enough of a forward-looking concern.

I'd say no. We should avoid variable number of returns as much as possible. Returning a tuple doesn't help at all; changing the tuple length would still be a serious backwards compat break.

I'll comment on gh-95, after looking at some of these PRs that's worth reconsidering.

[rgommers]

NumPy, MXNet, CuPy, and Torch set the default rcond value to 1e-15, while JAX and TF compute a default value based on the machine epsilon associated with the input array data type and the number of rows/cols. This PR follows JAX and TF in computing the default value (as 1e-15 does not make sense for non-float64 input, such as float32 or bfloat16) and requiring that rcond be a broadcast compatible array.

The default value here is given as 10.0 * max(M, N) * eps. I see a couple of potential issues:

• The factor 10.0 is missing from the lstsq proposal in Add specification for returning the least-squares solution to a linear matrix equation (linalg: lstsq) #119.
• This is a lot larger than 1e-15 for large inputs. Intuitively, rcond should not scale linearly with the size of the largest dimension - perhaps more like the square root of it.

[kgryte]

Re: 10.0 factor. Yeah, I am not sure the reasoning behind JAX and TF's use of the factor. That factor is absent from lstsq because the factor is not used by NumPy and was not clear to me whether the rcond defaults for pinv and lstsq should match.

[kgryte]

Re: namedtuple. Another alternative is to simply return a dictionary.

[kgryte]

Renamed rcond to tol in order to unify similar keyword arguments across pinv, lstsq, and matrix_rank. Removed 10.0 scaling factor; the default tolerances are now computed consistently across each of the three aforementioned APIs.

[leofang]

Just a quick comment 🙂

[kgryte]

Renamed tol to rtol to more explicitly indicate relative tolerance and pave the way for future specification evolution (e.g., atol).

[leofang]

btw CuPy will support batched pinv in the upcoming v9.0 (cupy/cupy#4686).

[kgryte]

Thanks, @leofang, for the review! This is ready for merge...