Naive question regarding shortest geodesic

[volkerkarle]

Dear creators of this library,

Before raising this question let me thank you for the amazing work you are doing.

For the shortest geodesic, you write in the documentary: "When there are multiple shortest geodesics, a deterministic choice will be taken." When I now look into the code in exp_log_geo.jl, I find the following,

geodesic(M::AbstractManifold, p, X) = t -> exp(M, p, X, t)
...
shortest_geodesic(M::AbstractManifold, p, q) = geodesic(M, p, log(M, p, q))

i.e. the shortest geodesic is the one with X=log(M,p,q). Hence, the deterministic choice relies on the unique choice of X, given by log, is that correct? Is it actually trivial that X=log(M,p,q) will always lead to the shortest geodesic?

Thanks!

v.

[kellertuer]

Thanks for the nice feedback :)

The story is a little more complicated – the TLDR: The log has the same problem, that it is only locally defined and there are points p,q (on certain manifolds) such that the log has multiple answers (is set valued – mathematically). Then we also require that a deterministic choice is returned for the log as well.

The log might also not be defined (for points p,q), then it might also be that it is not so easy to find the shortest geodesic (though it exists if the manifold is complete).

The example to best illustrate this is the sphere. On the sphere, given two points p,q, the shortest geodesic can even be constructed geometrically – at least if p and q are not opposing points. Take the plane that is span by the origin, p and q. The sphere intersects this plane even in a great circle (the same diameter as the sphere) – and since p and q are not opposite of each other one segment is shorter. This is the shortest geodesic. Now is p and q are opposite of each other – for example north and South Pole – the points N - origin - S lie on a line, so there is more than one plane. So all half-circles are geodesics. The same for the log (given two points, what is the direction (and length) „pointing“ from p to q. This is unique for the first case. For the second (p=N, q=S) any direction is an answer, as long as the length of the vector is $\pi$.

On the sphere, opposing points are the only “problematic” ones. For those we need that a deterministic choice is returned, that is calling log(Sphere(2), Northpole, Southpole) has to always return the same vector (among all of the ones that are valid answers).

[mateuszbaran]

Thank you for your feedback :slightly_smiling_face:

Ronny has a very good answer.

One more point I can add is that shortest_geodesic(M::AbstractManifold, p, q) = geodesic(M, p, log(M, p, q)) is just the default, I think more manifolds should actually implement specialized methods of shortest_geodesic because we can, on some manifolds, calculate it faster than through log. Then we would not require that the exact same geodesic as the one given by log is returned.

For example by default we could cache the solution to the BVP for geodesic in the object returned by shortest_geodesic(M::AbstractManifold, p, q), but the BVP solving is very new so I didn't get around to it yet.

[kellertuer]

Yes, that are two quite important implementation details. Calling first log then exp might be more expensive than a direct computation on a manifold.

The idea here is that (after implementing exp and log) geodesics „just work“ (but might not use the most efficient variant).

[volkerkarle]

So following from what you said, the problem of finding the shortest geodesic is directed towards finding the correct logarithm. I was raising this question because I am working with a complicated manifold defined by the roots of some non-polynomial function and I was wondering what I need to implement to get the shortest geodesic. In the end, I guess there is no free lunch and I cannot circumvent the problem of solving the geodesic equation on that manifold manually (i.e. numerically) and inferring the log and exp map from the results. Thanks again for your comments!

[kellertuer]

Interesting and good luck!!

Maybe a little more precise – there is actually two ways to obtain geodesics: 1) as an ODE with boundary values (boundary value problem) - or more intuitively: p and q are given find the shortest geodesic 2) initial values are given, that is the start point p and the derivative of the curve (a tangent vector) – this would actually be the exponential map as well – and if the tangent vector is long enough (or in other words if you follow the geodesic for long enough) it is a geodesic but not a shortest one anymore. For the sphere again, just imaging running 3/4 around the equator. You are following a geodesic but it is not the shortest one between your start and end point (after 3/4 that is the remaining 1/4 that is shortest ;))

So you do not need the log necessarily but a good solver for the ODE for sure. Actually we have a solver for the second case at least https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/connection.html#Manifolds.solve_exp_ode-Tuple{Any,%20Any,%20Any} if you provide the metric or the connection (locally).

[mateuszbaran]

Actually we have a solver for the second case at least https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/connection.html#Manifolds.solve_exp_ode-Tuple{Any,%20Any,%20Any} if you provide the metric or the connection (locally).

I wouldn't recommend that particular solver, I don't like it. The new solvers (which I like) are:

1. For the boundary value problem: https://juliamanifolds.github.io/Manifolds.jl/latest/features/atlases.html#Manifolds.solve_chart_log_bvp-Tuple{AbstractManifold,%20Any,%20Any,%20AbstractAtlas,%20Any}
2. For the initial value problem: https://juliamanifolds.github.io/Manifolds.jl/latest/features/atlases.html#Manifolds.solve_chart_exp_ode-Tuple{AbstractManifold,%20Any,%20Any,%20AbstractAtlas,%20Any}

The second one also supports switching charts. If you need a shortest geodesic between two points, then the first function does exactly that. One thing you do need is a chart so solve that in and Christoffel symbol of the second kind (which is defined by affine_connection in Manifolds.jl). This script demonstrates it for a torus: https://github.com/JuliaManifolds/Manifolds.jl/blob/master/test/manifolds/torus_example.jl . Ronny is working on turning that example into a nice Pluto notebook which has some additional comments: https://github.com/JuliaManifolds/Manifolds.jl/pull/534/files#diff-0d749ba9a5c20beb599d0eb8a00a179f7c18132e188a295721d2edf08e589356 (working-in-charts.jl).

So this is something you can start with.

[volkerkarle]

Ok, thanks a lot for your comments! That clarifies it all.