Expose algorithm kwargs for Grassmann methods

[Gertian]

fixed tiny bug

[lkdvos]

These test failures seem unrelated, we should probably check if the master even runs...

As a side-comment, this is not really a bug. The SDD algorithm should be a lot faster for large matrices, and this change really does have some performance implications. Nevertheless, I agree that in our case, as the stability of SDD turns out to be insufficient for the cases at hand, we should probably default to the (slower) more stable SVD algorithm

[Jutho]

We should probably adopt a SafeSVD version that does

try
   # use SDD
catch
   # fall back to SVD
end

A singular value decomposition is probably sufficiently costly in most cases that the try ... catch overhead is negligible.

[Gertian]

@lkdvos I'm quite certain that the failing tests are indeed unrelated since I've ran many hours of Gradient Descent with this change without any issues.

@Jutho This was my original patch, the problem is that SDD already changes the input tensor so we should really do :
Make hardcopy(tensor) try catch
but then there is a substantial overhead which is why I finally just accepted to SDD.

[Jutho]

In this case without the exclamation mark, tsvd should not change the input parameter. I agree that it is more tricky with the mutating method. So maybe SafeSVD should only be valid for tsvd but not for tsvd!.

[Gertian]

Hi all,

I updated the PR.

From the TensorKitManifolds perspective the only changes compared to master are now that all methods (I think) now take the alg keyword which can be SVD or SDD.

I simultaneously made a PR in MPSKit which ads an extra argument svd_alg to the GradientGrassman struct. This can be set to SVD or SDD and is passed to the underlying TensorKitManifolds algorithms.

[lkdvos]

I added some changes to make the tests run again:

• getproperty is treated specially in the Julia compiler, in the sense that it has a default aggressive constprop. This ensures that it is type stable with the different symbols. Adding a keyword argument breaks this, so I had to revert that change. In principle, we could consider manually adding a Base.@constprop :aggressive, but somehow I don't think getproperty with a keyword is very Julian.
• I refactored the default algorithm out, and changed this from SDD to SVD. I think the stability issues are severe enough to warrant this.

[Gertian]

This is probably a ridiculous edge case but I'm wondering what might happen if a user sets the TensorKitManifolds default svd algorithm to SDD and then uses the dispatching in MPSKit to use SVD everywhere.

In this case the SDD call inside getproperty might still lead to a crash ?

Or do you think users should simply not use SDD anymore ?

[lkdvos]

The current approach does not actually have an option to set the default algorithm, it is a compile time constant.

Note that your approach also did not really expose the alg keyword in the getproperty function, this is the function that gets called when you access fields of a tangent: y.U is completely equivalent to getproperty(y, :U). (This is also why this function is treated differently).

Thus, the current implementation has the stable one as default, and if you would want to change that, it would require quite a bit of rewrite. I am more or less okay with anything, so I would wait for @Jutho 's opinion

[Gertian]

Ok, that makes sense thanks !

[Jutho]

My apologies for not getting to this sooner, which will probably result in double work. I have the following requests / suggestions:

• I would define a default_svd_alg(A::TensorMap) = SVD() in the main module file ("src/TensorKitManifolds.jl")
• I would change the default value of the alg keyword of projectisometric! (also in the TensorKitManifolds.jl file) to defualt_svd_alg(A). Also use this value in any other tsvd and tsvd! call that you can find, I think they are only in "src/grassman.jl".
• I would then remove the alg value (reset it to nothing) in Grassmann.retract and friends (invretract, relativegauge). This alg argument is actually aimed at different retraction schemes (which exist for example for the Stiefel case).

This makes all svd calls in TensorKitManifolds.jl completely hard coded, but you can still change it by redefining TensorKitManifolds.default_svd_alg(::TensorMap) = SDD().

Ultimately, we should have algorithms structures at the level of retract and transport and friends, where then the algorithm to be used for the svd would be a field in that structure. But that should be part of a larger overhaul, possibly where we consider to buy in to the Manifolds.jl ecosystem.

[Gertian]

Thanks @Jutho for the suggestions and @lkdvos for implementing them faster then I could even read them 😓

Does this mean that

using TensorKitManifolds
TensorKitManifolds.default_svd_alg(::TensorMap) = SDD()


using MPSKit

rest of my code

will result in using the stable version of svd inside all MPSKit algorithms ?

[lkdvos]

No, SDD() is the unstable/fast algorithm, so by default the stable version will already be selected, you won't have to do anything.

[Gertian]

Oh yes, indeed. Thanks for pointing that out !

[Jutho]

Great. I will merge when lights turn green. There seems one svd! call hidden in "src/auxiliary.jl" in _stiefellog, but that is unused I assume so let's keep that for another time.

Also, we should probably switch to the LinearAlgebra name of these algorithms in TensorKit.jl at some point. default_svd_alg is also defined in LinearAlgebra, but it might actually be useful to have independent default_svd_alg routines in the different packages for more granular control, until we have a batter strategy implemented.

[lkdvos]

I saw the talk about ScopedValues on JuliaCon last week, it seems like these kinds of defaults and configurations are precisely what is being targeted with that, and from what I understand, this would solve our issues. Just something to keep in the back of our heads when we revisit this.