Interface discussion

[Jutho]

@lkdvos , maybe we can discuss the interface here.

For the family of methods like qr, eigen and svd (possibly trunced) there needs to be support for:

• backends (LAPACK, CUSOLVER, KrylovKit, ...)

• algorithms supported by these backends (QRIteration, DivideAndConquer, RobustRepresentations, Jacobi,...)

• possible options/keywords to configure those algorithms

• possible keywords to control truncation? specific to svd?

Then there is something like polar? Do we want a specific method for polar decomposition. That could be computed via SVD (and thus need the previous levels of control), but there are also other packages (=backends?) that can do this without SVD.

If we also define methods like leftorth and friends here (this could also be restricted to TensorKit), they can be computed using e.g. either svd or qr. I don't know if there are other relevant strategies? So the question is how to convey that information? A strategy or kind arg/kwarg?

Something similar exists for matrix functions, like exponential, which can use things like Pade approximations, Schur-Parlett, diagonalisation in the case of hermitian matrices, ….

[Jutho]

I also forgot GenericLinearAlgebra.jl as backend. Am I correct that what you were proposing, @lkdvos , is to combine algorithms and backends in one name, something like LAPACK_QRIteration, … . That reminds me of what is done in LinearSolve.jl: https://docs.sciml.ai/LinearSolve/stable/solvers/solvers/

[lkdvos]

I think that indeed is very similar to what I had in mind. To me, the difference between backend and algorithm seems rather minimal, so I would just merge the two concepts. It almost sounds like what you describe as backend is simply just a family of algorithms that are found in the same package, but then it might as well be a type union or abstract type, if this is even needed.
I am very intrigued by the LinearSolve link, in the sense that it almost looks like what we need, except they don't actually expose the algorithms at the level of the factorizations. Scrolling through the available algorithms, it seems like there as well, there was demand for FastLAPACK - in-place results, which is interesting.
I quite like that the full is just a field in their SVD algorithm as well, that seems quite reasonable.

Just as an aside, this is what I had laying around locally, which was from before I looked at your most recent changes: https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/tree/ld-interfaceproposal
It's definitely not finished, nor am I perfectly happy with the design, but it does have some of the aspects that I was thinking about:

• Having a struct to bundle the output, whose constructor simultaneously serves as "pre-allocating" the results
• Having structs to define the algorithms
• Having an interface consisting of svd(A; kwargs...), svd(A, alg), svd!(A, [F]; kwargs...) and svd!(A, F, alg), somewhat reminiscent of KrylovKit, with a default algorithm selector (which is more or less equivalent to the default_backend, which I had somehow missed)
• After talking to you: it is unclear how factorization structs would interact with ChainRules, as this sounds like a mess, so before we can decide on that, this needs to be checked.

Something that I had been thinking about, specifically for svd, is whether or not it would make sense to have something like:

struct TruncatedSVD{A,T}
    svd_alg::A
    trunc_alg::T
end

and just treat truncated svd on the exact same footing as svd. Of course, depending on the algorithm being Krylov-based or LAPACK, the implementation might vary quite a bit, but at least the interface would be the same.

[lkdvos]

Something that I just thought about: If the backends really are "overarching themes" for a couple algorithms, we could just put them in different modules, which would lead to LAPACK.QRIteration, ...

[Jutho]

One problem with all of this is that all algorithm/backend types need to be defined in the main module, as extensions cannot create new types. This seems to be what LinearSolve.jl also does. I have not always been fan of some of the design decisions in the SciML ecosystem (and I've learned recently I am not the only one), but for the case of LinearSolve.jl I don't know yet. One certainty though, SVDFactorisation is just linguistic nonsense! 😄

Regarding the output structs, I am not necessarily opposed to this, but I am not sure if I see the great benefit? The one example of the two-step process F = svd(...); truncate!(F) only makes sense for the dense truncation, and would not work for the truncation based on Krylov or randomized methods. The only case where I think I have ever used the base factorisation types in a way other than directly destructuring them, is for inverting a matrix by first factorising it (and choosing the factorisation myself).

[Jutho]

I am confused by the discussion page, you can write a reply in a message, or below as a new comment. That is too many options.

Anyway, some more thoughts:

• regarding the output structs: what if we for now define the functions as svd_compact!(A, (U, S, Vᴴ), backend; ...)? That way, if we go with the structs, you can insert the struct as well and things would still work.
• regarding the argument order: I would keep A the first argument for various reasons:
  • it is easier to insert the default values for the output arguments
  • it is often also modified or destroyed
  • it might at some point be useful to be able to insert implicit operators in the form of functions
• regarding algorithm types, I have no strong preference, but I think admitting keyword arguments as part of the interface to further finetune certain algorithms is necessary. Whether these keyword arguments then become part of the algorithm structure (maybe via a general <:NamedTuple field) or are passed along separately is both ok for me; I haven't given it much thought
• I have a mild dislike for making the shape of the output depend on a keyword like full, as it is somewhat similar to having the actual number of output values or their type depend on keyword arguments (as would be the case in Matlab). Also, attaching _full! and friends to the name made it clear that these are different functions from the ones in LinearAlgebra, and would thus avoid name clashes. There could still be an unexported svd!( ; full=...) which manually dispatches to svd_compact! and svd_full!. However, I am very much open to change my mind based on good arguments.

[lkdvos]

• For the output structs, I am definitely happy with that form.
• Argument order also ok for me, default arguments at the end is definitely a good reason
• For the algorithm types, I definitely agree that kwargs will be necessary to control the implementations. I would prefer to keep the signatures somewhat clean, as with KrylovKit.jl: either you have svd!(input, [output]; kwargs...) or you have the expert mode svd(input, [output], alg/backend). This is mostly since if you have both alg and kwargs, presumably the kind of kwargs depend on the specific algorithm, and that makes it slightly annoying to deal with that from a package development perspective