Finding geodesics joining given points

Finding a geodesic joining two given points in a complete path-connected Riemannian manifold requires much more effort than determining a geodesic from initial data. This is because it is much harder to solve boundary value problems than initial value problems. Shooting methods attempt to solve boundary value problems by solving a sequence of initial value problems, and usually need a good initial guess to succeed. The present paper finds a geodesic γ:[0,1]→M\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma :[0,1]\rightarrow M$$\end{document} on the Riemannian manifold M with γ(0) = x0 and γ(1) = x1 by dividing the interval [0,1] into several sub-intervals, preferably just enough to enable a good initial guess for the boundary value problem on each subinterval. Then a geodesic joining consecutive endpoints (local junctions) is found by single shooting. Our algorithm then adjusts the junctions, either (1) by minimizing the total squared norm of the differences between associated geodesic velocities using Riemannian gradient descent, or (2) by solving a nonlinear system of equations using Newton’s method. Our algorithm is compared with the known leapfrog algorithm by numerical experiments on a 2-dimensional ellipsoid Ell(2) and on a left-invariant 3-dimensional special orthogonal group SO(3). We find Newton’s method (2) converges much faster than leapfrog when more junctions are needed, and that a good initial guess can be found for (2) by starting with Riemannian gradient descent method (1).


Introduction
Geodesics on Riemannian manifolds are generalizations of straight lines in Euclidean space, and are fundamental in many areas of mathematics, engineering, and computer science, to name a few. For instance, geodesic regression is used to relate a real-valued independent variable and a manifold-valued dependent data variable [1,2]. Another example is an extension of principal component analysis namely principal geodesic analysis, which is used to study variability of data on a Riemannian manifold [3,4]. A third example is the detection of object boundaries, where geodesic active contours are used to explore relationships between active contours and geodesics via curve evolution theory [5].
Let M be a p-dimensional (1 ⩽ p < ∞) smooth path-connected manifold with a Riemannian metric g. When M is complete with respect to the Riemannian metric, the Hopf-Rinow theorem says that any two points x 0 ,x 1 ∈ M can be joined by a minimal geodesic, namely, a constant-speed curve ∶ [0, 1] → M of minimal length with γ(0) = x 0 and γ(1) = x 1 . For some well-studied examples of Riemannian manifolds, geodesics can be found in closed form, but such cases are rare and in practice numerical methods are usually needed. Even in cases where closed form expressions for geodesics are known, it may be difficult to use these expressions to find geodesics joining two given points. We refer to [31][32][33] for background in differential geometry.
By the fundamental theorem of Riemannian geometry, every Riemannian manifold (M,g) admits a unique Levi-Civita connection ∇, which relates to g by the Koszul formula [6], where X,Y,Z are vector fields on M and [⋅,⋅] is the Lie bracket of vector fields. It turns out that a geodesic joining two given endpoints is precisely a zero-acceleration curve with respect to the Levi-Civita connection ∇, namely the solution of the following 2-point boundary value problem for a system of ODEs For any system of ODEs it is usually much harder to solve a 2-point boundary value problem than an initial value problem. Solutions (often non-unique) to boundary value problems can sometimes be found using shooting methods [7]. In single shooting an unknown initial quantity is guessed, then improved using the terminal error of a solution to an initial value problem. In multiple shooting, multiple initial data are estimated on some range, and then improved. The quality of an initial guess is crucial to the performance of both single and multiple shooting.
In our Riemannian situation a good initial guess can be made when the endpoints x 0 and x 1 are reasonably close. In such cases, the 2-point boundary value problem (2) is solved efficiently using single shooting. In more general nonlocal cases, where x 0 and x 1 are distant, one possible strategy is to divide the interval [0,1] into small subintervals [t i− 1 ,t i ] and update on each subinterval. In [8] the so-called leapfrog 1 3 algorithm 1 is used to find geodesics joining given points (the method has also been used for optimal control problems [9,10]). The basic idea is to treat the local problem as effectively solved, then update the junctions y i := γ(t i ) using a minimal geodesic from y i− 1 to y i+ 1 . Leapfrog always converges to a geodesic (Theorem 5.1 in [8]) and generally performs well in the first few iterations, but convergence can be slow if many iterations are needed. The present paper achieves fast convergence by using single shooting to update junctions in a different way, as follows.
Recall that for the first variation of the energy E ∶= 1 2 ∫ 1 0 g (̇,̇)dt , we have assuming that (for the second equality) each | [t i ,t i+1 ] is a geodesic. So, for piecewise geodesic γ, E vanishes for all variations precisely when all differences ̇(t + i ) −̇(t − i ) are zero for 1 ≤ i ≤ n − 1. This reduces our infinite-dimensional boundary value problem to the following finite dimensional problem.
Setting t i := i/n, we denote by v + i the right-side velocity tangent to M at y i determined by the geodesic i ∶ [0, 1] → M from y i to y i+ 1 . Similarly v − i is the left-side velocity determined by the geodesic from y i− 1 to y i . We may write v + i = log y i y i+1 and v − i = − log y i y i−1 , found by solving (2) with boundary conditions γ(0) = y i ,γ(1) = y i+ 1 and γ(0) = y i− 1 ,γ(1) = y i respectively. The numerical values of these Riemannian logarithms are usually not available from closed form expressions. Then, for the corresponding piecewise geodesic γ we have and E vanishes precisely when (y 1 ,y 2 ,…,y n− 1 ) is a singularity of the vector field F on M × M ×… × M given locally by Note that (y 1 ,y 2 ,…,y n− 1 ) is also the global minimizer of f(y 1 ,y 2 ,…,y n ) for whose minimum value is 0. In the present paper we exploit the fast convergence of Newton's method to find the singularity of the vector field (4). When a good initialization for Newton's method is needed, we use Riemannian gradient descent to minimize the cost function (5). Other variants of gradient descent method, such as Nesterov's accelerated gradient descent [11,12] and accelerated higher-order gradient methods [13] might also be used to minimize (5).
The paper is organized as follows. In Section 2 we review the general framework, of finding a singularity of (4) by Newton's method and minimizing (5) by

General framework for finding geodesics
There may be multiple geodesics joining points x 0 ,x 1 in a general Riemannian manifold M, and the geodesic found by our method may depend on an initial choice for the junctions y 1 ,y 2 ,…,y n− 1 ∈ M. Consecutive junctions should be reasonably close, namely y i ,y i+ 1 should be joined by a unique minimal geodesic with an easy initial guess for single shooting. Because single shooting performs well with a good initial guess, we should not take n to be extremely large: consecutive junctions should be close enough to enable a good initial guess, but no closer.

Riemannian gradient descent for finding geodesics
To use Riemannian gradient descent we need to evaluate the Riemanian gradient of the cost function (5) with respect to all y i 's. Throughout this paper, we denote 2 y i f , R y i f the standard Euclidean gradient and the Riemannian gradient of f with respect to y i respectively, which are related by where ⟨⋅, ⋅⟩ is the Euclidean metric for a coordinate chart containing y i , ⟨⋅, ⋅⟩ y i is the Riemannian metric at y i , and v is an arbitrary tangent vector at y i . As usual when g ∶ M → N is a smooth map, we denote its derivative by dg ∶ TM → TN , where TM and TN are the tangent bundles. The derivative at x ∈ M is denoted by Given junctions y 1 ,y 2 ,…,y n− 1 we set y 0 := x 0 , y n := x 1 , y := (y 1 ,y 2 ,⋯ , where the Riemannian adjoint of the logarithm map derivative is defined by (8) follows from the fact that In the next section we develop strategies to calculate the derivatives and adjoints of derivatives of logarithm maps in the setting of general Riemannian manifolds, codimension-one embedded submanifolds and Lie groups. The Riemannian gradient method for finding a geodesic joining given points x 0 ,x 1 is summarized as Algorithm 1.
In general, the step size η is chosen to make the algorithm converge at a suitable rate. There is usually some trial and error. Many iterations are needed when η is too small, but when η is large there can be problems with overshooting. For a general Riemannian manifold where x 0 ,x 1 are in the same coordinate chart, a reasonable choice for the initial curve ̃ would be a line segment in chart coordinates. Alternatively, when M is embedded in Euclidean space, ̃ could be the orthogonal projection of the Euclidean line segment from x 0 to x 1 .
To implement the shooting method we use either MATLAB's function bvp4c/ bvp5c or Mathematica's function NDSolve with method shooting in Lines 2 and 10 of Algorithm 1.
In Line 7, y i is updated via the exponential map exp as follows where different forms of geodesics equations are needed to solve when it comes to Riemannian manifolds with local coordinates available, or of codimension one in Euclidean space, or left-invariant Lie groups.

Newton's method for finding geodesics
Newton's method is a powerful tool for finding the zeros of nonlinear functions in Euclidean spaces, which has a wide range of applications in mathematics and other subjects. This has motivated studies of the generalization of Newton's method from the Euclidean setting to Riemannian manifolds [23][24][25][26]. For a smooth vector field defined on a Riemannian manifold, Fernandes et al. [19] showed that in a suitable neighborhood of the singularity of the vector field, the sequence generated by Newton's method converges quadratically when the covariant derivative of the vector field is invertible in a convex neighborhood containing the intial guess. However, Newton's method may diverge if the initial guess is not sufficiently close to the solution. Many strategies have been proposed to overcome this drawback, such as incorporating the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, the Levenberg-Marquardt algorithm, the trust region algorithm [27,28], by introducing a merit function, which is the squared norm of the vector field, to replace a divergent iteration generated by Newton's method or make up a failed iteration in Newton's method [20,21]. For simplicity, throughout this paper, we consider the basic Newton's method and leave other variants and improvements for future study. In our experience, the Riemannian gradient descent of Section 2.1 gives a suitable initial guess for Newton's method.
We calculate the derivative ∂F(y) of the vector field F, in terms of derivatives of logarithm maps in Section 2.1 as follows, ∂F(y) = where we have arranged F as a column vector. Then the covariant derivative ∇F(y) is related to ∂F(y) by (6).

3
Newton's method for finding the singularity of the vector field F is to update y i by where exp is the exponential map on manifold M and w = (w 1 ,w 2 ,⋯ ,w n− 1 ) is the solution of Then Newton's method for finding a geodesic joining given points x 0 and x 1 by looking for a singularity of (4) is summarized as Algorithm 2.

The derivative of the logarithm map
In Algorithms 1 and 2 a key task is to find the derivatives of v + i and v − i , namely derivatives of Riemannian logarithms. We write Log(x, y) ∶= log x y . Let v := Log(x,y), then Differentiating y = Exp(x,v) with respect to y along the direction w yields then (11) Differentiating y = Exp(x,v) with respect to x along the direction w gives then So derivatives of the logarithm map are found from derivatives of the exponential, using (13) and (14). When x = y, In practice, when x is extremely close to y (though not exactly the same), using (13) and (14) may cause some numerical issues, and we use (15) instead.
To calculate derivatives of the exponential map, we first consider a variation of the geodesic γ(t) = Exp(x,v) given by γ 1 (t,s) := Exp(Exp(x,su 1 ),tv(s)), where u 1 ∈ T x M defines a variation of the initial point along the geodesic Exp(x,su 1 ). We extend v ∈ T x M as a vector field v(s) along Exp(x,su 1 ) via parallel translation. We also consider a variation of γ(t) given by γ 2 (t,s) := Exp(x,tv + su 2 ), where u 2 ∈ T v (T x M) = T x M. These variations determine Jacobi fields, which are solutions J of the secondorder linear dynamical system [2,6,22] where ∇ is the Levi-Civita connection for the Riemannian metric g on M, and R is the Riemannian curvature tensor. Then where the field J 1 is found by numerically solving (16) with the initial conditions J 1 (0) = u 1 , ̇J 1 (0) = 0 . We also have where J 2 is found by numerically solving (16) with the initial conditions J 2 (0) = 0, ̇J 2 (0) = u 2 .
Next we consider in more detail how to solve the Jacobi equation (16) where M is a general Riemannian manifold, and in the important special case when M is a leftinvariant Lie group (a specialized technique can then be applied).

M is a Riemannian manifold with local coordinates available
Let M be a p-dimensional (1 ⩽ p < ∞) Riemannian manifold with a Riemannian metric g, then g induces a Levi-Civita connection ∇ by the Koszul formula (1). Given a local coordinate system {x i } on an open set U ⊂ M, { i ∶= x i } is a basis for vector fields defined on U. Suppose we have two vector fields X = X i ∂ i and Y = Y j ∂ j , then the expression for the connection ∇ in local coordinates {x i } on U is given by The Christoffel symbols are defined as the functions Γ k ij ∶ U → ℝ given by where 1 ⩽ i, j, k ⩽ p . By the fundamental theorem of Riemannian geometry, Christoffel symbols have the following form with respect to the metric g ij = g(∂ i ,∂ j ) [6], where the matrix g kl is the inverse of the symmetric p × p matrix g kl . We define the symmetric bilinear Christoffel form Γ Example 1 Given a,b,c > 0, we define Ell(2) ⊂ ℝ 3 to be the 2-dimensional ellipsoid given by with Riemannian metric defined as the restriction of the standard Euclidean inner product. This manifold admits the following parameterization where ∈ [0,2π), φ ∈ [0,π] act as global coordinates of Ell (2). We can take advantage of the coordinate system { ,φ} to choose the initial junctions y i in Algorithms 1 and 2.
With respect to the restriction of the Euclidean metric, we have If a = b, then by (21), non-vanishing Christoffel symbols are given by
In summary, Algorithm 3 finds the geodesic joining two consecutive points y i ,y i+ 1 and the logarithm map v + i . Algorithm 4 finds the derivatives of the logarithm maps , where e α is a zero column except 1 in α-th position. All solutions for the type of initial conditions J(0) = e α , ̇J = 0 is arranged as a matrix J 1 (1), and those for the other type is arranged as a matrix J 2 (1). Then, the derivative (d y i Log) y i ,y i+1 (w) = −J 2 (1) −1 J 1 (1)w and (d y i+1 Log) y i ,y i+1 (w) = J 2 (1) −1 w.
The relationship between Algorithm 1 and Algorithm 4 is that Algorithm 4 is a significant sub-algorithm of Algorithm 1 that calculate the derivatives d y i v + i , d y i+1 v + i , which will be used in calculating the Riemannian gradient of the cost function (see (8)).
In order to represent the Jacobi equation (16) in terms of the function h or the unit normal ν, we consider a variation s ∶ (− , ) is the variational vector field along γ(t). Therefore, we have

Proposition 2 The Jacobi equation (16) is equivalent to
As with the geodesic equation, the initial conditions (0),̇(0) ∈ T (0) M force the solution of (27) to stay on the tangent space of M, although ∶ [0, 1] → p+1 is viewed as a curve in Euclidean space.
In summary, Algorithm 6 3 finds the geodesic joining two consecutive junctions y i ,y i+ 1 and the logarithm map v + i . Algorithm 7 finds the derivatives of the logarithm maps d y i v + i , d y i+1 v + i . Algorithm 8 finds the exponential map on M.

M is a left-invariant Lie group
Given an inner product on the Lie algebra of a Lie group G, a left-invariant Riemannian metric on G is defined by left-translation: where (dL g ) h ∶ T h G → T gh G is the derivative of L g at h ∈ G, w 1 ,w 2 ∈ T g G. To shorten notation we may drop subscripts in the metric 〈⋅,⋅〉 if there is no ambiguity. The Lie group G equipped with this Riemannian metric is said to be left-invariant.
To simplify the second-order system (16), we use the technique of left Lie reduction [15]. Given a C ∞ curve x ∶ [0, 1] → G , let F be any C ∞ vector field along x. Then the left Lie reduction F of F is defined as where t ∈ [0,1]. We refer to [16-18, 29, 30] and references for applications of left Lie reduction in reducing high order nonlinear dynamical systems on Lie groups and Riemannian homogeneous spaces.
Let {E i } be an orthonormal basis of at the identity e, then {Ẽ i ∶= (dL g ) e (E i )} forms an orthonormal basis at g ∈ G. By the Koszul formula (1), the constant ⟨Ẽ i ,Ẽ j ⟩ and the left invariance of the metric and vectors, we have (see page 118 in [14]) from which where we have used the Jacobi identity in the last equality.
Lemma 3 Let x ∶ [0, 1] → G be a C ∞ curve, V the left Lie reduction of the velocity vector field ̇x of x, W 0 any C ∞ vector field along x and W its left Lie reduction. Then ) e E i . By (31)  (34)

Numerical experiments
Next we illustrate our algorithms by numerical experiments on the ellipsoid Ell (2) and on a left-invariant special orthogonal group SO (3). All experiments were executed in MATLAB 2020b on a computer with Intel Core i7-8700K CPU, 32 GB RAM and Windows 10 Enterprise. MATLAB codes are available on github https:// github. com/ beili ren/ Endpo int-geode sic-probl ems.

The ellipsoid
Let Ell(2) be the 2-dimensional ellipsoid defined in Example 1. Define the diagonal matrix A ∶= diag{ 1 a 2 , , whose tangent space is given by T x Ell(2) = {v ∈ ℝ 3 |v T Ax = 0} . The normal space   (27), the Jacobi equation is where x is a solution of (37). For two points x,y on Ell (2), let TProj(y − x) be the projection of the vector y − x to the tangent space T x Ell (2). Then, TProj(y − x) = y − x − r 0 Ax for some r 0 ∈ ℝ .
Solving the geodesic (37) with linear initial guess by shooting is unsuccessful (as warned by MATLAB that the maximum residual at endpoint is very large with affordable mesh points), which indicates the necessity of dividing the interval [0,1] into several sub-intervals. In this experiment, we find when the number of sub-intervals n ≥ 4, shooting with the initial guess (39) is able to find the geodesic on each sub-interval. Figure 1 shows the geodesic joining x 0 and x 1 generated by Newton's method (Algorithm 2) with n = 5 and = 10 − 4 , where magenta and blue points represent initial and final junctions respectively, green and red curves represent initial and final piecewise geodesics respectively.
To compare Riemannian gradient descent (Algorithm 1) and Newton's method with the leapfrog algorithm (Algorithm 12), we vary the number of sub-intervals n from 4 to 9, set the step size η = 0.001 in Algorithm 1 and the time limitation for all Ax. algorithms is 100 seconds. The runtime for different methods is reported in Table 1, from which we find Riemannian gradient descent is the slowest, leapfrog method is getting slower and slower as n increases, and Newton's method is the most efficient one in general. Note that for a function f defined on Ell (2), its Riemannian gradient R x f is related to its Euclidean gradient ∂ x f by R x f = (I 3 − xx T A) x f . If we denote the operator (d x Log) x,y and (d y Log) x,y on Ell(2) by matrices, then their adjoints are given by matrix transposes.

Fig. 2
The red and black frame represent rotation matrix X 0 and X 1 respectively. Geodesic or junctions on the manifold SO (3) are represented by endpoints of frames rotated by elements in SO (3). Magenta and blue points denote initial and final junctions respectively and green curve denotes the final piecewise geodesic it is necessary to solve the Jacobi equation. On a general Riemannian manifold: this equation can be rewritten in terms of Christoffel symbols if local coordinates are available or a embedding function if the codimension of the manifold is one. On a Lie group, the Jacobi equation can be reduced to a simpler differential equation in terms of a bilinear transformation of the Lie algebra. Finally, we test the effectiveness of the proposed methods by finding geodesics joining two given points in the 2-dimensional ellipsoid Ell (2) with Euclidean metric and in the special orthogonal group SO(3) with a left-invariant metric.

Appendix: Leapfrog algorithm for finding geodesics
The following algorithm is adapted from [8].
material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.