alg keyword for svd and svd!

[carstenbauer]

This PR adds an alg keyword to svd and svd! which allows the user to switch between a divide-and-conquer algorithm (default, LAPACK.gesdd!) and a simple algorithm (LAPACK.gesvd!).

In the future, Jacobi methods could also be added. (#31080)

See #31023 and https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603 for related discussions.

Closes #31023.

[simonbyrne]

Can you add some tests?

I'm not sure I like alg as a keyword: maybe use the full algorithm (LAPACK calls them "drivers", but that is far from obvious).

[ararslan]

FWIW sort also uses alg.

[carstenbauer]

Can you add some tests?

Sure. What should I test? Only that the keyword argument exists or also that the results differ (maybe for something as in #31023)?

BTW, I can't find tests for the full keyword argument.

[StefanKarpinski]

Only that the keyword argument exists or also that the results differ (maybe for something as in #31023)?

I would test first that the keywords are accepted and that the results look like correct SVDs, i.e.:

1. that U*S*V' ≈ X
2. that U and V are approximately orthogonormal
3. that the singular values are all non-negative
4. if there are particular characteristics that each algorithm should have, test those

[carstenbauer]

This is being discussed on Slack. @andreasnoack pointed out the he made a similar change in CuArrays.jl recently.

Maybe it's better to call the keyword method as people seem to prefer it and it would be consistent with aboves change.

The remaining questions seem to be:

• How verbose do we want to be, e.g. :divide vs :divideandconquer?
• Do we want to use Symbols, Enums, or even types for the kw argument?

FWIW, I'd say verbosity might be ok here. It makes the code more self-explanatory and, given that this feature was absent for this long, people will only use this option very rarely.

[StefanKarpinski]

Maybe it's better to call the keyword method as people seem to prefer it and it would be consistent with aboves change.

Possibly but keep in mind that "method" already has a fairly well-established meaning in programming languages which does not really agree with this meaning. It's fine to call these methods informally but "algorithm" seems both more precise and less confusing in this context.

[simonbyrne]

Someone pointed out on Slack (I forgot who), that we do use alg in sort, so there is precedence.

I'm not hugely keen on using symbols as enums though.

[ararslan]

I said that in a comment in this PR 5 days ago. 😐

[andreasnoack]

I think we should consider using a struct for this. Then the svd can be extended easily by the JacobiSVD package.

[carstenbauer]

Alright, I updated (and rebased) this one. I'm unsure about naming. I introduced

abstract type SVDAlgorithm end
struct DivideAndConquer <: SVDAlgorithm end
struct Simple <: SVDAlgorithm end

but in particular Simple seems very general and vague. Any suggestions?

[andreasnoack]

Some of these types might make sense across more than a single factorization e.g. there is also a divide and conquer version of the symmetric tridiagonal eigensolver so either we should just have a shared Algorithm type in LinearAlgebra or we'd have to put SVD in the struct names. I'd prefer the former.

[StefanKarpinski]

I'm also not too keen on the name "Simple". That doesn't really narrow it down enough. Does this simple algorithm not have a proper name? Any linalg folks know what this is called?

[simonbyrne]

Golub and Van Loan just call it "The SVD Algorithm".

From this paper on the history of the SVD

The singular value decomposition was introduced into numerical analysis by [Golub and Kahan [23, 1965], who proposed a computational algorithm. However, it was Golub [24, 1970] who gave the algorithm that has been the workhorse of the past two decades.

The referenced papers are:

• Golub and Kahan (1965)
• Golub and Reinsch (1970)

"Golub-Kahan" usually refers to the bidiagonalization step. So I guess we could call it GolubReinsch? However as @andreasnoack points out, SVD algorithms are in some sense equivalent to symmetric eigenvalue algorithms: if I understand it correctly, the equivalent is the QR algorithm, which as the wikipedia article points out, is a terrible name.

Maybe @higham might be able to help here.

[carstenbauer]

Alright, would be great if someone could review this. AFAICT, the two run failures are unrelated.

[carstenbauer]

bump :)

[carstenbauer]

Would be great if this would make it into 1.3. Could someone review this before feature freeze?

[tpapp]

FWIW, I have reviewed, but note that I am not a core dev so I am not sure how much it helps. I just really like the change and would also love to see it in 1.3.

[StefanKarpinski]

@andreasnoack, could you possibly find the time to review this?

[andreasnoack]

It mostly looks ready to go. I've added a few comments and suggestions.

[carstenbauer]

Implemented the comments by @andreasnoack. Tests pass locally.

[andreasnoack]

@fredrikekre I've forgotten how we add the compat annotations. Could you help us out here?

[carstenbauer]

Until when does this have to be merged to be in 1.3?

[ararslan]

I've forgotten how we add the compat annotations

In the docstring for the function, add

!!! compat "Julia 1.3"
    The `method` keyword argument requires Julia 1.3 or later.

[carstenbauer]

@ararslan the keyword argument is named alg

I added the compat annotation and fixed the NEWS.md conflict.

[Jutho]

Great addition; I had this in some local packages. 1.3 will be a great release.