Implementation of deterministic and stochastic frozen flow integrators

[AdrienAngeAndre]

Hi!
I have been working on stochastic frozen flow integrators, with the goal of creating new efficient intrinsic methods for solving SDEs on manifolds, implementing them, and adding them to the ManifoldDiffEq package in the future. While the analysis is going great, the implementation is more difficult. I am using your packages ManifoldsBase, Manifolds, and ManifoldDiff and the specific homogeneous manifolds of the spheres and the SPD matrices. It is very convenient but I am still struggling with the frozen flow solver you present in ManifoldDiffEq.
Given an orthonormal frame of right-invariant vector fields E_d on my homogeneous manifold M (the standard one that you implemented), the frozen flow integrators rely on the exact flow of y'(t)=\sum_d alpha^d E_d(y(t)), y(0)=p, for given constants alpha^d (see the series of papers by Brynjulf Owren + collaborators of the 90s). You implement it with the Riemannian geodesic exponential and the vector transport operator in

ManifoldDiffEq.jl/src/frozen_solvers.jl

Lines 71 to 86 in 88bc20d

	@doc raw"""
	CG2
	A Crouch-Grossmann algorithm of second order for problems in the
	[`ExplicitManifoldODEProblemType`](@ref) formulation.
	The Butcher tableau is identical to the Euclidean RK2:
	```math
	\begin{array}{c|cc}
	0 & 0 \\
	\frac{1}{2} & \frac{1}{2} & 0 \\
	\hline
	& 0 & 1
	\end{array}
	```
	"""

I struggle to see how these two approaches coincide. Why is the exact solution of the frozen flow given by the Riemannian exponential? I do not think this is true.

[mateuszbaran]

Hi!
Thanks for looking into this. It's been a couple of years since I've implemented it and I remember I've struggled quite a lot to figure out how to translate math to code.

IIRC the trick here was that the vector transport used to define the frozen flow is not parallel transport or anything trying to approximate it but a separate, flat (as in zero Riemann tensor) vector transport. It isn't strictly enforced by the interface and maybe some non-flat vector transports could also result in reasonable integrators -- I haven't investigated it. So, for example, on a Lie group the frame could be a basis of the Lie algebra transported to any point in a left- or right-invariant way (as in flat Cartan-Schouten connections).

Does this make sense to you?

[AdrienAngeAndre]

I see how the use of vector transport and a retraction definitely makes sense from a geometric point of view. I would actually call this a geodesic integrator.
However, I do not think that the order theory carries on. I think the integrators are all of order one as they do not use the exact flow of the frozen problem.
But maybe I am wrong and just do not understand. I would be happy to discuss this in an online meeting with you.

[mateuszbaran]

It's also possible that I've made a mistake somewhere, I would also be happy to have an online meeting to figure this out. You can reach me on Julia's Slack or zulip, or on the numerical differential geometry zulip: https://numdiffgeo.zulipchat.com/ .

So far I've only numerically verified that the frozen solvers in ManifoldDiffEq are indeed of higher order on R^n by comparing results to OrdinaryDiffEq.jl. If you know a good manifold test problem I could implement it here.

[AdrienAngeAndre]

I think that the implementation works on homogeneous manifolds with the (special) orthogonal group as associated Lie group, as the Riemannian and right-invariant one coincide. I think it does not work on the manifold of symmetric positive definite matrices. I also have to try numerical experiments.

[mateuszbaran]

Sure. I've taken another look at those papers and I'm not sure how to implement those integrators properly. Either we choose the frame such that its flow corresponds to a connection we can efficiently compute (via exponential map and vector transport) or we are stuck with numerically solving an ODE. The first class of schemes can already be realized using the current implementation while it's unclear to me how to make something practically useful out of the second class.

[AdrienAngeAndre]

Yes indeed. The whole difficulty lies in choosing the right frame. One would prefer to use high-order integrators relying only on parallel transport and geodesic exponentials, but this is already an open problem on general smooth manifolds. I will implement the new frozen-flow integrators on simple homogeneous manifolds and we can then discuss on how to extend them in general.

[mateuszbaran]

Cool, it would be great to add your contributions to the library 👍

[olivierverdier]

@AdrienAngeAndre funny to bump into you here! If you give more information about what you try to do, I might be able to help. I also use integration on manifolds with ManifoldDiffEq in my projects. Cheers!

[olivierverdier]

@AdrienAngeAndre you can also consider using this Python package for integration on homogeneous spaces. It supports several standard integration methods on manifolds and it's easy to add new ones.

[AdrienAngeAndre]

@olivierverdier Nice to see you here Olivier! Indeed, I am working with Crouch-Grossman/commutator-free/frozen flow methods (I leave RKMK for later). It requires solving exactly frozen ODEs with respect to a frame basis. It is quite easy on matrix Lie groups: the frame basis is Ei(y)=Ai y (matrix multiplication with Ai in the Lie algebra), then the solution of $y'=\sum_i \alpha_i E_i(y)$ is $y(t)=\exp(t\sum_i \alpha_i A_i) y_0$. On matrix homogeneous manifolds, same thing. However, on more general manifolds (such as the symmetric positive definite matrices), the implementation is not straightforward. I thought for some time that I could rephrase it with parallel transport, but it's quite different.