Implement second order approximation to geodesic for the Connection Manifold

[Nimrais]

Currently, the package implements the exponential map through an ODE solver (ODEExponentialRetraction), which uses OrdinaryDiffEq under the hood. This implementation is quite computationally intensive. It would be beneficial to have a more computationally efficient approximation.
One popular approach, especially in information geometry, is to use a second-order approximation. For reference, see equation 16 in section 5.3 of this paper: https://proceedings.mlr.press/v119/lin20d/lin20d.pdf. This approximation is relatively easy to implement and doesn't require any heavy dependencies.

it seems not to much of work and it will come just from a careful implementation of smt like this

function exp_secondorder(
   ...
    Γ,
    p0,
    v0
)   
    Δ = similar(p0)  # Preallocate Δ with same type/size as p0
    Manifolds.@einsum Δ[k] = -0.5 * Γ[k,i,j] * v0[i] * v0[j]
    return p0 + v0 + Δ
end

I want to propose implementing this approximation for the Manifolds.jl package if you find it a reasonable proposal. I can file a PR myself.

[mateuszbaran]

Good idea. In the context of this paper I think it would be appropriate to make a new <:AbstractRetractionMethod struct for this method. I'd suggest using the affine_connection! function instead of getting the coefficients of the Christoffel symbol, it can often be more efficient.

[kellertuer]

Hi,
this sounds like a very good idea! We would have to come up with a good name SecondOrderRetraction for example already has another meaning as well (see Absil, Mahony, Sepulchre p. 107).

One caveat here is, that the current formula you propose does not necessarily work for all manifolds,
namely in the paper above you need the Fisher-Rao metric (which we do not have as a type in Manifolds.jl),
and work in what they call natural parametrisation with respect to that.

[Nimrais]

One caveat here is, that the current formula you propose does not necessarily work for all manifolds,
namely in the paper above you need the Fisher-Rao metric.

Do you mean the formula above that I wrote in the issue (it seems a generic rule that they just computed for Fisher-Rao geometry), or are you talking about the equation (16) that indeed involves the Fisher information?

[kellertuer]

Ah, their Equation (16) seems to be more general, sure. But I am not sure what they would require especially for v0 (in Manifolds.jl usually denoted X)? I would think it at least requires an AAbstractBasis of T_pM and you would need som coordinates of X in which you would get the christoffel symbols and then how would you add that to p (p0)?

So from the theory, you might be right, this is maybe doable relatively generic, in practice, I am not fully convinced yet.

[Nimrais]

ahh, sure how would you add that to the p. Exponential family manifolds always have the ambient R^n around them, so we can always interpret it as just Euclidian +, so here you are right; the formula with the + works only for manifolds with some ambient space around them where this + indeed is defined.

[kellertuer]

And for a lot of manifolds this should work, but my most prominent example are either
the fixed rank matrices – represented by their three SVD factors and not their matrix in the embedding – or Grassmann, where we represent the equivalence class by one Stiefel element and + would be defined but not so clear what it is.

So in theory, this is super straight forward, if you have an embedding, in practice – we need a basis for the tangent space and might have to be a bit careful.

[mateuszbaran]

But we don't need the ambient space here, p0 + v0 + Δ is in a parametrization.

So in theory, this is super straight forward, if you have an embedding, in practice – we need a basis for the tangent space and might have to be a bit careful.

We specifically have InducedBasis for that.

[kellertuer]

I do not see how an induced Basis allows us to add the coordinates with respect to that basis (is that what you mean?) to a point.
For example on the Sphere, the coordinates in the induced basis are in d-dimensions, while p0 is in d+1 dimensions. The point p0 would never be in a parametrisation?

[mateuszbaran]

I mean you add parameters of p in the selected parametrization and coordinates of a tangent vector in the induced basis.

[kellertuer]

Ah, so we need at least a chart, if not an atlas and then additionally a basis (induced by the chart) of the tangent space? Yes, that would work, though I feel our support for (concrete) atlases is also still a bit limited.

edit: but right, I was not careful enough in thinking in parametrizations. Maybe the same way we use c for coordinates (I try to unify that every now and then) we should use some letter for parameters?

[mateuszbaran]

Yes, ideally an atlas would be needed.

I think our support for that is good enough, see for example EmbeddedTorus. We have an example of working with InducedBasis and in-chart computations in a tutorial. Let me know if there is anything specific you think is missing.

[kellertuer]

I mostly ment actually available concrete charts / atlases for manifolds.
The interface is fine, as far as I see, but we only have very few manifolds with an actual atlas available.

[mateuszbaran]

Sure, we need more atlases. But IMO this is the right way of implementing this feature.

[kellertuer]

Ah yes, sorry if I was unclear on that.
In parameters x for p, and in coordinates c of the induced basis for the tangent vector X, Δ being similar to c and the “Christoffelistic” being with respect to that basis – this is the completely generic right approach for this.
That is all the “be careful with the details” was referring to before and this should be careful enough. 👍

[mateuszbaran]

@Nimrais Do you need some help implementing the approach I've sketched? I can implement an example for sphere with StereographicAtlas so you could just focus on extending support to the manifolds and atlases you need.