A discrete framework to find the optimal matching between manifold-valued curves

TL;DR

This paper introduces a discrete Riemannian framework for optimal matching of manifold-valued curves, providing an efficient algorithm that converges to the continuous model and is applicable to various geometries like hyperbolic, Euclidean, and spherical spaces.

Contribution

It develops a novel discrete Riemannian structure for shape analysis of manifold-valued curves, with proven convergence to the smooth case and practical algorithms for optimal matching.

Findings

The discrete framework accurately approximates the continuous model.

The algorithm effectively matches curves in hyperbolic, Euclidean, and spherical spaces.

Comparison with dynamic programming shows competitive performance.

Abstract

The aim of this paper is to find an optimal matching between manifold-valued curves, and thereby adequately compare their shapes, seen as equivalent classes with respect to the action of reparameterization. Using a canonical decomposition of a path in a principal bundle, we introduce a simple algorithm that finds an optimal matching between two curves by computing the geodesic of the infinite-dimensional manifold of curves that is at all time horizontal to the fibers of the shape bundle. We focus on the elastic metric studied in the so-called square root velocity framework. The quotient structure of the shape bundle is examined, and in particular horizontality with respect to the fibers. These results are more generally given for any elastic metric. We then introduce a comprehensive discrete framework which correctly approximates the smooth setting when the base manifold has constant sectional curvature. It is itself a Riemannian structure on the product manifold of "discrete curves" given by a finite number of points, and we show its convergence to the continuous model as the size of the discretization goes to infinity. Illustrations of optimal matching between discrete curves are given in the hyperbolic plane, the plane and the sphere, for synthetic and real data, and comparison with dynamic programming is established.

Figures

Tables

Table 1. Table 1: Length of the geodesics shown in the first three lines of Figure 8 , displayed in the same order.

4.400	4.349	10.785	10.478
4.300	3.000	9.680	9.305
4.301	3.000	9.691	9.307

Equations

$$\mathcal{M}=\{c\in C^{\infty }([0,1],M)|c^{\prime }(t)\ne 0,\forall t\in [0,1]\},$$

$$T_c\mathcal{M}=\{w\in C^{\infty }([0,1],TM)|w(t)\in T_{c(t)}M,\forall t\in [0,1]\}.$$

$$c\mapsto \left(x:=c(0),q:=\frac{c^{\prime }}{\sqrt{|c^{\prime }|}}\right).$$

$$G_c(w,w)=|w(0){|}^2+{\int }_0^1|{\nabla }_s{q}^w(0,t){|}^2\mathrm{d}t.$$

$${\nabla }_s{q}^w(0,t)=|c^{\prime }{|}^{-1}({{\nabla }_{t}w}^N+\frac{1}{2}{{\nabla }_{t}w}^T).$$

$$G_c$$

$$G_c^{a,b}(w,w)=|w(0){|}^2+\int a^2\left(|{\nabla }_{\ell }{w}^N{|}^2+b^2|{\nabla }_{\ell }{w}^T{|}^2\right)\mathrm{d}\ell .$$

$$G_{c\circ \phi }(w\circ \phi ,z\circ \phi )=G_c(w,z),\forall w,z\in T_c\mathcal{M}.$$

$$d_G(c_0,c_1)=\inf \left\{L(c):c(0)=c_0,c(1)=c_1\right\},$$

$$s\mapsto \left(x(s):=c(s,0),q(s,\cdot ):=\frac{c_t(s,\cdot )}{|c_t(s,\cdot ){|}^{1/2}}\right),$$

$$L(c)={\int }_0^1\sqrt{|x_s(s){|}^2+{\int }_0^1|{\nabla }_{s}q(s,t){|}^2\mathrm{d}t}\mathrm{d}s.$$

$$E(c)={\int }_0^1\left(|x^{\prime }(s){|}^2+{\int }_0^1|{\nabla }_{s}q(s,t){|}^2\mathrm{d}t\right)\mathrm{d}s.$$

$${\nabla }_s{x}_s(s)+r(s,0)=$$

$${\nabla }_s^{2}q(s,t)+|q(s,t)|\left(r(s,t)+r(s,t{)}^T\right)=$$

$$r(s,t)={\int }_t^1\mathcal{R}(q,{\nabla }_{s}q)c_s(s,\tau {)}^{\tau ,t}\mathrm{d}\tau ,t\in [0,1].$$

$$J(s)={\frac{\partial }{\partial a}|}_{a=0}c(a,s).$$

$$T\mathcal{M}=\text{Ver}\oplus \text{Hor}$$

$${\text{Ver}}_c$$

$$\mathcal{C}$$

$${\text{Hor}}_c={\left({\text{Ver}}_c\right)}^{{\perp }_G}=\{h\in T_c{\mathcal{M}}_f:G_c(h,mv)=0\forall m\in \mathcal{C}\}.$$

$$G_c(w^{hor},z^{hor})={\bar{G}}_{\pi (c)}\left(T_c\pi (w),T_c\pi (z)\right),$$

$$\bar{d}\left(\overline{c_0},\overline{c_1}\right)=\inf \left\{d\left(c_0,c_1\circ \phi \right)|\phi \in {\text{Diff}}^{+}([0,1])\right\}.$$

$$\bar{d}\left(\overline{c_0\circ \phi },\overline{c_1\circ \psi }\right)=\bar{d}\left(\overline{c_0},\overline{c_1}\right),\forall \phi ,\psi \in {\text{Diff}}^{+}([0,1]).$$

$$\left((a/b{)}^2-1\right)⟨{\nabla }_{t}h,{\nabla }_{t}v⟩-⟨{\nabla }_t^{2}h,v⟩$$

$$+|c^{\prime }{|}^{-1}⟨{\nabla }_t{c}^{\prime },v⟩⟨{\nabla }_{t}h,v⟩=0.$$

$${\text{Hor}}_c=\{h\in T_c\mathcal{M}:\forall t\in [0,1],3⟨{\nabla }_{t}h,{\nabla }_{t}v⟩$$

$$-⟨{\nabla }_t^{2}h,v⟩+|c^{\prime }{|}^{-1}⟨{\nabla }_t{c}^{\prime },v⟩⟨{\nabla }_{t}h,v⟩=0\}.$$

$$m^{\prime \prime }-⟨{\nabla }_t{c}^{\prime }/|c^{\prime }|,v⟩m^{\prime }-4|{\nabla }_{t}v{|}^{2}m$$

$$=⟨{\nabla }_t^{2}w,v⟩-3⟨{\nabla }_{t}w,{\nabla }_{t}v⟩-⟨{\nabla }_t{c}^{\prime }/|c^{\prime }|,v⟩⟨{\nabla }_{t}w,v⟩.$$

$$G_c^{a,b}(h,mv)=⟨h(0),m(0)v(0)⟩$$

$$+{\int }_0^1(a^2⟨{\nabla }_t{h}^N,{\nabla }_t(mv{)}^N⟩+b^2⟨{\nabla }_t{h}^T,{\nabla }_t(mv{)}^T⟩)\frac{\mathrm{d}t}{|c^{\prime }|}$$

$$={\int }_0^1(a^{2}m⟨{\nabla }_{t}h,{\nabla }_{t}v⟩+b^2{m}^{\prime }⟨{\nabla }_{t}h,v⟩)\frac{\mathrm{d}t}{|c^{\prime }|}$$

$$\displaystyle=\int_{0}^{1}\!a^{2}m\langle\nabla_{t}h,\nabla_{t}v\rangle\frac{\mathrm{d}t}{|c^{\prime}|}-\!\int_{0}^{1}\!b^{2}m\frac{d}{dt}\big{(}\langle\nabla_{t}h,v\rangle|c^{\prime}|^{-1}\big{)}\mathrm{d}t$$

$$\displaystyle=\int_{0}^{1}m/|c^{\prime}|\Big{(}(a^{2}-b^{2})\langle\nabla_{t}h,\nabla_{t}v\rangle-b^{2}\langle\nabla_{t}^{2}h,v\rangle$$

$$\displaystyle\hskip 105.00015pt+b^{2}\langle\nabla_{t}c^{\prime},v\rangle\langle\nabla_{t}h,v\rangle|c^{\prime}|^{-1}\Big{)}\mathrm{d}t,$$

Full text

1 Introduction

1.1 Context

The study of curves and their shapes is a research area with numerous and varied applications, which is why it has known a great deal of activity over the past years. These curves can be closed or open, and take their values in a Euclidean space or more generally in a Riemannian manifold. To name a few examples, closed plane curves are central in shape analysis of objects younes12 ; the study of trajectories on the Earth requires to deal with open curves on the sphere zhang16 ; and in signal processing, locally stationary Gaussian processes can be represented by open curves in the hyperbolic plane, seen as the statistical manifold of Gaussian densities moi16 , moi17 . Here we are concerned with the study of open curves in a manifold $M$ of constant sectional curvature.

There are naturally many ways to go about comparing curves in a manifold. One way is to see the space of manifold-valued curves as an (infinite-dimensional) manifold $\mathcal{M}$ itself, and equip it with a Riemannian metric $G$. Then the geodesic between two curves in $\mathcal{M}$ describes how one optimally deforms into the other, while its length gives a measure of dissimilarity : the geodesic distance. The advantage of this strategy is that it provides all the convenient tools of the Riemannian framework. An interesting property for the metric $G$, from the point of view of the applications, is invariance with respect to re-parameterization: for closed curves, this amounts to considering only the shape of an object; for an open curve representing the evolution in time of a given process, this allows us to analyze it regardless of speed or pace. A popular strategy is to consider the quotient space $\mathcal{S}$ of curves modulo reparameterization, where two curves are considered identical if they pass through the same points of $M$ but at different speeds, or equivalently when one can be obtained by reparameterizing the other. This quotient space is often called the shape space. If the Riemannian metric $G$ defines the same scalar product at all points of $\mathcal{M}$ which project on the same "shape", then $G$ induces a Riemannian structure on the quotient space. Computing geodesics between two shapes for that metric can be done in two steps : (1) establishing an "optimal matching" through the choice of two parameterizations, one for each shape, that will put in correspondence the points on the first shape with points on the second shape, as in Figure 1, and (2) computing the geodesic between the two parameterized curves obtained. Both steps depend on the choice of the metric.

1.2 Previous work

Since the simplest metric one can think of, the $L^{2}$-metric (slightly modified to stay constant along the fibers), induces a vanishing distance on the quotient space michor05 , different classes of metrics have been studied to perform shape analysis. The large class of Sobolev metrics involves higher order derivatives to overcome the vanishing problem of the $L^{2}$-metric michor07 . A first-order Sobolev metric for plane curves was introduced in younes98 and used in younes08 , which can be mapped to the $L^{2}$-metric by a change of coordinates, namely by considering the complex square root of the speed of the curve. This metric was modified by the authors of mio07 to define the family of elastic metrics $G^{a,b}$, parameterized by two constants $a$ and $b$ which control the degree of bending and stretching of the curve. In srivastava11 , the authors show that for a certain choice of parameters, the elastic metric can again be mapped to the $L^{2}$-metric using the so-called square root velocity (SRV) coordinates, where a curve is represented by its speed renormalized by the square root of its norm. The SRV framework was generalized in bauer12 for any elastic metric with weights $a$ and $b$ satisfying a certain relation. A quotient structure for the metric used in srivastava11 is carefully developed in lahiri15 , where the authors prove that if at least one of two curves is piecewise-linear, then there exists a minimizing geodesic between the two, and give a precise algorithm to solve the matching problem. In bruveris15 , it is proven that in the same framework, there always exists an optimal reparameterization realizing the minimal distance between two $C^{1}$ plane curves. Another approach is proposed in tumpach16 , where the authors restrict to arc-length parameterized curves and characterize the horizontal space of the quotient structure for these curves in the elastic framework.

Concerning manifold-valued curves, the geodesic equations for Sobolev metrics in the space of curves and in the shape space were given in bauer11 in terms of the gradient of the metric with respect to itself. A generalization of the SRV framework to manifold-valued curves was introduced in zhang15 and used in zhang16 , while another one was proposed in moi17 . Extension to curves in a Lie group or a homogeneous space can also be found in celledoni16 , celledoni17 , su17 . Both metrics in zhang15 and moi17 coincide with the metric of srivastava11 in the flat case, however they define different Riemannian structures when the base space has curvature. The difference lies in the way computations are made – in zhang15 and zhang16 they are moved to the tangent spaces at the origins of the curves, resulting in simpler computations, whereas in moi17 they are done directly in the base manifold $M$, transporting data pointwise across $M$ from one curve to the other, thus making the comparison more sensitive to the local "geography" of $M$. When comparing the geodesic distances, this difference is measured by a curvature term (see moi17 , Prop. 3). In addition, unlike the metric of zhang15 , the one in moi17 can be written as an elastic metric (with $a=2b=1$). The work of zhang15 is applied to curves in the space of symmetric positive definite matrices, while the case of spherical trajectories is investigated in zhang16 , where the authors exhibit simplifications. In both cases, optimal matching between curves is achieved through dynamic programming. On the other hand, the Riemannian structure in moi17 is applied to curves in the hyperbolic plane, but the question of optimal matching is not studied. This is remedied here.

1.3 Contributions of this paper

The aim of this paper is three-fold : (1) study the quotient structure associated to the Riemannian metric $G$ of moi17 , and more generally to any elastic metric on manifold-valued curves, (2) exploit this knowledge of horizontality in an algorithm that computes an optimal matching between shapes – and whose range of application goes beyond elastic metrics – and (3) give a comprehensive discrete Riemannian framework on the finite-dimensional manifold of "discrete curves" that correctly approximates these procedures for the particular case of the SRV framework and constant sectional curvature. More specifically,

• $\bullet$

we characterize the horizontal subspace associated to the quotient shape space, for any elastic metric and in particular for the SRV metric $G$. Namely, we decompose any infinitesimal deformation of a smooth curve into a vertical part, which reparameterizes the curve without changing its shape, and a $G$-orthogonal horizontal part. This is done in a similar way as in tumpach16 but without restriction to arc-length parameterized curves.

• $\bullet$

We write any path in the space of smooth curves as a horizontal path composed with a path of reparameterizations. We use this decomposition to define an algorithm that, for a fixed parameterization of one of two curves, approximates the horizontal geodesic linking it to the fiber of the other curve, thereby yielding the "closest parameterization" of the latter with respect to the fixed parameterization of the former. We refer to this correspondence as an optimal matching. This algorithm can be used for any metric structure as long as one knows how to compute geodesics and characterize horizontality. Comparison with the popular dynamic programming approach srivastava16 is established in the simulations section.

• $\bullet$

We define a discrete version of $G$ that is a Riemannian metric on the finite-dimensional manifold $M^{n+1}$ of "discrete curves" given by $n+1$ points, when $M$ has constant sectional curvature $K\in\{-1,0,1\}$. We show that the energy of a path of discrete curves converges to the energy of the limit path as the size $n$ of the discretization goes to $\infty$. We give the geodesic equations for this metric, characterize the Jacobi fields, describe geodesic shooting, and show simulations on synthetic and real data in $\mathbb{R}^{2},\mathbb{R}^{3}$, the hyperbolic plane and the sphere.

1.4 Outline

After reminding the continuous model previously introduced in moi17 , Section 2 describes the horizontal space of the quotient structure and a way to compute horizontal geodesics. In Section 3, we introduce the discretization, and give the convergence result toward the continuous model, which is later proved in Section 5. Section 4 shows results of simulations in the three settings of positive, zero and negative curvature.

2 The continuous model

2.1 Notations

Let $\left(M,\langle\cdot,\cdot\rangle\right)$ be a Riemannian manifold. We first introduce a few notations. The norm associated to the Riemannian metric $\langle\cdot,\cdot\rangle$ is denoted by $|\cdot|$, the Levi-Civita connection by $\nabla$ and the curvature tensor by $\mathcal{R}$. If $t\mapsto c(t)$ is a curve in $M$ and $t\mapsto w(t)\in T_{c(t)}M$ a vector field along $c$, we denote by $c_{t}:=dc/dt=c^{\prime}$ the derivative of $c$ with respect to $t$ and by $\nabla_{t}w:=\nabla_{c_{t}}w$, $\nabla^{2}_{t}w:=\nabla_{c_{t}}\nabla_{c_{t}}w$ the first and second order covariant derivatives of $w$ along $c$. Parallel transport of a tangent vector $u\in T_{c(t_{1})}M$ from $c(t_{1})$ to $c(t_{2})$ along $c$ is denoted by $P_{c}^{t_{1},t_{2}}(u)$, or when there is no ambiguity on the choice of the curve $c$, $u^{t_{1},t_{2}}$, or even $u^{\parallel}$ to lighten notations in some cases. We associate to each curve $c$ its renormalized speed vector field $v:=c^{\prime}/|c^{\prime}|$, and to each vector field $t\mapsto w(t)$ along $c$, its tangential and normal components $w^{T}:=\langle w,v\rangle v$ and $w^{N}:=w-w^{T}$. Finally, for all $x\in M$ we denote by $\exp_{x}:T_{x}M\rightarrow M$ the exponential map on $M$ and by $\log_{x}:M\rightarrow T_{x}M$ its inverse map.

2.2 The space of smooth parameterized curves

2.2.1 The Riemannian structure

We represent open oriented curves in $M$ by smooth immersions, i.e. smooth curves with velocity that doesn’t vanish. The set $\mathcal{M}$ of smooth immersions in $M$ is an open submanifold of the Fréchet manifold $C^{\infty}([0,1],M)$ michor16 and its tangent space at a point $c$ is the set of infinitesimal deformations of $c$, which can be seen as vector fields along the curve $c$ in $M$

[TABLE]

Reparametrizations are represented by increasing diffeomorphisms $\varphi:[0,1]\rightarrow[0,1]$ (so that they preserve the end points of the curves), and their set is denoted by $\text{Diff}^{+}([0,1])$. We adopt the so-called square root velocity (SRV) representation, i.e. we represent each curve $c\in\mathcal{M}$ by the pair formed by its starting point $x$ and its speed vector field renormalized by the square root of its norm, via the bijection $\mathcal{M}\rightarrow M\times T\mathcal{M}$

[TABLE]

The inverse of this function is simply given by $M\times T\mathcal{M}\ni(x,q)\mapsto\pi_{\mathcal{M}}(q)\in\mathcal{M}$, if $\pi_{\mathcal{M}}$ is the canonical projection $T\mathcal{M}\rightarrow\mathcal{M}$. The renormalization of the speed vector field in $q$ allows us to define a reparameterization invariant metric, as we will see shortly. For any tangent vector $w\in T_{c}\mathcal{M}$, consider a path of curves $s\mapsto c^{w}(s)\in\mathcal{M}$ such that $c^{w}(0)=c$ and $c^{w}_{s}(0):=\partial c^{w}/\partial s(0)=w$. We denote by $q^{w}:=c^{w}_{t}/|c^{w}_{t}|^{1/2}$ the square root velocity representation of $c^{w}$. With these notations, we equip $\mathcal{M}$ with a Riemannian metric $G$, defined by

[TABLE]

This definition does not depend on the choice of $c^{w}$ and we can reformulate this scalar product in terms of (covariant) derivatives of $w$. Indeed, note that $\nabla_{s}q^{w}(0,t)=\nabla_{c^{w}_{s}(0,t)}(c^{w}_{t}/|c^{w}_{t}|^{\frac{1}{2}})=|c^{w}_{t}(0,t)|^{-\frac{1}{2}}(\nabla_{c^{w}_{s}(0,t)}c^{w}_{t}-1/2(\nabla_{c^{w}_{s}(0,t)}c^{w}_{t})^{T})$, which gives after inverting the derivatives according to $s$ and $t$,

[TABLE]

The scalar product can then be rewritten

[TABLE]

where $\mathrm{d}\ell=|c^{\prime}(t)|\mathrm{d}t$ and $\nabla_{\ell}=\frac{1}{|c^{\prime}(t)|}\nabla_{t}$ respectively denote integration and covariant derivation according to arc length. This metric belongs to the class of so-called elastic metrics parameterized by any $a,b\in\mathbb{R}$, which can be defined for manifold-valued curves as

[TABLE]

With formulation (2) it is clear that $G=G^{1,\frac{1}{2}}$ is invariant under the action of reparameterizing the curve and its tangent vectors by any increasing diffeomorphism $\varphi\in\text{Diff}^{+}([0,1])$,

[TABLE]

This reparameterization invariance property will allow us to induce a Riemmannian structure on the quotient space as we will see in Section 2.3.

2.2.2 Geodesics between parameterized curves

Two curves $c_{0},c_{1}\in\mathcal{M}$ can be compared using the geodesic distance induced by $G$, i.e. by computing the length of the shortest path of curves $[0,1]\ni s\mapsto c(s)\in\mathcal{M}$ from $c_{0}$ to $c_{1}$

[TABLE]

where the length of a path $c$ can be written in terms of its SRV representation $(x,q):[0,1]\rightarrow M\times T\mathcal{M}$

[TABLE]

as

[TABLE]

Note that here - and in all that follows - we indifferently use the notations $c(s,t)=c(s)(t)$, $q(s,t)=q(s)(t)$ for all $s,t\in[0,1]$. Now we recall a result shown in moi17 , which characterizes the geodesic paths of $\mathcal{M}$, i.e. those which achieve the infimum in (4), by searching for the critical points of the energy functional $E:\mathcal{M}\rightarrow\mathbb{R}_{+}$,

[TABLE]

Proposition 1 (Geodesic equations)

A geodesic path $[0,1]\ni s\mapsto c(s)\in\mathcal{M}$ for $G$, or more specifically its SRV representative $s\mapsto(x(s),q(s))$ (5), verifies the equations

[TABLE]

for all $s,t\in[0,1]$, where the curvature term $r(s,t)$ integrates the vector field $\mathcal{R}(q,\nabla_{s}q)c_{s}$ parallel transported back to $t$ along $c(s,\cdot)$,

[TABLE]

Remark 1

In the flat case $M=\mathbb{R}^{d}$, the curvature term $r$ vanishes and we obtain $\nabla_{s}x_{s}(s)=0$, $\nabla_{s}^{2}q(s,t)=0$ for all $s$ and $t$. This means that the geodesic between two curves $(x_{0},q_{0})$ and $(x_{1},q_{1})$ in the SRV coordinates is composed of a straight line $s\mapsto x(s)$ and an $L^{2}$-geodesic $s\mapsto q(s,\cdot)$. In other words, the geodesic links the starting points with a straight line and linearly interpolates between the renormalized speeds. Notice that if this linear interpolation goes through zero, the geodesic does not exist in $\text{Imm}([0,1],M)$. This can be avoided by considering a larger space of curves, such as the set of absolutely continuous curves lahiri15 .

A possibility to construct the geodesics of $\mathcal{M}$ is to use geodesic shooting. By solving the geodesic equations (7) we can construct the geodesic path starting from a given curve $c_{0}$ at a given speed $w\in T_{c_{0}}\mathcal{M}$ - this is the exponential map on $\mathcal{M}$. Given two curves $c_{0},c_{1}$, geodesic shooting allows us to iteratively find the appropriate initial speed $w$ which will make the geodesic land on $c_{1}$. The idea is to "shoot" from $c_{0}$ in a certain direction using the exponential map, estimate the gap between the end point of the obtained geodesic and the target point $c_{1}$, "bring back" this information at $c_{0}$ using a Jacobi field and finally use this information to correct the shooting direction.

Jacobi fields are vector fields that describe the way that geodesics spread out in the Riemannian manifold: for any geodesic $s\mapsto c(s)$ in $\mathcal{M}$ and Jacobi field $s\mapsto J(s)$ along $c$, there exists a family of geodesics $(-\delta,\delta)\ni a\mapsto c(a,\cdot)$ such that for all $s$, $c(0,s)=c(s)$ and

[TABLE]

At a given step of the geodesic shooting algorithm, we consider the Jacobi field that measures the difference between the geodesic obtained by shooting and the desired geodesic between $c_{0}$ and $c_{1}$ : it has initial value $J(0)=0$ since both have same starting point $c_{0}$, and its end value can be estimated by $J(1)=\log^{L2}_{c_{0}}c_{1}$, where $\log^{L2}_{c}$ denotes the inverse of the exponential map for the $L^{2}$-metric on $\mathcal{M}$, simply given by $\log_{c_{0}}^{L2}(c_{1})(t)=\log_{c_{0}(t)}(c_{1}(t))$ for $t\in[0,1]$. The shooting direction can then be corrected by the derivative $\nabla_{s}J(0)$ of the Jacobi field at the origin, as shown in Figure 2. For more details, we refer the reader to moi17 .

Algorithm 1** (Geodesic shooting)**

Input: $c_{0},c_{1}\in\mathcal{M}$.

Initialization: $w=\log_{c_{0}}^{L2}(c_{1})$.

Repeat until convergence :

compute the geodesic $s\mapsto c(s)$ starting from $c_{0}$ at speed $w$ by solving the geodesic equations (7), 2. 2.

evaluate the difference $j:=\log_{c(1)}^{L^{2}}(c_{1})$ between the target curve $c_{1}$ and the extremity $c(1)$ of the obtained geodesic, 3. 3.

compute the initial derivative $\nabla_{s}J(0)$ of the Jacobi field $s\mapsto J(s)$ along the path of curves $c$ verifying $J(0)=0$ and $J(1)=j$, 4. 4.

correct the shooting direction $w=w+\nabla_{s}J(0)$.

Output: geodesic $c$.

This algorithm requires the characterization of the Jacobi fields for $G$ on $\mathcal{M}$, and a way to deduce the initial derivative $\nabla_{s}J(0)$ of a Jacobi field from its initial and final values $J(0)$, $J(1)$. Concerning these two points, we refer the reader to moi17 : the Jacobi fields of $\mathcal{M}$ are shown to be solutions of a linear PDE, which can be solved to obtain the final value $J(1)$ of a Jacobi field $J$ along a path of curves $c$ knowing its initial conditions $J(0)$ and $\nabla_{s}J(0)$. If we consider only Jacobi fields with initial value $J(0)=0$, then the function $\varphi:T_{c(0)}\mathcal{M}\rightarrow T_{c(1)}\mathcal{M}$, $\nabla_{s}J(0)\mapsto J(1)$ is a linear bijection between two vector spaces and its inverse map can be computed by considering the image of a basis of $T_{c(0)}\mathcal{M}$. The equations characterizing the Jacobi fields in the discrete setting will be given in Section 3.

2.3 The space of unparameterized curves

2.3.1 The quotient structure

In order to compare curves regardless of parameterization, we consider the quotient $\mathcal{S}=\mathcal{M}/\text{Diff}^{+}([0,1])$ of the space of curves by the diffeomorphisms group. This quotient is not a manifold, as it has singularities, i.e. points with non trivial isotropy group michor16 . That is why from now on we assume that $\mathcal{M}$ denotes the set of free immersions, i.e. immersions on which the diffeomorphism group acts freely. That way, the space of curves $\mathcal{M}$ and the quotient shape space $\mathcal{S}$ are respectively the total and base spaces of a principal bundle, the fibers of which are the sets of all the curves that are identical modulo reparameterization, i.e. that project on the same "shape". We denote by $\pi:\mathcal{M}\rightarrow S$ the projection of the fiber bundle and by $\bar{c}:=\pi(c)\in\mathcal{S}$ the shape of a curve $c\in\mathcal{M}$. The tangent bundle can then be decomposed

[TABLE]

into a vertical subspace consisting of all vectors tangent to the fibers of $\mathcal{M}$ over $\mathcal{S}$, that is, those which have an action of reparameterizing the curve without changing its shape

[TABLE]

and a horizontal subspace defined as the orthogonal of the vertical subspace according to $G$

[TABLE]

If $G$ is constant along the fibers, i.e. verifies property (3), then there exists a Riemannian metric $\bar{G}$ on the shape space $\mathcal{S}$ such that $\pi$ is a Riemannian submersion from $(\mathcal{M},G)$ to $\mathcal{(}S,\bar{G})$,

[TABLE]

where $w^{hor}$ and $z^{hor}$ are the horizontal parts of tangent vectors $w$ and $z$, as well as the horizontal lifts of $T_{c}\pi(w)$ and $T_{c}\pi(z)$, respectively. This expression defines $\bar{G}$ in the sense that it does not depend on the choice of the representatives $c$, $w$ and $z$ (michor08 , §29.21). If a geodesic for $G$ has a horizontal initial speed, then its speed vector stays horizontal at all times - we say it is a horizontal geodesic - and projects on a geodesic of the shape space for $\bar{G}$ (michor08 , §26.12). To compute the distance between two shapes $\overline{c_{0}}$ and $\overline{c_{1}}$ in the quotient space we choose a representative $c_{0}$ of $\overline{c_{0}}$ and compute the distance (in $\mathcal{M}$) to the closest representative of $\overline{c_{1}}$, as shown schematically in Figure 3,

[TABLE]

By definition, the distance in the quotient space allows us to compare curves regardless of parameterization

[TABLE]

We now characterize the horizontal subspace for any elastic metric $G^{a,b}$ and in particular for our metric $G^{1,\frac{1}{2}}$, and give the decomposition of a tangent vector.

Proposition 2 (Horizontal part of a vector)

Let $c\in\mathcal{M}$ be a smooth immersion. Then $h\in T_{c}\mathcal{M}$ is horizontal for the elastic metric $G^{a,b}$ if and only if

[TABLE]

In particular, for $a=2b=1$ we obtain

[TABLE]

Any tangent vector $w\in T_{c}\mathcal{M}$ can be decomposed in horizontal and vertical components $w=w^{hor}+w^{ver}$ given by $\,\,w^{ver}=mv,\,\,w^{hor}=w-mv$, where the real function $m\in\mathcal{C}$ verifies the ordinary differential equation

[TABLE]

Proof

Let $h\in T_{c}\mathcal{M}$ be a tangent vector. It is horizontal if and only if it is orthogonal to any vertical vector, that is any vector of the form $mv$ with $m\in C^{\infty}([0,1],\mathbb{R})$ such that $m(0)=m(1)=0$. We have $\nabla_{t}(mv)=m^{\prime}v+m\nabla_{t}v\,$ and since $\langle\nabla_{t}v,v\rangle=0$ we get $\nabla_{t}(mv)^{N}=m\nabla_{t}v$ and $\nabla_{t}(mv)^{T}=m^{\prime}v$. The scalar product can then be written

[TABLE]

by integration by parts. The vector $h$ is horizontal if and only if $G_{c}(h,mv)=0$ for all such $m$, and so multiplying by $|c^{\prime}|/b^{2}$ gives the desired equation. Now consider a tangent vector $w$ and a real function $m:[0,1]\rightarrow\mathbb{R}$ such that $m(0)=m(1)=0$. Then according to the above, $w-mv$ is horizontal if and only if it verifies

[TABLE]

i.e., since $\langle\nabla_{t}v,v\rangle=0$, $\langle\nabla_{t}^{2}v,v\rangle=-|\nabla_{t}v|^{2}$ and $\nabla_{t}^{2}(mv)=m^{\prime\prime}v+2m^{\prime}\nabla_{t}v+m\nabla_{t}^{2}v$, if

[TABLE]

which is what we wanted.

2.3.2 Computing geodesics in the shape space

Recall that the geodesic path $s\mapsto\bar{c}(s)$ between the shapes of two curves $c_{0}$ and $c_{1}$ is the projection of the horizontal geodesic - if it exists - $s\mapsto c_{h}(s)$ linking $c_{0}$ to the fiber of $c_{1}$ in $\mathcal{M}$, i.e. such that $c_{h}(0)=c_{0}$, $c_{h}(1)\in\pi^{-1}(\overline{c_{1}})$ and $\partial_{s}c_{h}(s)\in\text{Hor}_{c_{h}(s)}$ for all $s\in[0,1]$,

[TABLE]

The end point of $c_{h}$ then gives the optimal reparameterization $c_{1}\circ\varphi$ of the target curve $c_{1}$ with respect to the initial curve $c_{0}$, i.e. such that

[TABLE]

The parameterized curves $(c_{0},c_{1}\circ\varphi)$ define what we call an optimal matching between the shapes $\overline{c_{0}}$ and $\overline{c_{1}}$, in the sense that each point $c_{0}(t)$ on $\overline{c_{0}}$ is matched with the point $c_{1}(\varphi(t))$ on $\overline{c_{1}}$. Here we propose a method to approach the horizontal geodesic $c_{h}$, and thereby the corresponding optimal matching. To that end we decompose any path of curves $s\mapsto c(s)$ in $\mathcal{M}$ into a horizontal path composed with a path of reparameterizations, $c(s)=c^{hor}(s)\circ\varphi(s)$, or equivalently

[TABLE]

where the path $[0,1]\ni s\mapsto c^{hor}(s)\in\mathcal{M}$ is such that $c^{hor}_{s}(s)\in\text{Hor}_{c^{hor}(s)}$ for all $s\in[0,1]$, and $[0,1]\ni s\mapsto\varphi(s)\in\text{Diff}^{+}([0,1])$ is a path of increasing diffeomorphisms. The horizontal and vertical parts of the speed vector of $c$ can be expressed in terms of this decomposition. Indeed, by taking the derivative of (9) with respect to $s$ and $t$ we obtain

[TABLE]

and so with $v^{hor}(s,t):=c^{hor}_{t}(s,t)/|c^{hor}_{t}(s,t)|$, since $\varphi_{t}>0$, (10b) gives

[TABLE]

We can see that the first term on the right-hand side of Equation (10a) is horizontal. Indeed, for any path of real functions $m:[0,1]\rightarrow C^{\infty}([0,1],\mathbb{R})$, $s\mapsto m(s,\cdot)$ such that $m(s,0)=m(s,1)=0$ for all $s$, since $G$ is reparameterization invariant we have

[TABLE]

with $\tilde{m}(s)=m(s)\circ\varphi(s)^{-1}$. Since $\tilde{m}(s,0)=\tilde{m}(s,1)=0$ for all $s$, the vector $\tilde{m}(s)\cdot v^{hor}(s)$ is vertical and its scalar product with the horizontal vector $c^{hor}_{s}(s)$ vanishes. On the other hand, the second term on the right hand-side of Equation (10a) is vertical, since it can be written

[TABLE]

with $m(s)=|c_{t}(s)|\varphi_{s}(s)/\varphi_{t}(s)$ verifying $m(s,0)=m(s,1)=0$ for all $s$. Finally, the vertical and horizontal parts of the speed vector $c_{s}(s)$ are given by

[TABLE]

Definition 1

We call $c^{hor}$ the horizontal part of the path $c$ with respect to $G$.

Proposition 3

The horizontal part of a path of curves $c$ is at most the same length as $c$

[TABLE]

Proof

Since the metric $G$ is reparameterization invariant, the squared norm of the speed vector of the path $c$ at time $s\in[0,1]$ is given by

[TABLE]

where $\|\cdot\|^{2}_{G}:=G(\cdot,\cdot)$. This gives $\|c^{hor}_{s}(s)\|_{G}\leq\|c_{s}(s)\|$ for all $s$ and so $L_{G}(c^{hor})\leq L_{G}(c)$.

Now we will see how the horizontal part of a path of curves can be computed.

Proposition 4 (Horizontal part of a path)

Let $s\mapsto c(s)$ be a path in $\mathcal{M}$. Then its horizontal part is given by $c^{hor}(s,t)=c(s,\varphi(s)^{-1}(t))$, where the path of diffeomorphisms $s\mapsto\varphi(s)$ is solution of the partial differential equation

[TABLE]

and where $m(s):[0,1]\rightarrow\mathbb{R}$, $t\mapsto m(s,t)$ is solution for all $s$ of the ordinary differential equation

[TABLE]

Proof

We have seen in Equation (11a) that the vertical part of $c_{s}(s)$ can be written as $m(s)\cdot v(s)$ where $m(s)=|c_{t}(s)|\varphi_{s}(s)/\varphi_{t}(s)$, and as the norm of the vertical part of $c_{s}(s)$, $m(s)$ is solution of the ODE (8) for all $s$.

If we take the horizontal part of the geodesic linking two curves $c_{0}$ and $c_{1}$, we will obtain a horizontal path linking $c_{0}$ to the fiber of $c_{1}$ which will no longer be a geodesic path. However this path reduces the distance between $c_{0}$ and the fiber of $c_{1}$, and gives a "better" representative $c_{1}\circ\varphi(1)$ of the target curve. By computing the geodesic between $c_{0}$ and this new representative, we are guaranteed to reduce once more the distance to the fiber. The algorithm that we propose simply iterates these two steps, as illustrated in Figure 4.

Algorithm 2** (Optimal matching)**

Input: $c_{0},c_{1}\!\in\mathcal{M}$.

Initialization: $n=0$, $c_{1}^{n}=c_{1}$.

Repeat until convergence:

construct the geodesic $s\mapsto c^{n}(s)$ between $c_{0}$ and $c_{1}^{n}$, 2. 2.

compute the horizontal part $s\mapsto[c^{n}]^{hor}(s)$ of $c^{n}$, 3. 3.

set $c^{n+1}_{1}=[c^{n}]^{hor}(1)$, 4. 4.

set $n=n+1$.

Output: horizontal geodesic $c^{n}$.

Let us specify why the obtained geodesic is horizontal at the limit. The series of lengths $\left(L(c^{n})\right)_{n\geq 0}$ and $\left(L([c^{n}]^{hor})\right)_{n\geq 0}$ are non negative decreasing and verify at each step $n$

[TABLE]

which means that they converge to the same limit. The same is true for the energies $E(c^{n})$ and $E([c^{n}]^{hor})$, and since the $s$-derivative of $[c^{n}]^{hor}$ is equal to the horizontal part of the $s$-derivative of $c^{n}$, we get

[TABLE]

and so $\|c^{n}_{s}(s)^{ver}\|_{G}$ converges to zero for almost all $s$. Since it is enough that a geodesic be horizontal at one given time for it to be horizontal for all time $s$ (michor08 , §26.12), we have the following result.

Proposition 5 (Horizontality of the solution)

At the limit, the geodesic between the fibers computed in Algorithm 2 is horizontal

[TABLE]

Remark 2

In this work, we will carry out step 1 using geodesic shooting. However it is important to stress that Algorithm 2 is a general method that can be applied to any metric structure (not only elastic metrics) for which one knows how to compute geodesics and characterize the horizontal subspace of the shape bundle. It can be seen as an alternative method to the popular dynamic programming approach srivastava16 , with which we establish comparisons in Section 4. But before that, let us first introduce a formal discretization of the continuous model presented so far.

3 The discrete model

Applications usually give access to a finite number of observations of a continuous process and provide series of points instead of continuous curves. It is therefore important to discretize the framework presented above and to consider the finite-dimensional space of "discrete curves". From now on we restrict to base manifolds $M$ of constant sectional curvature $K$. This allows us to get an explicit formula for the Jacobi fields of $M$ (Lemma 1) and thus derive a precise approximation of the smooth Riemannian structure of Section 2. Generalization to any Riemannian manifold for which the Jacobi fields are computable should not be problematic. For more complex manifolds, a faster and more approximate solution would be to directly discretize the smooth equations, at the cost of the precision of the discrete approximation.

3.1 The Riemannian structure

We consider the product manifold $M^{n+1}$ of "discrete curves" given by $n+1$ points, for a fixed $n\in\mathbb{N}^{*}$. Its tangent space at a given point $\alpha=(x_{0},\ldots,x_{n})$ is

[TABLE]

Assuming that there exists a connecting geodesic between $x_{k}$ and $x_{k+1}$ for all $k$ – which seems reasonable considering that the points $x_{k}$ should be "close" since they correspond to the discretization of a continuous curve – and that two consecutive points are always distinct ($x_{k}\neq x_{k+1}$), we use the following notations

[TABLE]

as well as ${w_{k}}^{T}=\langle w_{k},v_{k}\rangle v_{k}$ and $w_{k}^{N}=w_{k}-{w_{k}}^{T}$ to refer to the tangential and normal components of a tangent vector $w_{k}\in T_{x_{k}}M$. Given a tangent vector $w\in T_{\alpha}M^{n+1}$, we consider a path of piecewise geodesic curves $[0,1]^{2}\ni(s,t)\mapsto c^{w}(s,t)\in M$ such that $c^{w}(0,\frac{k}{n})=x_{k}$ for $k=0,\ldots,n$, $c^{w}(s,\cdot)$ is a geodesic of $M$ on the interval $[\frac{k}{n},\frac{k+1}{n}]$ for all $s\in[0,1]$ and $k$ – and in particular $c^{w}_{t}(0,\frac{k}{n})=n\tau_{k}$ – and such that $c^{w}_{s}(0,\frac{k}{n})=w_{k}$, as illustrated in Figure 5. Then we define the squared norm of $w$ by

[TABLE]

This definition is a discrete analog of (1), and just as in the continuous case, it does not depend on the choice of $c^{w}$. Indeed, we can also obtain a discrete analog of (2).

Proposition 6

The metric $G^{n}$ can also be written

[TABLE]

where the map $D_{\tau}:T_{\alpha}M^{n+1}\rightarrow T_{\alpha}M^{n+1}$, w\mapsto D_{\tau}w=\big{(}(D_{\tau}w)_{0},\ldots,(D_{\tau}w)_{n}\big{)} is defined by

[TABLE]

if $w_{k+1}^{\parallel}$ denotes the parallel transport of $w_{k+1}$ from $x_{k+1}$ to $x_{k}$ along the geodesic, and the coefficients $a_{k}$ and $b_{k}$ take the following values depending on the sectional curvature $K$ of the base manifold $M$

[TABLE]

Remark 3

Notice that in the flat case our definition gives $(D_{\tau}w)_{k}=w_{k+1}-w_{k}$. In the non-flat case, when the discretization gets "thinner", i.e. $n\rightarrow\infty$ and $|\tau_{k}|\rightarrow 0$ while $n|\tau_{k}|$ stays bounded for all $0\leq k\leq n$, we get $(D_{\tau}w)_{k}\underset{n\rightarrow\infty}{=}{w_{k+1}}^{\parallel}-w_{k}+o(1)$.

Before we prove this proposition, let us recall a well-known result about Jacobi fields that will prove useful to derive the equations in the discrete case.

Lemma 1

Let $\gamma:[0,1]\rightarrow M$ be a geodesic of a manifold $M$ of constant sectional curvature $K$, and $J$ a Jacobi field along $\gamma$. Then the parallel transport of $J(t)$ along $\gamma$ from $\gamma(t)$ to $\gamma(0)$ is given by

[TABLE]

for all $t\in[0,1]$, where

[TABLE]

Proof (Proof of Lemma 1)

For the sake of completeness, the proof is reminded in the appendix.

Proof (Proof of Proposition 6)

Let $\alpha\in M^{n+1}$ be a "discrete curve" and $w\in T_{\alpha}M^{n+1}$ a tangent vector at $\alpha$. Consider a path of piecewise geodesic curves $s\mapsto c^{w}(s)$ that verifies all the conditions given above to define $G^{n}(w,w)$, and set $(D_{\tau}w)_{k}:=\frac{1}{n}\nabla_{t}c^{w}_{s}(0,\frac{k}{n})$. Then by definition, the vector field $J_{k}(u)=c^{w}_{s}(s,\frac{k+u}{n})$, $u\in[0,1]$ is a Jacobi field along the geodesic linking $x_{k}$ to $x_{k+1}$, verifying $J_{k}(0)=w_{k}$, $J_{k}(1)=w_{k+1}$ and $\nabla_{u}J_{k}(0)=(D_{\tau}w)_{k}$. Applying Lemma 1 gives

[TABLE]

This gives $({w_{k+1}}^{\parallel})^{T}={w_{k}}^{T}+{(D_{\tau}w)_{k}}^{T}$ and $({w_{k+1}}^{\parallel})^{N}=a_{k}{w_{k}}^{N}+b_{k}{(D_{\tau}w)_{k}}^{N}$ and so

[TABLE]

Finally, we observe that the covariant derivative involved in the definition of $G^{n}$ can be written

[TABLE]

i.e.

[TABLE]

Injecting this into (14) gives the desired formula.

Now we present the main result of this section, that is, the convergence of the discrete model toward the continuous model.

Definition 2

Let $\alpha=(x_{0},\ldots,x_{n})\in M^{n+1}$ be a discrete curve, and $t\mapsto c(t)\in M$ a smooth curve. We say that $\alpha$ is the discretization of size $n$ of $c$ when $c(\frac{k}{n})=x_{k}$ for all $k=0,\ldots,n$. If $s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))\in M^{n+1}$ is a path of discrete curves and $s\mapsto c(s)\in\mathcal{M}$ a path of smooth curves, then $\alpha$ is the discretization of size $n$ of $c$ when $\alpha(s)$ is the discretization of $c(s)$ for all $s\in[0,1]$, i.e. when $x_{k}(s)=c(s,\frac{k}{n})$ for all $s$ and $k$. We will still use this term if $c$ is not smooth, and speak of the only path of piecewise-geodesic curves of which $\alpha$ is the discretization.

Let $[0,1]\ni s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))\in M^{n+1}$ be a path of discrete curves. Defining $\tau_{k}(s)$ and $q_{k}(s)$ as in (13) for all $s\in[0,1]$, the path $\alpha$ can be represented by its SRV representation $[0,1]\rightarrow M\times T_{\alpha}{M^{n+1}}$,

[TABLE]

To compute the squared norm of its speed vector $\alpha^{\prime}(s)$, consider the path of piecewise geodesic curves $[0,1]^{2}\ni(s,t)\mapsto c(s,t)\in M$ such that $c(s,\frac{k}{n})=x_{k}(s)$ and $c_{t}(s,\frac{k}{n})=n\tau_{k}(s)$ for all $s$ and $k$. Then, notice that we have

[TABLE]

and so the squared norm of the speed vector of $\alpha$ can be expressed in terms of the SRV representation

[TABLE]

In the following result, we show that if $s\mapsto\alpha(s)$ is the discretization of a path $s\mapsto c(s)\in\mathcal{M}$ of continuous curves, then its energy with respect to $G^{n}$,

[TABLE]

gets closer to the energy (6) of $c$ with respect to $G$ as the size of the discretization grows.

Theorem 3.1 (Convergence of the discrete model to the continuous model)

Let $s\mapsto c(s)$ be a $C^{1}$-path of $C^{2}$-curves with non vanishing derivative with respect to $t$. This path can be identified with an element $(s,t)\mapsto c(s,t)$ of $C^{1,2}([0,1]\times[0,1],M)$ such that $c_{t}\neq 0$. Consider the $C^{1}$-path in $M^{n+1}$, $s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))$, that is the discretization of size $n$ of $c$. Then there exists a constant $\lambda>0$ such that for $n$ large enough, the difference between the energies of $c$ and $\alpha$ is bounded by

[TABLE]

where $E$ and $E^{n}$ are the energies with respect to metrics $G$ and $G^{n}$ respectively and where

[TABLE]

if $|w|_{\infty}\!:=\!\!\!\!\underset{s,t\in[0,1]}{\sup}|w(s,t)|$ denotes the supremum over both $s$ and $t$ of a vector field $w$ along $c$.

Remark 4

Note that since we assume that $c$ is a $C^{1}$-path of $C^{2}$-curves, the following norms are bounded for $i=1,2$,

[TABLE]

Proof (Proof of Theorem 3.1)

The proof is put off to Section 5.

Now that we have established a formal Riemannian setting to study discrete curves defined by a series of points, and that we have studied its link to the continuous model, we need to derive the equations of the corresponding geodesics and Jacobi fields to apply the methods described in Section 2. For the sake of readability, we first introduce some notations.

3.2 Computing geodesics in the discrete setting

3.2.1 Notations

The purpose of the notations that we introduce here is to lighten the equations derived in the rest of the paper. For any discrete curve $\alpha=(x_{0},\ldots,x_{n})\in M^{n+1}$ we define for all $0\leq k\leq n$, using the coefficients $a_{k}$ and $b_{k}$ defined by (15) and (13), the functions $f_{k},g_{k}:T_{x_{k}}M\rightarrow T_{x_{k}}M$,

[TABLE]

and for $k=0,\ldots,n-1$, the functions $f_{k}^{(-)},g_{k}^{(-)}:T_{x_{k+1}}M\rightarrow T_{x_{k}}M$ by

[TABLE]

where $\gamma_{k}$ denotes the geodesic between $x_{k}$ and $x_{k+1}$, which we previously assumed existed. Notice that when the discretization gets "thinner", that is $n\rightarrow\infty$, $|\tau_{k}|\rightarrow 0$ while $n|\tau_{k}|$ stays bounded for all $0\leq k\leq n$, we get in the non flat setting, for any fixed $w\in T_{x_{k+1}}M$, $f_{k}(w)=w+o(1/n)$ and $g_{k}(w)=|q_{k}|(w+{w}^{T})+o(1/n)$ - in the flat setting, these are always equalities. Now if we consider a path $s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))$ of discrete curves, we can define for each $s$ the functions

[TABLE]

for $0\leq k\leq n$ and $0\leq k\leq n-1$ respectively, corresponding to the discrete curve $\alpha(s)$. It is of interest for the rest of this paper to compute the covariant derivatives of these maps with respect to $s$.

Lemma 2

The first and second order covariant derivatives of $f_{k}$ and $g_{k}$ with respect to $s$ are functions $T_{x_{k}(s)}M\rightarrow T_{x_{k}(s)}M$ defined by

[TABLE]

Proof

For any vector field $s\mapsto w(s)\in T_{x_{k}(s)}M$ along $s\mapsto x_{k}(s)$ we have by definition

[TABLE]

Noticing that $\nabla_{s}({w}^{T})=(\nabla_{s}w)^{T}+\langle w,\nabla_{s}v_{k}\rangle v_{k}+\langle w,v_{k}\rangle\nabla_{s}v_{k}$ and $\nabla_{s}(w^{N})=\nabla_{s}w-\nabla_{s}({w}^{T})$, the formulas given in Lemma 2 result from simple calculation.

Using these functions, we can deduce the covariant derivatives of $f_{k}^{(-)}$ and $g_{k}^{(-)}$. Denoting by $\gamma_{k}(s)$ the geodesic of $M$ linking $x_{k}(s)$ to $x_{k+1}(s)$ for all $s\in[0,1]$ and $0\leq k\leq n-1$, we have the following result.

Lemma 3

The covariant derivatives of the functions $f_{k}^{(-)}$ and $g_{k}^{(-)}$ with respect to $s$ are functions $T_{x_{k+1}(s)}M\rightarrow T_{x_{k}(s)}M$ given by

[TABLE]

where

[TABLE]

if $K$ is the sectional curvature of the base manifold.

Proof

The proof is given in Appendix A.

3.2.2 Geodesic equations and exponential map

With these notations, we can characterize the geodesics for metric $G^{n}$. The geodesic equations can be derived in a similar way as in the continuous case, that is by searching for the critical points of the energy (18). We obtain the following characterization in terms of the SRV representation (16).

Proposition 7 (Discrete geodesic equations)

A path $s\mapsto\alpha(s)=\left(x_{0}(s),\ldots,x_{n}(s)\right)$ in $M^{n+1}$ is a geodesic for metric $G^{n}$ if and only if its SRV representation s\mapsto\big{(}x_{0}(s),(q_{k}(s))_{k}\big{)} verifies the following differential equations

[TABLE]

for all $k=0,\ldots,n-1$, with the notations (13) and $R_{k}:=\mathcal{R}(q_{k},\nabla_{s}q_{k}){x_{k}}^{\prime}$.

Proof

The proof is given in Appendix B.

Remark 5 (Link with the continuous setting)

Let $[0,1]\ni s\mapsto c(s,\cdot)\in\mathcal{M}$ be a $C^{1}$ path of smooth curves and $[0,1]\ni s\mapsto\alpha(s)\in M^{n+1}$ the discretization of size $n$ of $c$. We denote as usual by $q:=c_{t}/|c_{t}|^{1/2}$ and $(q_{k})_{k}$ their respective SRV representations. When $n\rightarrow\infty$ and $|\tau_{k}|\rightarrow 0$ while $n|\tau_{k}|$ stays bounded for all $0\leq k\leq n$, the coefficients of the discrete geodesic equation (20) for $\alpha$ converge to the coefficients of the continuous geodesic equation (7) for $c$, i.e.

[TABLE]

for all $s\in[0,1]$ and $k=0,\ldots,n-1$, where $r_{n-1}=0$ and for $k=1,\ldots,n-2$,

[TABLE]

with the exception that the sum starts at $\ell=0$ for $r_{0}$. More details on this can be found in Appendix B.

Remark 6 (Euclidean case and existence of geodesics)

Just as in the continuous case, when $M=\mathbb{R}^{d}$, the curvature terms $R_{k}$’s vanish and we obtain

[TABLE]

i.e. the geodesics are composed of straight lines in the SRV coordinates. We can again avoid the problem of the $q_{k}(s)$’s going through zero by allowing two consecutive components $x_{k}$ and $x_{k+1}$ to be equal, and setting $q_{k}=0$ when that happens. In that case, we get a complete finite-dimensional manifold, which is by the Hopf-Rinow theorem geodesically complete, i.e. any two curves can be linked by a minimizing geodesic. Indeed, since the SRV coordinates of geodesics are straight lines, a sequence in $(\mathbb{R}^{d})^{n+1}$ converges if and only if its SRV coordinates $x_{0},q_{0},\ldots,q_{n-1}$ converge in $\mathbb{R}^{d}$ (the sequence subscript is omited), and so the completeness of $(\mathbb{R}^{d})^{n+1}$ follows from that of $\mathbb{R}^{d}$. The question of whether this property still holds in the non flat case is postponed to future work.

Using equations (20) we can now build the exponential map, that is, an algorithm allowing us to approximate the geodesic of $M^{n+1}$ starting from a point $(x^{0}_{0},\ldots,x^{0}_{n})\in M^{n+1}$ at speed $(u_{0},\ldots,u_{n})$ with $u_{k}\in T_{x_{k}}M$ for all $k=0,\ldots,n$. In other words, we are looking for a path $[0,1]\ni s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))$ such that $x_{k}(0)=x^{0}_{k}$ and ${x_{k}}^{\prime}(0)=u_{k}$ for all $k$, and that verifies the geodesic equations (20). Assume that we know at time $s\in[0,1]$ the values of $x_{k}(s)$ and ${x_{k}}^{\prime}(s)$ for all $k=0,\ldots,n$. Then we propagate using

[TABLE]

In the following proposition, we see how we can compute the acceleration $\nabla_{s}{x_{k}}^{\prime}$ for each $k$.

Proposition 8 (Discrete exponential map)

Let $[0,1]\ni s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))$ be a geodesic path in $M^{n+1}$. For all $s\in[0,1]$, the coordinates of its acceleration $\nabla_{s}\alpha^{\prime}(s)$ can be iteratively computed in the following way

[TABLE]

for $k=0,\ldots,n-1$, where the $R_{k}$’s are defined as in Proposition 7, the symbol $\cdot^{\parallel}$ denotes the parallel transport from $x_{k+1}(s)$ back to $x_{k}(s)$ along the geodesic linking them, the maps $\nabla_{s}f_{k}$ and $\nabla_{s}g_{k}$ are given by Lemma 2, $Y_{k}$ is given by Equation (Appendix A) and

[TABLE]

Proof

The proof is given in Appendix B.

The equations of Proposition 8 allow us to iteratively construct a geodesic $s\mapsto\alpha(s)$ in $M^{n+1}$ for metric $G^{n}$ from the knowledge of its initial conditions $\alpha(0)$ and $\alpha^{\prime}(0)$. The next step is to construct geodesics under boundary constraints, i.e. to find the shortest path between two elements $\alpha_{0}$ and $\alpha_{1}$ of $M^{n+1}$.

3.2.3 Jacobi fields and geodesic shooting

As explained in Section 2.2.2 for the continuous model, we solve the boundary value problem using geodesic shooting. To do so, recall that we need to characterize the Jacobi fields for the metric $G^{n}$, since these play a role in the correction of the shooting direction at each iteration of the algorithm. Recall also that for any geodesic $s\mapsto\alpha(s)$ in $M^{n+1}$ and Jacobi field $s\mapsto J(s)$ along $\alpha$, there exists a family of geodesics $(-\delta,\delta)\ni a\mapsto\alpha(a,\cdot)$ such that $\alpha(0,s)=\alpha(s)$ for all $s$ and

[TABLE]

Proposition 9 (Discrete Jacobi fields)

Let $[0,1]\ni s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))$ be a geodesic path in $M^{n+1}$, $[0,1]\ni s\mapsto J(s)=(J_{0}(s),\ldots,J_{n}(s))$ a Jacobi field along $\alpha$, and $(-\delta,\delta)\ni a\mapsto\alpha(a,\cdot)$ a corresponding family of geodesics, in the sense just described. Then $J$ verifies the second order linear ODE

[TABLE]

for all $0\leq k\leq n-1$, where $R_{k}:=\mathcal{R}(q_{k},\nabla_{s}q_{k}){x_{k}}^{\prime}$ and the various covariant derivatives according to $a$ can be expressed as functions of $J$ and $\nabla_{s}J$,

[TABLE]

with the notation conventions $f_{k+1}^{(-)}\circ\ldots\circ f_{k-1}^{(-)}:=Id$, $\sum_{\ell=n}^{n-1}:=0$ and with the maps

[TABLE]

and

[TABLE]

Proof

The proof is given in Appendix B.

The equations of Proposition 9 allow us to iteratively compute the Jacobi field $J$ along a geodesic $\alpha$ - and in particular, its end value $J(1)$ - from the knowledge of the initial conditions $\{J_{k}(0),0\leq k\leq n\}$ and $\{\nabla_{s}J_{k}(0),0\leq k\leq n\}$. Indeed, if at time $s\in[0,1]$ we have $J_{k}(s)$ and $\nabla_{s}J_{k}(s)$ for all $k=0,\ldots,n$, then we can propagate using

[TABLE]

where $\nabla_{s}^{2}J_{k}(s)$ is deduced from $\nabla_{s}^{2}J_{k-1}(s)$ using Proposition 9. We can now apply Algorithm 1, where we replace the smooth geodesic equations (7) by the discrete geodesic equations (20) and we solve them using the exponential map described in Proposition 8. Notice that in $M^{n+1}$, the $k^{th}$ component of the $L^{2}$-logarithm map between two elements $\alpha_{0}=(x_{0}^{0},\ldots,x_{n}^{0})$ and $\alpha_{1}=(x_{0}^{1},\ldots,x_{n}^{1})$ is given by $\log_{x^{0}_{k}}(x^{1}_{k})$.

Algorithm 3** (Discrete geodesic shooting)**

Input: $\alpha_{0}=(x^{0}_{0},\ldots,x^{0}_{n}),\alpha_{1}=(x^{1}_{0},\ldots,x^{1}_{n})$.

Initialization: $w=\log_{\alpha_{0}}^{L^{2}}(\alpha_{1})$.

Repeat until convergence :

compute the geodesic $s\mapsto\alpha(s)$ starting from $\alpha_{0}$ at speed $w$ using Proposition 8, 2. 2.

evaluate the difference $j:=\log_{\alpha(1)}^{L^{2}}(\alpha_{1})$ between the target curve $\alpha_{1}$ and the extremity $\alpha(1)$ of the obtained geodesic, 3. 3.

compute the initial derivative $\nabla_{s}J(0)$ of the Jacobi field $s\mapsto J(s)$ along $\alpha$ verifying $J(0)=0$ and $J(1)=j$ using Proposition 9, 4. 4.

correct the shooting direction $w=w+\nabla_{s}J(0)$.

Output : geodesic $\alpha(s)$.

Recall that the map $\varphi:T_{\alpha(0)}\mathcal{M}\rightarrow T_{\alpha(1)}\mathcal{M}$, $\nabla_{s}J(0)\mapsto J(1)$ associating to the initial derivative $\nabla_{s}J(0)$ of a Jacobi field with initial value $J(0)=0$ its end value $J(1)$, is a linear bijection between two vector spaces which can be obtained using Proposition 9. Its inverse map can be computed by considering the image of a basis of $T_{c(0)}\mathcal{M}$.

3.3 A discrete analog of unparameterized curves

The final step in building our discrete model is to introduce a discretization of the quotient shape space. There seems to be no natural, intrinsic definition of the shape of a discrete curve, as by definition we are lacking information : we only have access to a finite number $n+1$ of points. Therefore we will make the assumption that we know the equations of the underlying curves, that is, that for each discrete curve $\alpha$, we have access to the shape $\bar{\alpha}$ of the smooth curve of which $\alpha$ is the discretization. In practice, we can set $\bar{\alpha}$ to be the shape of an optimal interpolation. The goal, for two elements $\alpha_{0},\alpha_{1}$ of shapes $\overline{\alpha_{0}},\overline{\alpha_{1}}$, is to redistribute the $n+1$ points on $\overline{\alpha_{1}}$ to minimize the discrete distance to the $n+1$ points $\alpha_{0}$ on $\overline{\alpha_{0}}$, and obtain

[TABLE]

where $d_{n}$ is the geodesic distance associated to the discrete metric $G^{n}$. We approximate $\alpha_{1}^{opt}$ using Algorithm 2, i.e. by iteratively computing the "horizontal part" of the geodesic linking $\alpha_{0}$ to an iteratively improved discretization of $\overline{\alpha_{1}}$. Since there is no "discrete shape bundle", we simply define the vertical and horizontal spaces in $\alpha$ as the discrete analogs of the ones of the smooth case

[TABLE]

where $v=(v_{k})_{k}$ is still defined by (13). Similarly to the continuous case, we can show the following result.

Proposition 10 (Discrete horizontal space)

Let $\alpha\in M^{n+1}$ and $h\in T_{\alpha}M^{n+1}$. Then $h\in\text{Hor}^{n}_{\alpha}$ if and only if it verifies

[TABLE]

with the notation $\lambda_{k}:=\langle v_{k+1}^{\parallel},v_{k}\rangle$. Any tangent vector $w\in T_{\alpha}M^{n+1}$ can be uniquely decomposed into a sum $w=w^{ver}+w^{hor}$ where $w^{ver}=mv\in\text{Ver}_{\alpha}^{n}$, $w^{hor}=w-mv\in\text{Hor}_{\alpha}^{n}$ and the components $(m_{k})_{k}$ verify $m_{0}=m_{1}=0$ and the following recurrence relation

[TABLE]

with coefficients

[TABLE]

[TABLE]

Proof

Let $h\in T_{\alpha}\mathcal{M}$ be a tangent vector. It is horizontal if and only if it is orthogonal to any vertical vector, that is any vector of the form $mv$ with $m=(m_{k})_{k}\in\mathbb{R}^{n+1}$ such that $m_{0}=m_{n}=0$. Recall that by definition

[TABLE]

and so with the notation $\lambda_{k}:=\langle{v_{k+1}}^{\parallel},v_{k}\rangle$, we get

[TABLE]

The scalar product between $h$ and $mv$ is then

[TABLE]

Changing the indices in the first sum and taking into account that $m_{0}=m_{n}=0$, we obtain

[TABLE]

Since this is true for all such $m$ the summand is equal to zero for all $k$ and we get the desired equation. The decomposition of a tangent vector $w$ into a vertical part $mv$ and a horizontal part $w-mv$ with $m=(m_{k})_{k}\in\mathbb{R}^{n+1}$ such that $m_{0}=m_{n}=0$, is then simply characterized by the fact that $w-mv$ verifies this equation.

In the discrete case, soving the ODE (8) to find the horizontal part of a vector simply boils down to solving the recurrence relation (22), allowing us to compute the coefficients of the PDE (12). Now we present an algorithm to solve a discrete version of this PDE and compute the discrete analog of the horizontal part of a (discrete) path of discrete curves.

Algorithm 4** (Horizontal part of a path)**

Input : $\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))$, $s=0,\frac{1}{m},\ldots,1$.

Initialization : for $k=0,\ldots,n$,

[TABLE]

For $j=0,\ldots,m-1$,

set $s=j/m$, $m_{0}(s)=m_{n+1}(s)=0$ and solve

[TABLE] 2. 2.

for $k=0,\ldots,n$,

[TABLE] 3. 3.

interpolate between the $\{x_{k}(s+\tfrac{1}{m}),k=0,\ldots,n\}$ to obtain $N+1$ values $\{y_{\ell}(s+\tfrac{1}{m}),\ell=0,\ldots,N\}$ and interpolate between the $\{\varphi(s+\tfrac{1}{m})(\tfrac{k}{n}),k=0,\ldots,n\}$ to obtain $N+1$ values $\{\psi(s+\tfrac{1}{m})(\tfrac{\ell}{n}),\ell=0,\ldots,N\}$, 4. 4.

for $k=0,\ldots,n$,

[TABLE]

Output: $\alpha^{hor}(s)=(x_{0}^{hor}(s),..,x_{n}^{hor}(s)),s=0,\frac{1}{m},\ldots,1$.

Step 1 computes the coefficients of the PDE. The $A_{k},B_{k},C_{k},D_{k}$ are computed using the definitions of Proposition 10 for $w=\alpha^{\prime}(s)$, i.e. $(D_{\tau}\alpha^{\prime}(s))_{k}=\nabla_{s}\tau_{k}(s)$. Step 2 solves the PDE. The increment $\Delta\varphi$ is a discretization of the $t$-derivative, and so it is crucial to make it depend on the sign of $m$ in order to follow the characteristic curves. Steps 3 and 4 simply inverse the reparameterization $\varphi$ obtained after step 2 in order to deduce the horizontal part of $\alpha$. It is important to note that for $s=1$, interpolation between the points $\{x_{k}(s),k=0,\ldots,n\}$ should be achieved so as to remain on the initial shape $\overline{\alpha(1)}$, obtained e.g. by spline interpolation of the points of $\alpha(1)$. Finally, we can perform optimal matching in the discrete case.

Algorithm 5** (Discrete optimal matching)**

Input: $\alpha_{0}=(x_{0}^{0},\ldots,x_{n}^{0}),\alpha_{1}=(x_{0}^{1},\ldots,x_{n}^{1})$.

Initialization : $\tilde{\alpha}_{1}=\alpha_{1}$.

Repeat until convergence :

compute geodesic $s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))$ between $\alpha_{0}$ and $\tilde{\alpha}_{1}$ using Algorithm 3, 2. 2.

compute $s\mapsto\alpha^{hor}(s)=(x_{0}^{hor}(s),\ldots,x_{n}^{hor}(s))$ the horizontal part of $\alpha$ using Algorithm 4, 3. 3.

set $\tilde{\alpha}_{1}=\alpha^{hor}(1)$.

Output : $\alpha_{1}^{opt}=\tilde{\alpha}_{1}$.

4 Simulations

We test the optimal matching (OM) algorithm in several settings : for curves in the hyperbolic half-plane $\mathbb{H}^{2}$, for curves in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$, and for curves on the sphere $\mathbb{S}^{2}$. Regarding the geometry of $\mathbb{H}^{2}$ and the useful algorithms such as the exponential map and the logarithm map, we refer the reader to moi17 . Concerning the geometry of $\mathbb{S}^{2}$, we have used the same formulas as those given in appendix in zhang16 .

We start by comparing Algorithm 5 to the popular dynamic programming (DP) method, presented e.g. in srivastava16 . This alternative approach to optimal matching between two curves $c_{1}$ and $c_{2}$ considers a discrete grid of $[0,1]\times[0,1]$ and tries to find the optimal non-decreasing path $\varphi(t)$ in that grid that puts into correspondence $c_{1}(t)$ and $c_{2}(\varphi(t))$. To do so, at each point $(t_{1},t_{2})$ of the grid, it tests all the possible matchings between the pieces of curves $c_{1}([0,t_{1}])$ and $c_{2}([0,t_{2}])$, and keeps only the one that gives the shortest geodesic. This boils down to testing all the non-decreasing paths that lead to this point. Since this is computationally costly, the algorithm only tests the paths going through the points located at the bottom-left in a square of a certain size, as shown in Figure 6. Even though the computations are additive, this method requires to compute a large number of geodesics. Results of comparison of Algorithm 5 to this method are shown in Figure 8 for a pair of curves in $\mathbb{R}^{2}$ and three pairs of curves in $\mathbb{R}^{3}$. The DP method is carried out for a square of side $s=7$. The first row displays the geodesics between the initial parameterized curves (which are all arc-length parameterized except for the second example) and the second and third rows the horizontal geodesics obtained with the optimal matching and dynamic programing algorithms respectively. We can see that the two methods give very similar results : in the first example, at each extremity, a portion of the circle is matched to almost a single point of the segment, in order to maximize the portions in correspondence where the speeds are positively colinear. In the second example, we consider two curves in $\mathbb{R}^{3}$ that are identical modulo translation and parameterization. This difference in the parameterization induces an artificial deformation of the geodesic (first row), which is "straightened out" in the horizontal geodesics given by both methods. The last two examples illustrate the fact that arc-length parameterization does not always yield a relevant matching : in the third example, it seems more natural to put into correspondence the portions of the curves that are "before the turn", and to match together the ones that are "after the turn", as given by both the OM algorithm and DP. The fourth row shows the corresponding optimal matchings, i.e. the optimal reparameterizations $\varphi$ such that $c_{1}(t)$ is matched to $c_{2}(\varphi(t))$, in red (OM) and black (DP), and we can check in the fifth row the horizontality of the solutions by looking at the ratio, in norm, of the vertical part of the speed vector of the geodesic divided by its horizontal part. We can see that this ratio is largely reduced from the initial geodesic (dashed black line) to the horizontal ones (full red and black lines). Finally, the lengths of the geodesics of the first three lines of Figure 8 are given in Table 1 in the corresponding order, showing that the horizontal geodesics are indeed shorter. The lengths on the first row give the distance between the initial parameterized curves, while the lengths on the second and third row yield a distance between the shapes.

To summarize, it seems that both methods give very similar results when tested on the same metric : they both tend to put into correspondence the parts of the curves that have same shape and orientation. However the dynamic programing approach requires the computation of a large number of geodesics between pieces of curves, whereas in the examples shown here (Figures 9 to 11), the number of iterations required range from 4 to 12, resulting in the same amount of geodesic computations. It should also be noted that it is usually the first iteration that gives the biggest "jump" on the fiber, i.e. that results in the most important reparameterization of the target curve. It could therefore also be used to get a good initialization for some approximate faster method.

We then consider examples where the base manifold has negative curvature. Figure 9 shows horizontal geodesics obtained using the OM algorithm as well as the corresponding optimal matchings in $[0,1]\times[0,1]$, for plane curves (first row) and the same curves taking their values in the hyperbolic half-plane (second row). We can see that the geometry of the base manifold significantly influences the optimal matching between two given curves. To evaluate the performance of the geodesic shooting algorithm used to perform OM, we display in Figure 7 the evolution of the norm of the speed of the geodesic obtained for the two curves on the bottom-right corner of Figure 9 (the vertical and horizontal segments of $\mathbb{H}^{2}$) as the number of points used to compute the geodesic increases (from 20 to 500). The evolutions of the maximum and minimum values $\max_{s\in[0,1]}|c_{s}(s)|$, $\min_{s\in[0,1]}|c_{s}(s)|$ are shown in dashed lines, and that of the mean value with a full line. We can see that the more refined the discretization, the closer we get to a geodesic. It is to be noted that these two segments are the same as those considered in su17 to illustrate an extension of the SRV framework to curves in Lie groups. Since this metric structure coincides with the one studied here in the planar case, it is not surprising to find similar results for the plane curves. However, the results in $\mathbb{H}^{2}$ are quite different.

Finally, we show illustrations on the sphere. Figure 10 shows geodesics before (top row) and after (bottom row) optimal matching, for two pairs of curves representing hurricane tracks from the NASA Tropical Storm Tracks Database111https://ghrc.nsstc.nasa.gov/storms/.. Once again we can see that the optimal matching algorithm seems to "flatten out" the deformation between two curves. In the last figure, we show a set of plane trajectories between Paris and Caracas downloaded from the IAGOS-MOZAIC database 222http://iagos.sedoo.fr/.. It is visually clear that there are two clusters, probably due to different weather conditions. These clusters are easily retrieved using agglomerative hierarchical clustering based on the horizontal geodesic distance. We can find the centers of these clusters by computing the Fréchet means, i.e. for each cluster, the curve that minimizes the sum of the squared distances to all the curves of the cluster. This can be achieved using a Karcher flow algorithm, summarized as follows.

Algorithm 6** (Karcher flow)**

Input: $\alpha_{1},\ldots,\alpha_{N}$.

Initialization : $\hat{\alpha}=\alpha_{1}$.

Repeat until convergence :

for $i=1,\ldots,N$,

• •

compute the horizontal geodesic $s\mapsto\alpha^{hor}(s)$ from $\hat{\alpha}$ to $\alpha_{i}$ using Algorithm 5,

• •

set $w_{i}=(\alpha^{hor})^{\prime}(0)$,

• •

update $\alpha_{i}=\alpha^{hor}(1)$. 2. 2.

set $w=\frac{1}{N}\sum_{i=1}^{N}w_{i}$, 3. 3.

update $\hat{\alpha}=\exp_{\hat{\alpha}}w$.

Output : $\hat{\alpha}$.

The obtained mean curves are shown on the right-hand side of Figure 11.

5 Proof of Theorem 3.1

We conclude this paper with the proof of Theorem 3.1. Let us first remind the result.

Theorem 1

Let $s\mapsto c(s)$ be a $C^{1}$-path of $C^{2}$-curves with non vanishing derivative with respect to $t$. This path can be identified with an element $(s,t)\mapsto c(s,t)$ of $C^{1,2}([0,1]\times[0,1],M)$ such that $c_{t}\neq 0$. Consider the $C^{1}$-path in $M^{n+1}$, $s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))$, that is the discretization of size $n$ of $c$. Then there exists a constant $\lambda>0$ such that for $n$ large enough, the difference between the energies of $c$ and $\alpha$ is bounded by

[TABLE]

where $E$ and $E^{n}$ are the energies with respect to metrics $G$ and $G^{n}$ respectively and where

[TABLE]

and $|w|_{\infty}:=\sup_{s,t\in[0,1]}|w(s,t)|$ denotes the supremum over both $s$ and $t$ of a vector field $w$ along $c$.

Proof (Proof of Theorem 1)

To prove this result, we introduce the unique path $\hat{c}$ of piecewise geodesic curves of which $\alpha$ is the $n$-discretization. It is obtained by linking the points $x_{0}(s),x_{1}(s),\ldots,x_{n}(s)$ of $\alpha$ by pieces of geodesics for all times $s\in[0,1]$

[TABLE]

for $k=0,\ldots,n$. Then the difference between the energy of the path of curves $E(c)$ and the discrete energy of the path of discrete curves $E^{n}(\alpha)$ can be controlled in two steps :

[TABLE]

Step 1. We first consider the difference between the continuous energies of the smooth and piecewise geodesic curves

[TABLE]

Note that the parallel transports $\nabla_{s}q(s,t)^{t,\frac{k}{n}}$ and $\nabla_{s}\hat{q}(s,t)^{t,\frac{k}{n}}$ are performed along different curves – $c(s,\cdot)$ and $\hat{c}(s,\cdot)$ respectively. Let us fix $s\in[0,1]$, $0\leq k\leq n$ and $t\in\left[\frac{k}{n},\frac{k+1}{n}\right]$. From now on we will omit "$s$" in the notation $w(s,t)$ to lighten notations. Using the notation $w^{\parallel}(t):=w(t)^{t,\frac{k}{n}}$ to denote the parallel transport of a vector field $w$ from $t$ to $\frac{k}{n}$ along its baseline curve, the difference we need to control is

[TABLE]

Since $|w-\frac{1}{2}w^{T}|\leq|w|$ for any vector $w$, we can write

[TABLE]

Let us first consider the difference $|c_{t}^{\parallel}-\hat{c}_{t}^{\parallel}|$. Since $\hat{c}_{t}(t)^{t,\frac{k}{n}}=\hat{c}_{t}(\frac{k}{n})$, we can write

[TABLE]

The first term is smaller than $1/n\cdot|\nabla_{t}c_{t}|_{\infty}$. To bound the second term, we place ourselves in a local chart $(\varphi,U)$ centered in $c(\tfrac{k}{n})=c(s,\tfrac{k}{n})$, such that $c([0,1]\times[0,1])\subset U$. After identification with an open set of $\mathbb{R}^{d}$ – where $d$ is the dimension of the manifold $M$– using this chart, we get

[TABLE]

Since a geodesic locally looks like a straight line (see e.g. emery84 ) there exists a constant $\lambda_{1}$ such that

[TABLE]

and so

[TABLE]

The second derivative in $t$ of the coordinates of $c$ in the chart $(U,\varphi)$ can be written ${c_{tt}}^{\ell}={\nabla_{t}c_{t}}^{\ell}-\Gamma_{ij}^{\ell}{c_{t}}^{i}{c_{t}}^{j}$ for $\ell=1,\ldots,d$, and so there exists a constant $\lambda_{2}$ such that $|c_{tt}|\leq\lambda_{2}\left(|\nabla_{t}c_{t}|_{\infty}+|c_{t}|_{\infty}^{2}\right)$, and

[TABLE]

This means that for $n$ large enough, we can write e.g.

[TABLE]

From (24) we can also deduce that

[TABLE]

Let us now consider the difference $|{\nabla_{s}c_{t}}^{\parallel}-{\nabla_{s}\hat{c}_{t}}^{\parallel}|$. Since $c_{s}(s,\tfrac{k}{n})=\hat{c}_{s}(s,\tfrac{k}{n})$, we get

[TABLE]

and so

[TABLE]

We can decompose $\nabla_{t}^{2}\hat{c}_{s}(s,t)=\nabla_{t}\nabla_{s}\hat{c}_{t}(s,t)=\nabla_{s}\nabla_{t}\hat{c}_{t}(s,t)+\mathcal{R}(\hat{c}_{t},\hat{c}_{s})\hat{c}_{t}(s,t)$, and since $\nabla_{t}\hat{c}_{t}(s,t)=0$ and $|\mathcal{R}(X,Y)Z|\leq|K|\cdot(|\langle Y,Z\rangle||X|+|\langle X,Z,\rangle||Y|)\leq 2|K|\cdot|X|\cdot|Y|\cdot|Z|$ by Cauchy Schwarz, we get using Equation (25)

[TABLE]

To bound $|\hat{c}_{s}|$ we apply Lemma 1 to the Jacobi field $J:[0,1]\ni u\mapsto\hat{c}_{s}(s,\frac{k+u}{n})$ along the geodesic $\gamma(u)=\hat{c}(s,\frac{k+u}{n})$, that is

[TABLE]

where, since $\gamma^{\prime}(0)=\frac{1}{n}\hat{c}_{t}(s,\frac{k}{n})=\tau_{k}(s)$, the coefficients are defined by

[TABLE]

This gives $J(1)^{1,0}=J(0)^{T}+a_{k}(1)J(0)^{N}+\nabla_{t}J(0)^{N}+b_{k}(1)\nabla_{t}J(0)^{N}$ and so

[TABLE]

Injecting this into (29), we obtain since $u=nt-k$ and $\hat{c}_{s}(s,\tfrac{k}{n})=c_{s}(s,\frac{k}{n})$,

[TABLE]

When $n\to\infty$, $a_{k}(1)\to 1$, $b_{k}(1)\to 1$, and since $0\leq nt-k\leq 1$, $a_{k}(nt-k)\to 1$, $b_{k}(nt-k)\to 1$ also. Therefore, for $n$ large enough we can see that $|\hat{c}_{s}|\leq\lambda_{4}|c_{s}|_{\infty}$ for some constant $\lambda_{4}$. Injecting this into (28) gives

[TABLE]

and so the difference (27) can be bounded by

[TABLE]

for some constant $\lambda_{5}$. Injecting (25), (26) and (31) in Equation (23) we obtain

[TABLE]

for some constant $\lambda_{6}$. To conclude this first step, let us bound the sum

[TABLE]

Taking the derivative according to $t$ on both sides of (30), we get since $n|\tau_{k}(s)|=|\hat{c}_{t}(s,\frac{k}{n})|$,

[TABLE]

since derivation of the coefficients give $b_{k}^{\prime}(u)=a_{k}(u)$ and $a_{k}^{\prime}(u)=|\tau_{k}|e_{k}(u)=\frac{1}{n}|\hat{c}_{t}(\frac{k}{n})|e_{k}(u)$, where

[TABLE]

Since the coefficients $e_{k}(nt-k)$, $a_{k}(nt-k)/b_{k}(1)$ and $a_{k}(1)$ are bounded for $n$ large enough, and since $|\hat{c}_{t}|\leq\frac{3}{2}|c_{t}|_{\infty}$, we can write for some constant $\lambda_{7}$,

[TABLE]

Inserting this into (33) gives

[TABLE]

Finally, we are able to bound the difference between the energies of the smooth and piecewise-geodesic paths by combining Equations (5) and (35)

[TABLE]

Step 2. Let us now consider the difference of energy between the path of piecewise geodesic curves and the path of discrete curves. Since $\nabla_{s}q_{k}(s)=\nabla_{s}\hat{q}(s,\tfrac{k}{n})$ for all $s\in[0,1]$ and $0\leq k\leq n$, we can write

[TABLE]

We fix once again $s\in[0,1]$, $0\leq k\leq n$ and $t\in\left[\frac{k}{n},\frac{k+1}{n}\right]$. As in step 1, we will omit "$s$" in most notations. Since $|\hat{c}_{t}(t)|=|\hat{c}_{t}(\frac{k}{n})|$, we get

[TABLE]

Considering once again the Jacobi field

[TABLE]

along the geodesic $\gamma(u)=\hat{c}(\frac{k+u}{n})$, Equation (29) gives

[TABLE]

Recall that $b_{k}^{\prime}(u)=a_{k}(u)$ and $a_{k}^{\prime}(u)=|\tau_{k}|e_{k}(u)$, and so taking the derivative with respect to $t$ and decomposing $\nabla_{t}\hat{c}_{s}(\tfrac{k}{n})^{T}=\nabla_{t}\hat{c}_{s}(\tfrac{k}{n})-\nabla_{t}\hat{c}_{s}(\tfrac{k}{n})^{N}$, we obtain

[TABLE]

Noticing that $\frac{e_{k}(nt-k)}{(nt-k)|\tau_{k}|}\to 1$ and $\frac{a_{k}(nt-k)-1}{(nt-k)|\tau_{k}|}\to 0$ when $n\to\infty$, we can deduce that for $n$ large enough,

[TABLE]

This gives

[TABLE]

Recall from (34) and (35) that

[TABLE]

and so

[TABLE]

Finally, we obtain

[TABLE]

which completes the proof.

Acknowledgments

This research was supported by Thales Air Systems and the french MoD DGA. MOZAIC/CARIBIC/IAGOS data were created with support from the European Commission, national agencies in Germany (BMBF), France (MESR), and the UK (NERC), and the IAGOS member institutions (http://www.iagos.org/partners). The participating airlines (Lufthansa, Air France, Austrian, China Airlines, Iberia, Cathay Pacific, Air Namibia, Sabena) supported IAGOS by carrying the measurement equipment free of charge since 1994. The data are available at http://www.iagos.fr thanks to additional support from AERIS.

Appendix A

Lemma 1 * Let $\gamma:[0,1]\rightarrow M$ be a geodesic of a manifold $M$ of constant sectional curvature $K$, and $J$ a Jacobi field along $\gamma$. Then the parallel transport of $J(t)$ along $\gamma$ from $\gamma(t)$ to $\gamma(0)$ is given by*

[TABLE]

for all $t\in[0,1]$, where

[TABLE]

Proof

As a Jacobi field along $\gamma$, $J$ satisfies the well-known equation (see e.g. jost )

[TABLE]

If $M$ is flat, we get $\nabla_{t}^{2}J(t)=0$ and so $J(t)=J(0)+t\nabla_{t}J(0)$. If not, we can decompose $J$ in the sum $J=J^{T}+J^{N}$ of two vector fields that parallel translate along $\gamma$, $J^{T}=\langle J,v\rangle v$ with $v=\gamma^{\prime}/|\gamma^{\prime}|$, and $J^{N}=J-J^{T}$. Since $\langle\nabla_{t}^{2}J,\gamma^{\prime}\rangle=0$ and $\gamma^{\prime}$ is parallel along $\gamma$, we get by integrating twice that

[TABLE]

Since

[TABLE]

the normal component $J^{N}$ is also a Jacobi field, that is it verifies

[TABLE]

And since $M$ has constant sectional curvature $K$, for any vector field $w$ along $\gamma$ we have

[TABLE]

and the differential equation verified by $J^{N}$ can be rewritten $\nabla_{t}^{2}J^{N}(t)=-K\,|\gamma^{\prime}(t)|^{2}\,J^{N}(t)$. Since the speed of the geodesic $\gamma$ has constant norm, the solution to that differential equation is of the form

[TABLE]

when $K=-1$ and

[TABLE]

when $K=1$. Using the initial conditions $J^{N}(0)$ and $\nabla_{t}J^{N}(0)$ to find the constants $\lambda,\mu$, we get for $K=-1$

[TABLE]

and for $K=1$, the same formula with cosine and sine functions instead of hyperbolic cosine and sine.

Lemma 3* The covariant derivatives of the functions $f_{k}^{(-)}$ and $g_{k}^{(-)}$ with respect to $s$ are functions $T_{x_{k+1}(s)}M\rightarrow T_{x_{k}(s)}M$ given by*

[TABLE]

where

[TABLE]

*if $K$ is the sectional curvature of the base manifold. *

Proof

Fix $0\leq k\leq n$ and let $w_{k+1}:s\mapsto w_{k+1}(s)$ be a vector field along the curve $x_{k+1}:s\mapsto x_{k+1}(s)$. By definition,

[TABLE]

Consider the path of gedesics $s\mapsto\gamma_{k}(s,\cdot)$ such that for all $s\in[0,1]$, $\gamma_{k}(s,0)=x_{k}(s)$, $\gamma_{k}(s,1)=x_{k+1}(s)$ and $t\mapsto\gamma_{k}(s,t)$ is a geodesic. We denote by ${w_{k+1}}^{\parallel}$ the vector field along the curve $x_{k}$ obtained by parallel transporting back the vector $w_{k+1}(s)$ along the geodesic $\gamma_{k}(s,\cdot)$ for all $s\in[0,1]$, i.e. ${w_{k+1}}^{\parallel}(s)=P_{\gamma_{k}(s,\cdot)}^{1,0}(w_{k+1}(s))$. We have

[TABLE]

and so we need to compute $\nabla_{s}({w_{k+1}}^{\parallel})$. Let $V(s,t):=P_{\gamma_{k}(s,\cdot)}^{1,t}(w_{k+1})$ so that $\nabla_{s}V(s,1)=\nabla_{s}w_{k+1}$ and $\nabla_{s}V(s,0)=\nabla_{s}({w_{k+1}}^{\parallel})$, then

[TABLE]

since $\nabla_{t}V=0$, and where $\partial_{t}\gamma_{k}(s,t)^{t,0}=\tau_{k}(s)$. We get, since $\nabla R=0$,

[TABLE]

To find an expression for ${\partial_{s}\gamma_{k}}^{t,0}$, we consider the Jacobi field $J(t):=\partial_{s}\gamma_{k}(s,t)$ along the geodesic $t\mapsto\gamma_{k}(s,t)$. The vector field $J$ verifies

[TABLE]

where the last equality results from the inversion $\nabla_{t}\partial_{s}\gamma_{k}(s,0)=\nabla_{s}\partial_{t}\gamma_{k}(s,0)$ and $\partial_{t}\gamma_{k}(s,0)=\tau_{k}(s)$. Applying Lemma 1 gives, for all $k=0,\ldots,n-1$,

[TABLE]

with the coefficients

[TABLE]

[TABLE]

Integrating this and injecting it in (37) gives

[TABLE]

where $Y_{k}$ is defined by

[TABLE]

and injecting this in (36) finally gives,

[TABLE]

which is what we wanted. The covariant derivative \nabla_{s}\big{(}g_{k}^{(-)}(w_{k+1})\big{)} can be computed in a similar way.

Appendix B

Proposition 6 (Discrete geodesic equations)* A path $s\mapsto\alpha(s)=\left(x_{0}(s),\ldots,x_{n}(s)\right)$ in $M^{n+1}$ is a geodesic for metric $G^{n}$ if and only if its SRV representation s\mapsto\big{(}x_{0}(s),(q_{k}(s))_{k}\big{)} verifies the following differential equations*

[TABLE]

*for all $k=0,\ldots,n-1$, with the notations (13) and $R_{k}:=\mathcal{R}(q_{k},\nabla_{s}q_{k}){x_{k}}^{\prime}$. *

Proof

We consider a variation $(-\delta,\delta)\ni a\mapsto\alpha(a,\cdot)=(x_{0}(a,\cdot),\ldots,x_{n}(a,\cdot))$ of this curve which coincides with $\alpha$ for $a=0$, i.e. $\alpha(0,s)=\alpha(s)$ for all $s\in[0,1]$, and which preserves the end points of $\alpha$, i.e. $\alpha(a,0)=\alpha(0)$ and $\alpha(a,1)=\alpha(1)$ for all $a\in(-\delta,\delta)$. The energy of this variation with respect to metric $G^{n}$ can be seen as a real function of the variable $a$ and is given by

[TABLE]

and its derivative $(E^{n})^{\prime}(a)$ with respect to $a$ is

[TABLE]

where we integrate by parts to obtain the third line from the second. The goal is to express $\partial_{a}x_{k}$ in terms of $\partial_{a}x_{0}$ and $\nabla_{a}q_{\ell}$, $\ell=0,\cdots,k$. That way, the only elements that depend on $a$ once we take $a=0$ are $(\partial_{a}x_{0},\nabla_{a}q_{0},\cdots,\nabla_{a}q_{n-1})$ which can be chosen independently to be whatever we want. Let us fix $0\leq k\leq n-1$ and $s\in[0,1]$ and consider the path of geodesics $a\mapsto\gamma_{k}(a,\cdot)$ such that $\gamma_{k}(a,0)=x_{k}(a,s)$, $\gamma_{k}(a,1)=x_{k+1}(a,s)$ and $\partial_{t}\gamma_{k}(a,0)=\tau_{k}(a,s)=\log_{x_{k}(a,s)}(x_{k+1}(a,s))$. Then by definition, for each $a\in[0,1]$, $t\mapsto J(a,t):=\partial_{a}\gamma_{k}(a,t)$ is a Jacobi field along the geodesic $t\mapsto\gamma_{k}(a,t)$ of $M$, and so Lemma 1 gives

[TABLE]

where ${\partial_{a}x_{k+1}}^{\parallel}$ denotes the parallel transport of $\partial_{a}x_{k+1}$ from $x_{k+1}(s)$ to $x_{k}(s)$ along the geodesic. Differentiation of $q_{k}=\sqrt{n}\,\tau_{k}/|\tau_{k}|$ gives

[TABLE]

and taking the tangential part on both sides yields ${\nabla_{s}q_{k}}^{T}=\sqrt{n}\,|\tau_{k}|^{-1/2}\frac{1}{2}{\nabla_{s}\tau_{k}}^{T}$, and so finally

[TABLE]

Injecting this in (39) and noticing that $\langle f_{k}(w),z\rangle=\langle w,f_{k}(z)\rangle$ and $\langle g_{k}(w),z\rangle=\langle w,g_{k}(z)\rangle$ for any pair of vectors $w,z$ gives

[TABLE]

for any tangent vector $w_{k+1}\in T_{x_{k+1}}M$. From equation (41) we can deduce, for $k=1,\ldots,n$,

[TABLE]

With the notation $R_{k}:=\mathcal{R}(q_{k},\nabla_{s}q_{k}){x_{k}}^{\prime}$ we get

[TABLE]

and we can then write the derivative of the energy $(E^{n})^{\prime}(0)$ for $a=0$ in the following way

[TABLE]

where in the first sum we use the notation convention $f_{0}\circ\cdots\circ f_{-1}:=\text{Id}$. Noticing that the double sum can be rewritten

[TABLE]

we obtain for $(E^{n})^{\prime}(0)$

[TABLE]

where in the last sum we use the convention $\sum_{\ell=n}^{n-1}=0$. Since this quantity has to vanish for any choice of $(\partial_{a}x_{0}(0,\cdot),\nabla_{a}q_{0}(0,\cdot),\ldots,\nabla_{a}q_{n-1}(0,\cdot))$, the geodesic equations for the discrete metric are

[TABLE]

for all $k=0,\ldots,n-1$, with the conventions $\sum_{\ell=n}^{n-1}=0$ and $f_{0}\circ\cdots\circ f_{-1}:=\text{Id}$.

Remark 4* Let $[0,1]\ni s\mapsto c(s,\cdot)\in\mathcal{M}$ be a $C^{1}$ path of smooth curves and $[0,1]\ni s\mapsto\alpha(s)\in M^{n+1}$ the discretization of size $n$ of $c$. We denote as usual by $q:=c_{t}/|c_{t}|^{1/2}$ and $(q_{k})_{k}$ their respective SRV representations. When $n\rightarrow\infty$ and $|\tau_{k}|\rightarrow 0$ while $n|\tau_{k}|$ stays bounded for all $0\leq k\leq n$, the coefficients of the discrete geodesic equation (20) for $\alpha$ converge to the coefficients of the continuous geodesic equation (7) for $c$, i.e.*

[TABLE]

for all $s\in[0,1]$ and $k=0,\ldots,n-1$, where $r_{n-1}=0$ and for $k=1,\ldots,n-2$,

[TABLE]

*with the exception that the sum starts at $\ell=0$ for $r_{0}$. *

Proof

This is due to three arguments : (1) at the limit, $f_{k}(w)=w+o(1/n)$ and $g_{k}(w)=|q_{k}|(w+{w}^{T})+o(1/n)$, (2) parallel transport along a piecewise geodesic curve uniformly converges to the parallel transport along the limit curve, and (3) the discrete curvature term $R_{k}(s)$ converges to the continuous curvature term $\mathcal{R}(q,\nabla_{s}q)c_{s}(s,\tfrac{k}{n})$ for all $k$. Indeed, let $\hat{c}$ be the unique piecewise geodesic curve of which $\alpha$ is the discretization, i.e. c\big{(}\frac{k}{n}\big{)}=\hat{c}\big{(}\frac{k}{n}\big{)}=x_{k} for all $k=0,\ldots,n$ and $\hat{c}$ is a geodesic on each segment \big{[}\tfrac{k}{n},\tfrac{k+1}{n}\big{]}. Defining

[TABLE]

the geodesic equations can be written in terms of the vectors $\hat{r}_{k}$

[TABLE]

We can show that for any $0\leq k\leq\ell\leq n-2$ and any vector $w\in T_{x_{\ell+1}}M$,

[TABLE]

Since $|a_{j}-1|/|\tau_{k}|^{2}\to 0$ when $n\to\infty$ and $n|\tau_{k}|$ stays bounded, we have for all $0\leq j\leq n$ and $n$ large enough $|a_{j}-1|\leq\frac{1}{n^{2}}$, and using the fact that parallel transport along a piecewise geodesic curve uniformly converges to the parallel transport along the limit curve, we get

[TABLE]

when $n\to\infty$. Now, denoting by

[TABLE]

the curvature term involved in the continuous geodesic equations, we have since ${x_{k}}^{\prime}(s)=c_{s}(s,\tfrac{k}{n})$ and $|\mathcal{R}(X,Y)Z|\leq|K|\cdot(|\langle Y,Z\rangle||X|+|\langle X,Z,\rangle||Y|)\leq 2|K|\cdot|X|\cdot|Y|\cdot|Z|$ by Cauchy Schwarz,

[TABLE]

Let us show that both summands of this upper bound tend to [math] when $n\to\infty$.

[TABLE]

and since the portion of $c(s,\cdot)$ on the segment $[\tfrac{k}{n},\tfrac{k+1}{n}]$ is close to a geodesic at the limit, $|n\tau_{k}-c_{t}(\tfrac{k}{n})|\to 0$ when $n\to\infty$, and so does $|q_{k}(s)-q(\tfrac{k}{n})|$. Similarly,

[TABLE]

where once again $||n\tau_{k}|^{-1/2}-|c_{t}(\tfrac{k}{n})|^{-1/2}|\to 0$ and $|n\nabla_{s}\tau_{k}|$ is bounded. The last term can be bounded, for $n$ large enough, by

[TABLE]

since $\nabla_{s}\tau_{k}=(D_{\tau}\alpha^{\prime})_{k}=({x_{k+1}}^{\parallel}-x_{k})^{T}+b_{k}^{-1}({x_{k+1}}^{\prime}-{x_{k}}^{\prime})^{N}$ and $b_{k}^{-1}\to 1$. Finally, we can see that

[TABLE]

goes to [math] when $n\to\infty$. We can show in a similar way that $|g_{k}^{(-)}(\hat{r}_{k})-|q_{k}|(r_{k}+{r_{k}}^{T})|\to 0$ when $n\to\infty$.

Proposition 7 (Discrete exponential map)* Let $[0,1]\ni s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))$ be a geodesic path in $M^{n+1}$. For all $s\in[0,1]$, the coordinates of its acceleration $\nabla_{s}\alpha^{\prime}(s)$ can be iteratively computed in the following way*

[TABLE]

for $k=0,\ldots,n-1$, where the $R_{k}$’s are defined as in Proposition 7, the symbol $\cdot^{\parallel}$ denotes the parallel transport from $x_{k+1}(s)$ back to $x_{k}(s)$ along the geodesic linking them, the maps $\nabla_{s}f_{k}$ and $\nabla_{s}g_{k}$ are given by Lemma 2, $Y_{k}$ is given by Equation (Appendix A) and

[TABLE]

Proof

For all $s\in[0,1]$, we initialize $\nabla_{s}{x_{k}}^{\prime}(s)$ for $k=0$ using the first geodesic equation in (20); the difficulty lies in deducing $\nabla_{s}{x_{k+1}}^{\prime}(s)$ from $\nabla_{s}{x_{k}}^{\prime}(s)$. Just as we have previously obtained (40), we can obtain by replacing the derivatives with respect to $a$ by derivatives with respect to $s$

[TABLE]

and by differentiating with respect to $s$

[TABLE]

We have already computed (38) the covariant derivative of a vector field $s\mapsto{w_{k+1}(s)}^{\parallel}\in T_{x_{k}(s)}M$ and so we can write

[TABLE]

where $Y_{k}$ is defined by Equation (Appendix A). Together with Equation (44), this gives the desired equation for ${\nabla_{s}{x_{k+1}}^{\prime}}^{\parallel}$. Finally, $\nabla_{s}\tau_{k}=(D_{\tau}\alpha^{\prime})_{k}$ results directly from (43), $\nabla_{s}^{2}q_{k}$ is deduced from the second geodesic equation and the remaining equations follow from simple computation.

Proposition 8 (Discrete Jacobi fields)* Let $[0,1]\ni s\mapsto\alpha(s)=(x_{0}(s),\ldots,x_{n}(s))$ be a geodesic path in $M^{n+1}$, $[0,1]\ni s\mapsto J(s)=(J_{0}(s),\ldots,J_{n}(s))$ a Jacobi field along $\alpha$, and $(-\delta,\delta)\ni a\mapsto\alpha(a,\cdot)$ a corresponding family of geodesics, in the sense just described. Then $J$ verifies the second order linear ODE*

[TABLE]

for all $0\leq k\leq n-1$, where $R_{k}:=\mathcal{R}(q_{k},\nabla_{s}q_{k}){x_{k}}^{\prime}$ and the various covariant derivatives according to $a$ can be expressed as functions of $J$ and $\nabla_{s}J$,

[TABLE]

[TABLE]

with the notation conventions $f_{k+1}^{(-)}\circ\ldots\circ f_{k-1}^{(-)}:=Id$, $\sum_{\ell=n}^{n-1}:=0$ and with the maps

[TABLE]

[TABLE]

and

[TABLE]

Proof

For all $a\in(-\delta,\delta)$, $\alpha(a,\cdot)$ verifies the geodesic equations (20). Taking the covariant derivative of these equations according to $a$ we obtain

[TABLE]

Since for $a=0$, $\nabla_{a}\nabla_{s}\partial_{s}x_{0}=\nabla_{s}^{2}J_{0}+\mathcal{R}(J_{0},\partial_{s}x_{0})\partial_{s}x_{0}$, we get

[TABLE]

and the differentiation

[TABLE]

gives the desired equation for $\nabla_{s}^{2}J_{0}$. Now we will try to deduce $\nabla_{s}^{2}J_{k+1}$ from (46). If ${J_{k+1}}^{\parallel}(s)$ denotes the parallel transport of the vector $J_{k+1}(s)$ from $x_{k+1}(s)$ back to $x_{k}(s)$ along the geodesic that links them, we know from (40) that

[TABLE]

We also know from (38) that

[TABLE]

and by iterating

[TABLE]

Developping and injecting Equation (47) in the latter gives

[TABLE]

Developping the covariant derivatives \nabla_{s}^{2}\big{(}f_{k}(J_{k})\big{)} and \nabla_{s}^{2}\big{(}g_{k}(\nabla_{a}q_{k})\big{)} gives the desired formula. Now let us explicit the different terms involved in these differential equations. Since $\nabla\mathcal{R}=0$ and $\nabla_{a}\partial_{s}x_{k}=\nabla_{s}\partial_{a}x_{k}$, we have

[TABLE]

By taking the inverse of (47) we get

[TABLE]

and taking the derivative according to $s$ on both sides and injecting Equation (48) gives

[TABLE]

To obtain $\nabla_{s}^{2}\nabla_{a}q_{k}$, notice that

[TABLE]

and injecting Equation (46) with

[TABLE]

gives us the desired formula. $\nabla_{s}Y_{k}$ results from simple differentiation, and differentiating the maps $f_{k}^{(-)}$ and $g_{k}^{(-)}$ with respect to $a$ is completely analogous to the the computations of Lemma 3. Finally, the inverse of $g_{k}$ is given by ${g_{k}}^{-1}:T_{x_{k}}M\rightarrow T_{x_{k}}M$,

[TABLE]

and since

[TABLE]

it is straightforward to verify that

[TABLE]

gives

[TABLE]