Add SVD function specification

[kgryte]

This PR

• specifies the interface for performing singular value decomposition (SVD).
• is derived from comparing signatures across array libraries.

Notes

• NumPy and Torch order the returned tuple (u, s, v), while TF orders as (s, u, v). This proposal follows NumPy.

• NumPy et al provide a compute_uv keyword to indicate whether to return the left and right singular vectors along with the singular values. This proposal omits such a keyword, in favor of a separate svdvals specification (see gh-160). This is to promote consistency with eig and eigvals and minimize polymorphism within the set of linagl APIs.

• By default, NumPy computes full matrixes, while TF and Torch do not. This proposal follows NumPy and sets the full_matrices keyword argument default as True.

• MXNet does not currently support any keyword arguments.

• Dask currently supports a coerce_signs keyword argument to indicate whether or not to apply sign coercion to singular vectors in order to maintain deterministic results. It is alone in this regard.

• This proposal follows NumPy, Torch, MXNet, JAX, and TF in supporting stacks of matrices. CuPy and Dask do not currently support providing stacks.

• NumPy currently supports a hermitian keyword argument for speeding up computation; however, it is alone in doing this.

• TF (and LAPACK) returns both the left and right singular vectors in columns. NumPy returns the left singular vectors in rows (i.e., NumPy returns the adjoint). This proposal follows NumPy.

[leofang]

CuPy and Dask do not currently support providing stacks.

I checked that we have this capability (added in cupy/cupy#3247), so it's straightforward for us to support it:
https://github.com/cupy/cupy/blob/a6c75b901caa3be19ce8c4f717f2f780e457559d/cupy/cusolver.py#L71
It's just that for some reason we ended up staying with the QR-based SVD (gesvd, without the suffix j) from cuSOLVER, which does not have a batched version.

[rgommers]

This proposal follows TF (and LAPACK) in requiring that both the left and right singular vectors be returned in columns. NumPy returns the left singular vectors in rows.

tensorflow.experimental.numpy doesn't have an svd function yet, PyTorch matches NumPy, not sure about the rest (probably matches NumPy as well) - so this choice doesn't seem quite right. It's also hard to interpret the specification here. How about adding a note on this, mentioning how to reconstruct the input:

s, u, v = svd(x)
y = dot(u, dot(s, v))   # or does this need a transpose somewhere??
assert_allclose(x, y)

[rgommers]

There's an incompatibility between NumPy and PyTorch missed in number of values returned for compute_uv = False - see #95 (comment)

[leofang]

CuPy and Dask do not currently support providing stacks.

I checked that we have this capability (added in cupy/cupy#3247), so it's straightforward for us to support it:
https://github.com/cupy/cupy/blob/a6c75b901caa3be19ce8c4f717f2f780e457559d/cupy/cusolver.py#L71
It's just that for some reason we ended up staying with the QR-based SVD (gesvd, without the suffix j) from cuSOLVER, which does not have a batched version.

So it turns out not as straightforward, but is still possible on both CUDA and HIP, see cupy/cupy#4628.

btw, we could also note the support for complex numbers will land in the next version (like what we discussed in the eigensolver PR #113). The nice point of SVD is that u and v are real if the matrix is real, so unlike eigensolvers here it is not blocked by the lack of complex support.

[kgryte]

Updated this PR based on the above feedback and meeting discussions. Reordered the output values to follow NumPy (u,s,v) and specified that v be the adjoint (i.e., right singular vectors returned in rows).

[kgryte]

Updated this PR to return an array, rather than a tuple, when compute_uv is False.

[IvanYashchuk]

Why does compute_uv = True/False exists in NumPy and should NumPy be followed in this case? There is no compute_v = True/False for numpy.linalg.eig and numpy.linalg.eigh, there are separate functions for that: numpy.linalg.eigvals and numpy.linalg.eigvalsh.

Should this flag be dropped for svd and functionality replaced by a separate function svdvals, similarly to SciPy and Julia?

[rgommers]

Good question @IvanYashchuk. No/fewer boolean keywords would make the design better, that is probably how we'd design it if starting from scratch. I do see that scipy.linalg.svd has a compute_uv keyword. That's a bit unnecessary because it has svdvals too, it was probably added for compatibility with numpy.

I like the idea of dropping compute_uv, the comparison with eig/eigvals is nice.

[rgommers]

The signature svd(x, /, *, full_matrices=True) would be fully compatible with what existing libraries already have, so this may be a case where we can make the design more consistent without introducing extra issues/work for libraries that want to support this API standard in their current/main namespace.

[leofang]

Why does compute_uv = True/False exists in NumPy and should NumPy be followed in this case?

I think in reality it depends on how things are implemented at the low level. Taking CUDA as an example, this function cusolverDn<t>gesvd() is one of the routines implementing SVD, and if you take a closer look it provides both compute_uv and full_matrices options, which CuPy (and at least JAX I think) takes advantage of.

In fact, I think this interface dates way back to Lapack's <t>gesvd, and I assume this was why NumPy had this API design in the first place. So regardless how we design/split the API, under the hood they all call the same routine. So, it might lead to some code duplication if we split, though it's an implementation detail that I am not sure we should care here.

But, if we take into account decades of familiarity in numerical routines I don't think changing the API is a good idea, though I don't have strong opinion.

[kgryte]

I'm also in favor of @IvanYashchuk's proposal to have a separate API which only returns singular values.

Regarding whether an implementation reuses the same routine for both svd and svdvals I don't think should be a strong factor in how we specify the APIs. In general, APIs which always return the same object shapes (array or tuple, but not either) should be desired given that code is (a) easier to optimize (non-polymorphism) and (b) easier to reason about.

[kgryte]

I've updated this proposal (and OP) to no longer include the compute_uv keyword in favor a separate svdvals to support only returning the singular values (see gh-160).

[kgryte]

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