An alternative algorithm for the logarithm map of Stiefel manifold.

[dnguyend]

Hello, Per Professor Bergmann's suggestion, I have implemented the first cut of a method to compute the Log function on the Stiefel manifold with the Submersion metric proposed in my paper [1].

I forked a version of Manifolds.jl to a repository and added the implementation. Examples and some test results are in a notebook.

I also implemented the function expm_frechet in Julia, I understand it is not available in the LinearAlgebra package. I follow the open-source code in scipy.linear_alg. At some point, I would like to make it part of LinearAlgebra.

Thank you all for creating a very nice, professional framework to study manifolds. I will be very happy if I could contribute to the effort. I am new to this, I try my best to follow the contributing guidelines but I am sure there will be mistakes - please feel free to send feedback. I am happy to make changes to conform with your high standard.

Du [1] Nguyen, D. Closed-form Geodesics and Optimization for Riemannian Logarithms of Stiefel and Flag Manifolds. J Optim Theory Appl 194, 142–166 (2022). https://doi.org/10.1007/s10957-022-02012-3

[sethaxen]

Thanks! I haven't looked through the algorithm in the paper yet, but I really appreciate you benchmarking the performance vs our implementation of the shooting method. Sounds like we should consider this method as a new default for the logarithm. A few notes:

• Yes, LinearAlgebra doesn't implement the Frechet derivative of exp, and this intentional. Previous discussions (see https://github.com/JuliaLang/julia/issues/5840#issuecomment-1025033449) have concluded that implementations of matrix functions belong in external packages unless that would be type piracy. One result of that discussion is that a standard implementation of the Frechet derivative of exp was added to ChainRules (see https://github.com/JuliaDiff/ChainRules.jl/blob/2cc27e2124fc7746afdd7753fbc9f19c1c24ac13/src/rulesets/LinearAlgebra/matfun.jl#L103-L280). It's executed like this:

  using ChainRules, LinearAlgebra
  using ChainRules.ChainRulesCore
  A, ΔA = randn(100, 100), randn(100, 100)
  Y, ΔY = frule((NoTangent(), ΔA), LinearAlgebra.exp!, A)

  This might be moved to its own package though, since ForwardDiff needs something like this and does not have ChainRules as a dependency (https://github.com/JuliaDiff/ForwardDiff.jl/issues/174). We also don't have ChainRules as a dependency, but I suspect it isn't a heavy one.

• My main concern using this method as a default is adding Optim as a dependency. We should check if that has a measurable increase on load time. If so, it might be better to offer it as an alternative method.
• Do you use a factorization similar to our "Stiefel factorization"? (https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/stiefel.html#Internal-types-and-functions).

[dnguyend]

Thanks, Seth for your thoughtful response.

• No worries about alternative v.s. the main method, I hope people will try it out.
• Yes, I use something similar to Stiefel factorization. Inspired by [2,3] I express the logarithm vector $X = p A + Q R$, where $Q$ forms the complement basis of the span of $[p, q]$, where $p$ is the initial point and $q$ is the terminal point. If the manifold is in $R^{n\times k}$ then usually $Q$ has $k$ columns, however, it could have less. I believe $X = p A + Q R$ is essentially the Stiefel factorization. I determine $Q$ from a SVD in the very first step. In the special case $q = p \exp(A)$ for a $k\times k$ matrix $A$, $Q$ has zero column, and we need to solve for $A$ only. I solve for $(A, R)$ using the optimizer. Stiefel factorization is a very nice trick! In the paper I argue a bit that this is similar to the case of the great circle, the span of $[p, q]$ intersects the Stiefel manifold along a smaller Stiefel manifold, and we only need to work on this small Stiefel manifold, but didn't think about being able to do it in code - so it was nice seeing that implementation! [2] Edelman, A., Arias, T.A., Smith, S.T.: The geometry of algorithms with orthogonality constraints. SIAM J. Matrix Anal. Appl. 20(2), 303–353 (1999). https://doi.org/10.1137/S0895479895290954 [3] Zimmermann, R.: A matrix-algebraic algorithm for the Riemannian logarithm on the Stiefel manifold under the canonical metric. SIAM J. Matrix Anal. Appl. 38(2), 322–342 (2017). https://doi.org/10.1137/16M1074485

[dnguyend]

And thanks for the tip about ChainRules - will look it up!

[kellertuer]

Hi, thanks for the nice comparison. You implementation seems to be quite fast.

My main concern would (as well) be the dependency on Optim, since this is quite a large dependency.

[dnguyend]

Thank you! I didn't see this until now. In the python version, I implemented a customized LBFGS. I can also port it to Julia. It will take another weekend. Du

On Mon, Oct 24, 2022 at 8:38 AM Ronny Bergmann @.***> wrote:

Hi, thanks for the nice comparison. You implementation seems to be quite fast.

My main concern would (as well) be the dependency on Optim, since this is quite a large dependency.

— Reply to this email directly, view it on GitHub https://github.com/JuliaManifolds/Manifolds.jl/issues/546#issuecomment-1288972118, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEOOGFOQBJR5P3H74HQYKZ3WEZ7LLANCNFSM6AAAAAARMPIAWM . You are receiving this because you authored the thread.Message ID: @.***>

-- Thanks, Du

[dnguyend]

I managed to remove the dependency on Optim, using a customized version of LBFGS. The result is saved in my forked version of Manifolds.jl, with an updated workbook containing tests and explanation. I also moved the old workbook and test data out to a different repository to keep Manifolds.jl clean.

The new version is slightly slower and less robust than the version with Optim. This may be because the line search I use out of convenience is not as good as the More'-Thuente search. It converges less well than if we use Optim on the Rosenbrock function. I hope this is good enough for now. The shooting method does perform well for large distances and large manifolds, I think a hybrid version, using the Frechet derivative-based gradient method to bring the cost down to 1e-3, then using the shooting method for the remaining iterations probably will take advantage of both approaches. Potentially we can try this on other manifolds with closed-form geodesics.

Hope you all have a nice weekend!

[kellertuer]

Thanks for the update. I will have to take a look at the new version, but this week I am a little busy (more busy than usual) with teaching.

[dnguyend]

No rush at all!

[kellertuer]

Sorry for the further delay, here a semester is coming to its end and that brings quite a bit of administration with it. I think overall my main problem is, that the algorithm is quite long and complicated, to I fear it will take me a few days to check it, which I find difficult to find time for. Can you maybe open a PR and write tests for it (and maybe an example)?

[dnguyend]

Sure! I totally understand the time constraints, and it is quite a bit of code. The workbook has some tests but I will write tests in the format of the package. Please take your time.

If we include expm_frechet at some point, there are other interesting things we can do - for example we can get Jacobi fields of Stiefel manifolds under these metrics for almost free. With the Riemannian curvature in this article https://www.heldermann.de/JLT/JLT32/JLT322/jlt32027.htm we can have the Jacobi field equation numerically checked, presumably that's another topic.

On Fri, Nov 25, 2022 at 7:02 AM Ronny Bergmann @.***> wrote:

Sorry for the further delay, here a semester is coming to its end and that brings quite a bit of administration with it. I think overall my main problem is, that the algorithm is quite long and complicated, to I fear it will take me a few days to check it, which I find difficult to find time for. Can you maybe open a PR and write tests for it (and maybe an example)?

— Reply to this email directly, view it on GitHub https://github.com/JuliaManifolds/Manifolds.jl/issues/546#issuecomment-1327390709, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEOOGFKRZ7SNXEZRFX4UYN3WKCTD3ANCNFSM6AAAAAARMPIAWM . You are receiving this because you authored the thread.Message ID: @.***>

-- Thanks, Du

[kellertuer]

I think that would be a next topic, yes; also I do not have access to that journal.

[dnguyend]

I got some complains from reviewdog. I reran JuliaFormatter, but can I set up reviewdog locally so I can test before resubmitting the PR ? Appreciate very much if you can give me some pointers! If there is any other thing I should fix - please let me know.

About the paper the final version, shared privately, has much more than the arxiv version https://arxiv.org/abs/2105.01834.