# Numerical Accuracy of Ladder Schemes for Parallel Transport on Manifolds

Nicolas Guigui, Xavier Pennec

► To cite this version:

Nicolas Guigui, Xavier Pennec. Numerical Accuracy of Ladder Schemes for Parallel Transport on Manifolds. Foundations of Computational Mathematics, 2022, 22, pp.757-790. 10.1007/s10208-021-09515-x. hal-02894783v3

HAL Id: hal-02894783

https://inria.hal.science/hal-02894783v3

Submitted on 21 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical Accuracy of Ladder Schemes for Parallel Transport on Manifolds

Nicolas Guigui<fnref n="1" />,² · Xavier Pennec<fnref n="1" />,³

Received: 9 July 2020 / Revised: 8 March 2021 / Accepted: 22 March 2021

© The Author(s) 2021

## Abstract

Parallel transport is a fundamental tool to perform statistics on Riemannian manifolds. Since closed formulae do not exist in general, practitioners often have to resort to numerical schemes. *Ladder* methods are a popular class of algorithms that rely on iterative constructions of geodesic parallelograms. And yet, the literature lacks a clear analysis of their convergence performance. In this work, we give Taylor approximations of the elementary constructions of Schild’s ladder and the pole ladder with respect to the Riemann curvature of the underlying space. We then prove that these methods can be iterated to converge with quadratic speed, even when geodesics are approximated by numerical schemes. We also contribute a new link between Schild’s ladder and the Fanning scheme which explains why the latter naturally converges only linearly. The extra computational cost of ladder methods is thus easily compensated by a drastic reduction of the number of steps needed to achieve the requested accuracy. Illustrations on the 2-sphere, the space of symmetric positive definite matrices and the special Euclidean group show that the theoretical errors we have established are measured with a high accuracy in practice. The special Euclidean group with an anisotropic left-invariant metric is of particular interest as it is a tractable example of a non-symmetric space in general, which reduces to a Riemannian symmetric space in a particular case. As a secondary contribution, we compute the covariant derivative of the curvature in this space.

**Keywords** Riemannian geometry · Parallel transport · Numerical scheme

**Mathematics Subject Classification** 53a35 · 53B21 · 65D30

Communicated by Hans Munthe-Kaas.

<figure><img src="image_2.png" /></figure>

Nicolas Guigui
nicolas.guigui@inria.fr

Xavier Pennec
xavier.pennec@inria.fr

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹ Université Côte d'Azur and Inria, Epione Team, 2004 route des Lucioles, 06902 Sophia Antipolis, France</footnote>
# 1 Introduction

In many applications, it is natural to model data as points that lie on a manifold. Consequently, there has been a growing interest in defining a consistent framework to perform statistics and machine learning on manifolds [31]. A fruitful approach is to locally linearize the data by associating to each point a tangent vector. The parallel transport of tangent vectors then appears as a natural tool to compare tangent vectors across tangent spaces. For example, [2] uses it on the manifold of symmetric positive definite (SPD) matrices to centralize batches when training a neural network. In [33], parallel transport is used again on SPD matrices for domain adaptation. In [17], it is a key ingredient to spline-fitting in a Kendall shape space. In computational anatomy, it allows to compare longitudinal, intra-subject evolution across populations [3,19,20,32]. It is also used in computer vision in, e.g., [6,15].

However, there is usually no closed-form solution to compute parallel transport and one must use approximation schemes. Two classes of approximations have been developed. The first intends to solve ordinary or partial differential equations (ODE/PDE) related to the definition of parallel transport itself or to the Jacobi fields [22,34]. For instance, [17] derived a homogeneous first-order differential equation stemming from the structure of quotient space that defines Kendall shape spaces. [22] leveraged the relation between Jacobi fields and parallel transport to derive a numerical scheme that amounts to integrating the geodesic equations. They prove that a convergence speed of order one is reached. This scheme is particularly appealing as only the Hamiltonian of the metric is required, and both the main geodesic and the parallel transport are computed simultaneously when no-closed form solution is available for the geodesics.

The second class of approximations, referred to as *ladder methods* [19,20], consists in iterating elementary constructions of geodesic parallelograms (Fig. 1). The most famous, Schild's ladder (SL), was originally introduced by Alfred Schild in 1970 although no published reference exist.<fnref n="1" /> Its first appearance in the literature is in [25] where it is used to introduce the Riemannian curvature. A proof that the construction of a geodesic parallelogram—i.e., one *rung* of the ladder—is an approximation of parallel transport was first given in [16]. However, there is currently no formal proof that it converges when iterated along the geodesic along which the vector is to be transported, as prescribed in every description of the method, and because only the first order is given in [16], Schild's ladder is considered to be a first-order method in the literature.

The first aim of this paper is at filling this gap. Using a Taylor expansion of the *double exponential* introduced by [8,10], we compute an expansion of the SL construction up to order three, with coefficients that depend on the curvature tensor and its covariant derivatives. Understanding this expansion then allows to give a proof of the convergence of the numerical scheme with iterated constructions, and to improve the scheme to a second-order method. We indeed prove that an arbitrary speed of convergence in $\frac{1}{n\alpha}$ can be obtained where $1 \le \alpha \le 2$ and $n$ is the number of parallelograms, or rungs of the ladder. This improves over the assumption made in [22] that only a first-order convergence speed can be reached with SL. These results are observed with

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<footnote n="1">¹ [5] is often cited but no mention of the scheme is made in this work.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 1** Schild's ladder on the sphere. The main geodesic is in blue, while the sides of the intermediate geodesic parallelograms are in dotted green and the diagonals in dotted red. The tangent vectors are represented by (scaled) black arrows. The construction is detailed in Sect. 2.1 (Color figure online)</figcaption></figure>

a high accuracy in numerical experiments performed using the open-source Python package `geomstats`<fnref n="2" /> [24].

A slight modification of this scheme, pole ladder (PL), was introduced in [20]. It turns out that this scheme is exact in symmetric spaces, and in general that an elementary construction of PL is even more precise by an order of magnitude than that of SL [28]. With the same method as for SL, we give a new proof of this result. By applying the previous analysis to PL, we demonstrate that a convergence speed of order two can be reached. Furthermore, we introduce a new construction which consists in averaging two SL constructions and show that it is of the same order as the PL.

In most cases, however, the exponential and logarithm maps are not available in closed form, and one also has to resort to numerical integration schemes. As for the Fanning scheme (FS) [22], we prove that the ladder methods converge when using approximate geodesics and that all geodesics of the construction may be computed in one pass—i.e., using one integration step (e.g., Runge-Kutta) per parallelogram construction, thus reducing the computational cost. We study the FS under the hood of ladder methods and show that its implementation make it very close to SL. However, their elementary steps differ at the third order, so that the FS cannot be improved to a second-order method like SL or the PL.

To observe the impact of the different geometric structures on the convergence, we study the Lie group of isometries of $\mathbb{R}^d$, the special Euclidean group $SE(d)$. Endowed with a left-invariant metric, this space is a Riemannian symmetric space if the metric is isotropic [35]. However, we show that it is no longer symmetric when using anisotropic metrics and that geodesics may be computed by integrating the Euler-Poincaré equations. The same implementation thus allows to observe both geometric structures. Furthermore, the curvature and its derivative may be computed explicitly [23] and confirm our predictions. We treat this example in detail to demonstrate the impact of curvature on the convergence, and the code for the computations, and all

<footnote n="2">² available at http://geomstats.ai</footnote>
the experiments of the paper are available online at https://github.com/nguigs/ladder-methods.

The first part of this paper is dedicated to Schild's ladder, while the second part applies the same methodology to the modified constructions. All the results are illustrated by numerical experiments along the way. The main proofs are in the text, and remaining details can be found in the appendices.

## 1.1 Notations and Assumptions

We consider a complete Riemannian manifold $(\mathcal{M}, g)$ of finite dimension $d \in \mathbb{N}$. The associated Levi-Civita connection defines a covariant derivative $\nabla$ and the parallel transport map. Denote $\exp$ the Riemannian exponential map, and log its inverse, defined locally. For $x \in \mathcal{M}$, let $T_x \mathcal{M}$ be the tangent space at $x$. The map $\exp_x$ sends $T_x \mathcal{M}$ to (a subset of) $\mathcal{M}$, and we will often write $x_w = \exp_x(w) \in \mathcal{M}$ for $w \in T_x \mathcal{M}$.

Let $\gamma : [0, 1] \to \mathcal{M}$ be a smooth curve with $\gamma(0) = x$. For $v \in T_x \mathcal{M}$, the parallel transport $\Pi_{\gamma(0)v}^t v \in T_{\gamma(t)}\mathcal{M}$ of $v$ along $\gamma$ is defined as the unique solution at time $t \le 1$ to the ODE $\nabla_{\dot{\gamma}(s)}X(s) = 0$ with $X(0) = v$. In general, parallel transport depends on the curve followed between two points. In this work, however, we focus on the case where $\gamma$ is a geodesic starting at $x$, and let $w$ be its initial velocity, i.e., $\gamma(t) = \exp_x(tw)$. Thus, the dependence on $\gamma$ in the notation $\Pi_{\gamma}$ will be omitted, and we instead write $\Pi_x^y$ for the parallel transport along the geodesic joining $x$ to $y$ when it exists and is unique. The methods developed below can be extended straightforwardly to piecewise geodesic curves.

We denote by $\|\cdot\|$ the norm defined on each tangent space by the metric $g$ and $R$ the Riemann curvature tensor that maps for any $x \in \mathcal{M}$, $u, v, w \in T_x \mathcal{M}$ to $R(u, v)w \in T_x \mathcal{M}$. Throughout this paper, we consider $\gamma$ to be contained in a compact set $K \subset \mathcal{M}$ of diameter $\delta > 0$, and thus, $R$ and all its covariant derivatives $\nabla^n R$ can be uniformly bounded.

## 1.2 Double Exponential and Neighboring Logarithm

We now introduce the main tool to compute a Taylor approximation of the ladder constructions, the double exponential (also written $\exp$) [8,10]. It is defined for $x \in \mathcal{M}$, $(v, w) \in (T_x \mathcal{M})^2$ by first applying $\exp_x$ to $v$, then $\exp_{x_v}$ to the parallel transport of $w$ along the geodesic from $x$ to $x_v$:

$$ \exp_x(v, w) = \exp_{x_v}(\Pi_x^{x_v} w). $$

As the composition of smooth maps, it is also a smooth map. As the exponential map is locally one-to-one and the parallel transport is an isomorphism, we may define the function $h_x : U_x \to T_x \mathcal{M}$ on an open neighborhood $U_x$ of $(0, 0) \in T_x \mathcal{M} \times T_x \mathcal{M}$ such that $\exp_x(U_x)$ is contained in a convex neighborhood of $\mathcal{M}$ and (Fig. 2):

$$ \exp_x(h_x(v, w)) = \exp_x(v, w) \quad \forall v, w \in T_x \mathcal{M}. $$

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 2 The double exponential and its inverse (left) and the neighboring logarithm (right) represented in a normal coordinate system centered at x. Geodesics are in black, and tangent vectors are colored (Color figure online)</figcaption></figure>

As it is a smooth map, we can apply the Taylor theorem to $(s, t) \mapsto h_x(sv, tw)$. As stated in [9], “The function $h_x(v, w)$ is defined invariantly; therefore, the terms of its Taylor series are invariants of the connection. Such invariants must be expressed in terms of the curvature tensor, torsion tensor, and their derivatives.”

In this work, we will require an approximation at order four [8, Section 8.], although higher-order terms can be found in [10]. For $s, t \in \mathbb{R}$ small enough:

$$ 
\begin{aligned}
h_x(sv, tw) &= \log_x(\exp_{x_{sv}}(\Pi_x^{x_{sv}} tw)) \\
&= sv + tw + \frac{st}{6} R(w, v)(sv + 2tw) \\
&\quad + \frac{st}{24} \left( (\nabla_v R)(w, sv)(5tw + 2sv) \right. \\
&\quad \left. + (\nabla_w R)(tw, v)(sv + 2tw) \right) + q_5(sv, tw),
\end{aligned} 
\quad (1) $$

where $q_5(sv, tw)$ contains homogeneous terms of degree five and higher, whose coefficients can be bounded uniformly in $K$ (as they can be expressed in terms of the curvature and its covariant derivatives). Because we consider in this paper the Levi-Civita connection of a Riemannian manifold, there is no torsion term appearing here. Thus, taking $s, t \sim \frac{1}{n}$ for some $n \in \mathbb{N}$ large enough, $q_5(sv, tw) = O(\frac{1}{n^5})$. To simplify the notations, we write $O(5)$ in that sense. Formally (1) is similar to the BCH formula in Lie groups.

Similarly, [29] introduced the neighboring log and computed its Taylor approximation. It is defined by applying the $\exp_x$ map to small enough $v, w \in T_x M$ to obtain $x_v, x_w$ then computing the log of $x_w$ from $x_v$, and finally parallel transporting this vector back to $x$ (see Fig. 2). This defines $l_x : U_x \to T_x M$ by:

$$ l_x(v, w) = \Pi_{x_v}^x \log_{x_v}(x_w), $$

which relates to the double exponential by solving

$$ w = h_x(v, l_x(v, w)). $$

Using (1), one can solve for the first terms of a Taylor expansion and obtains for $(v, w) \in U_x \subset T_x M \times T_x M$ [29]:

<figure><img src="image_9.png" /></figure>
$$ 
\begin{aligned}
l_x(v, w) &= \prod_{x_v}^{x} \log_{x_v}(x_w) \\
&= w - v + \frac{1}{6} R(v, w)(2w - v) \\
&\quad + \frac{1}{24} \left( (\nabla_v R)(v, w)(3w - 2v) \right. \\
&\qquad \left. + (\nabla_w R)(v, w)(2w - v) \right) + O(5),
\end{aligned} 
\quad (2) $$

where we have implicitly assumed $v, w$ small enough and taken the time variables $s = t = 1$. We will do so in the next section to simplify the notations.

## 2 Schild's Ladder

We now turn to the construction of Schild's ladder and to the analysis of this numerical scheme.

### 2.1 Elementary Construction

The construction to parallel transport $v \in T_x M$ along the geodesic $\gamma$ with $\gamma(0) = x$ and $\dot{\gamma}(0) = w \in T_x M$ (such that $(v, w) \in U_x$) is given by the following steps (Fig. 3):

1. Compute the geodesics from $x$ with initial velocities $v$ and $w$ until time $s = t = 1$ to obtain $x_v$ and $x_w$. These are the sides of the parallelogram.

2. Compute the geodesic between $x_v$ and $x_w$ and the midpoint $m$ of this geodesic, i.e.,

$$ m = \exp_{x_v} \left( \frac{1}{2} \log_{x_v} (x_w) \right). $$

This is the first diagonal of the parallelogram.

3. Compute the geodesic between $x$ and $m$, let $a \in T_x M$ be its initial velocity. Extend it beyond $m$ for the same length as between $x$ and $m$ to obtain $z$, i.e.,

$$ a = \log_x(m); \quad z = \exp_x(2a) = x_{2a}. $$

This is the second diagonal of the parallelogram.

<figure><img src="image_10.png" /><figcaption>**Fig. 3** Construction of a geodesic parallelogram in Schild's ladder, represented in a normal coordinate system centered at $x$. Geodesics are in black and can be identified with tangent vectors (colored) when going through $x$ (Color figure online)</figcaption></figure>

<figure><img src="image_11.png" /></figure>
4. Compute the geodesic between $x_w$ and $z$. Its initial velocity $u^w$ is an approximation of the parallel transport of $v$ along the geodesic from $x$ to $x_w$, i.e.,

$$u^w = \log_{x_w}(x_{2a}).$$

By assuming that there exists a convex neighborhood that contains the entire parallelogram, all the above operations are well defined. In the literature, this construction is then iterated along $\gamma$ without further precision.

## 2.2 Taylor Expansion

We can now reformulate this elementary construction in terms of successive applications of the double exponential and the neighboring logarithm maps. Let $b^v = \log_{x_v}(x_w) \in T_{x_v}M$ be the initial velocity of the geodesic computed in step 2 of the above construction, and $b$ its parallel transport to $x$. Therefore, we have $b = \Pi_{x_v}^x \log_{x_v}(x_w) = l_x(v, w)$. The midpoint is now computed by $m = \exp_{x_v}(\frac{b^v}{2})$, i.e.,

$$m = \exp_{x_v}(\Pi_{x_v}^x \frac{b}{2}). \qquad (3)$$

Thus, step 3 is equivalent to $a = h_x(v, \frac{1}{2}l_x(v, w))$. Combining the Taylor approximations (1) and (2), we obtain an expansion of $2a$. The computations are detailed in Appendix A, and we only report here the third order, meaning that all the terms of the form $\nabla.R(\cdot, \cdot)$ are summarized in the term $O(4)$:

$$2a = w + v + \frac{1}{6}R(v, w)(w - v) + O(4). \qquad (4)$$

We notice that this expression is symmetric in $v$ and $w$, as expected. Furthermore, the deviation from the Euclidean mean of $v, w$ (the parallelogram law) is explicitly due to the curvature. Accounting for this correction term is a key ingredient to reach a quadratic speed of convergence.

Now, in order to compute the error $e_x$ made by this construction to parallel transport $v$, $e_x = \Pi_x^x v - u^w$, we parallel transport it back to $x$: define $u = \Pi_{x_w}^x u^w$ and $e = \Pi_{x_w}^x e_x = v - u \in T_x M$. Now as $u^w = \log_{x_w}(x_{2a})$, we have $u = l_x(w, 2a)$. Combining (2) with (4) (see Appendix A), we obtain the first main result of this paper, a third-order approximation of the Schild's ladder construction:

**Theorem 1** Let $(M, g)$ be a finite-dimensional Riemannian manifold. Let $x \in M$ and $v, w \in T_x M$ sufficiently small. Then, the output $u$ of one step of Schild's ladder parallel transported back to $x$ is given by

$$u = v + \frac{1}{2}R(w, v)v + O(4). \qquad (5)$$

The fourth order and a bound on the remainder are detailed in Appendix A. This theorem shows that Schild's ladder is in fact a second-order approximation of parallel

<figure><img src="image_12.png" /></figure>
transport. Furthermore, this shows that splitting the main geodesic into segments of length $\frac{1}{n}$ and simply iterating this construction $n$ times will in fact sum $n$ error terms of the form $R(\frac{w_i}{n}, v_i)v_i$, hence by linearity of $R$, the error won't necessarily vanish as $n \to \infty$. To ensure convergence, it is necessary to also scale $v$ in each parallelogram. The second main contribution of this paper is to clarify this procedure and to give a proof of convergence of the numerical scheme when scaling both $v$ and $w$. This is detailed in the next subsection.

## 2.3 Numerical Scheme and proof of convergence

With the notations of the previous subsection, let us define schild($x, w, v$) = $u^w \in T_{xw} M$. We now divide the geodesic $\gamma$ into segments of length $\frac{1}{n}$ for some $n \in \mathbb{N}^*$ large enough so that the previous Taylor expansions can be applied to $\frac{w}{n}$ and $\frac{v}{n}$. As mentioned before, $v$ needs to be scaled as $w$. In fact, let $\alpha \ge 1$ and consider the sequence defined by (see Fig. 4)

$$v_0 = v \\ v_{i+1} = n^{\alpha} \cdot \text{schild}(x_i, \frac{w_i}{n}, \frac{v_i}{n^{\alpha}}), \qquad (6)$$

where $x_i = \gamma(\frac{i}{n}) = \exp_X(\frac{i}{n}w)$, $w_i = n \log_{x_i}(x_{i+1}) = \Pi_x^{x_i} w$. We now establish the following result, which ensures convergence of Schild's ladder to the parallel transport of $v$ along $\gamma$ at order at most two.

**Theorem 2** Let $(x_i, v_i, w_i)_{(i \le n)}$ be the sequence defined as above. Then, $\exists \tau > 0, \exists \beta > 0, \exists N \in \mathbb{N}, \forall n > N$,

$$\|v_n - \Pi_x^{x_n} v\| \le \frac{\tau}{n^\alpha} + \frac{\beta}{n^2}.$$

Moreover, $\tau$ is bounded by a bound on the sectional curvature, and $\beta$ by a bound on the covariant derivative of the curvature tensor.

**Proof** To compute the accumulated error, we parallel transport back to $x$ the results $v_n$ of the algorithm after $n$ rungs. The error is then written as the sum of the errors at each rung, parallel transported to $x$ (isometrically):

<figure><img src="image_13.png" /><figcaption>Fig. 4 Schild's ladder, consisting in iterative constructions of geodesic parallelograms</figcaption></figure>

<figure><img src="image_14.png" /></figure>
$$
\begin{equation}
\begin{aligned}
\frac{1}{n^\alpha} \|v_0 - \Pi_{x_n}^x v_n\| &= \left\| \sum_{i=0}^{n-1} \Pi_{x_i}^x \left( \frac{v_i}{n^\alpha} - \Pi_{x_i+1}^{x_i} \frac{v_{i+1}}{n^\alpha} \right) \right\| \\
&\le \sum_{i=0}^{n-1} \left\| \Pi_{x_i}^x \left( \frac{v_i}{n^\alpha} - \Pi_{x_i+1}^{x_i} \frac{v_{i+1}}{n^\alpha} \right) \right\| \le \sum_{i=0}^{n-1} \|e_i\|.
\end{aligned}
\tag{7}
\end{equation}
$$

By (5) and Lemma 4, the error at each step can be written

$$
e_i = \frac{v_i}{n^\alpha} - \Pi_{x_i+1}^{x_i} \frac{v_{i+1}}{n^\alpha} = \frac{1}{2n^{(1+2\alpha)}} R(v_i, w_i)v_i + r_4\left(\frac{w_i}{n}, \frac{v_i}{n^\alpha}\right). \quad (8)
$$

where $r_4$ contains the fourth-order residual terms and is given in Appendix A. We start
as in [21] by assuming that $v_i$ does not grow too much, i.e., there exists $s \in \{1, \dots, n\}$
such that $\forall i < s, \|v_i\| \le 2\|v\|$. We will then show that the control obtained on $v_n$
is tight enough so that when $n$ is large enough, $\forall k \le n, \|v_k\| \le 2\|v\|$. With this
assumption each term of the right-hand side (r.h.s.) can be bounded. As $\alpha \ge 1$, we
apply Lemma 4 (in Appendix A): $\exists \beta > 0$ such that for $n$ large enough:

$$
\forall i \le n, \|r_4(v_i, w_i)\| \le \frac{\beta}{n^{\alpha+3}}.
$$

Similarly, as $\forall i \le n$, $\|w_i\| = \|w\|$ and by assumption $\|v_i\| \le 2\|v\|$, we have

$$
\| R(v_i, w_i) v_i \| \le \| R \|_{\infty} \| v_i \|^2 \| w_i \| \le 4 \| v \|^2 \| R \|_{\infty} \| w \|
$$

Let $\tau = 2\|v\|^2 \|R\|_{\infty} \|w\|$. This gives

$$
\|e_i\| \le \frac{\tau}{n^{(1+2\alpha)}} + \frac{\beta}{n^{\alpha+3}}. \tag{9}
$$

As the r.h.s. does not depend on $n$, it can be plugged into (7) to obtain by summing
for $i = 0, \dots, s-1 \le n$:

$$
\frac{1}{n^{\alpha}} \|v_0 - \Pi_{x_s}^{x} v_s\| \le \sum_{i=0}^{s-1} \|e_i\| \le \frac{\tau}{n^{2\alpha}} + \frac{\beta}{n^{\alpha+2}},
$$

and finally

$$
\|v - \Pi_{x_s}^x v_s\| \le \frac{\tau}{n^\alpha} + \frac{\beta}{n^2}. \tag{10}
$$

Now suppose for contradiction that for arbitrarily large $n$, there exists $k \le n$ such that
$\|v_k\| > 2\|v\|$ and choose $u(n) \in 1, \dots, n$ minimal with this property, i.e., $\|v_{u(n)}\| >$
$2\|v\|$ so that (10) can be used with $s = u(n)$. Then, we have:

$$
\|v\| < \|v\| - \|v_{u(n)}\| \le \|v - \Pi_{x_{u(n)}}^x v_{u(n)}\| \le \frac{\tau}{n^\alpha} + \frac{\beta}{n^2}. \quad (11)
$$

<figure><img src="image_15.png" /></figure>
But the r.h.s. goes to 0 as $n$ goes to infinity, leading to a contradiction. Therefore, for $n$ large enough, $\forall i \le n$, $\|v_i\| \le 2\|v\|$, and the previous control on $\|v - \Pi_{x_s}^x v_s\|$ given by (10) is valid for $s = n$. The result follows by parallel transporting this error to $x_n$, which does not change its norm. $\square$

**Remark 1** The proof allows to grasp the origin of the bounds in Thm. 1. The term $O(n^{-\alpha})$ comes from $R(w, v)v$ in (5), as it is bilinear in $v$, so that the division by $n^\alpha$ is squared in the error, while we only need to multiply by $n^\alpha$ at the final step to recover the parallel transport of $v$. Similarly, as $R(w, v)v$ is linear in $w$, the division by $n$ and summing of $n$ terms does not appear in the error.

On the other hand, the $O(n^{-2})$ comes from $\nabla_w R(w, v)w$ (from Lemma 4). Indeed, this term is invariant by the division/multiplication by $n^\alpha$, unlike other error terms where $v$ appears several times, and the multilinearity w.r.t. $w$ implies that it is a $O(n^{-3})$ that is then summed $n$ times. We therefore use $\alpha = 2$ in practice.

**Remark 2** The bound on the curvature that we use can be related to a bound on the sectional curvature $\kappa$ as follows. Recall

$$ \kappa(X, Y) = \frac{\langle R(Y, X)X, Y \rangle}{\|X\|^2 \|Y\|^2 - \langle X, Y \rangle^2}. \qquad (12) $$

Suppose it is bounded on the compact set $K: \forall X, Y \in \Gamma(K)$, $|\kappa(X, Y)| < \|\kappa\|_\infty$.
Therefore, for two linearly independent $X, Y$, write $\delta_{X,Y} = \frac{\langle X, Y \rangle^2}{\|X\|^2 \|Y\|^2} \in [0, 1)$,

$$ \begin{aligned} |\langle R(Y, X)X, Y \rangle| &= |\kappa(X, Y)|(\|X\|^2 \|Y\|^2 - \langle X, Y \rangle^2) \\ &= |\kappa(X, Y)|\|X\|^2 \|Y\|^2 (1 - \delta_{X,Y}) \\ &\le \|\kappa\|_\infty \|X\|^2 \|Y\|^2. \end{aligned} $$

Now, the linear map $Y \mapsto R(Y, X)X$ is self-adjoint, so if $|\lambda_1(X)| \le \dots \le |\lambda_d(X)|$ are its eigenvalues, we have on the first hand $|\lambda_d(X)| \le \|\kappa\|_\infty \|X\|^2$ and on the second hand $\|R(Y, X)X\| \le |\lambda_d(X)|\|Y\|$. Thus,

$$ \|R(Y, X)X\| \le \|\kappa\|_\infty \|X\|^2 \|Y\|. \qquad (13) $$

We then could have used $\tau = \frac{\delta^3 \|\kappa\|_\infty}{2}$ in the proof of theorem 2.

We notice in this proof and Remark 1 that the terms of the Taylor expansion where $v$ appears several times vanish faster thanks to the arbitrary exponent $\alpha \ge 1$. Together with Lemma 4, this implies that the dominant term is the one where $v$ appears only once (given explicitly in Appendix A), which imposes here a speed of convergence of order two. We will use this key fact to compute the speed of convergence of the other schemes.

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>
## 2.4 Numerical Simulations

In this section, we present numerical simulations that show the convergence bounds in two simple cases: the sphere $S^2$ and the space of $3 \times 3$ symmetric positive-definite matrices $SPD(3)$.

### 2.4.1 The Sphere

We consider the sphere $S^2$ as the subset of unit vectors of $\mathbb{R}^3$, endowed with the canonical ambient metric $\langle \cdot, \cdot \rangle$. It is a Riemannian manifold of constant curvature, and the geodesics and parallel transport are available in closed form (these are given in Appendix C for completeness).

Because the sectional curvature $\kappa$ is constant, the fourth-order terms vanish and Theorem 2 becomes $\|v - \Pi_{x_n}^x v_n\| \le \frac{\tau}{n^\alpha}$ for $\alpha \ge 1$. This is precisely observed in Fig. 5.

Moreover, using Eq. (8), we have at each rung:

$$ \begin{aligned} \langle \Pi_{x_i}^x e_i, w \rangle &= \langle e_i, w_i \rangle = \frac{1}{2n^{(1+2\alpha)}} \langle R(v_i, w_i)v_i, w_i \rangle + O\left(\frac{1}{n^{(4+\alpha)}}\right) \\ &= \frac{-1}{2n^{(1+2\alpha)}} \kappa(v_i, w_i) = \frac{-1}{2n^{(1+2\alpha)}} + O\left(\frac{1}{n^{(4+\alpha)}}\right). \end{aligned} $$

Therefore, by summing as in the proof of Theorem 2 we obtain:

$$ \langle v_n - \Pi_x^{x_n} v, w_n \rangle = \frac{1}{2n^\alpha} + O\left(\frac{1}{n^3}\right). $$

Thus, the projection of the error onto $w$, which we call longitudinal error, allows to verify our theory: we can measure the slope of the decay with respect to $n^{-\alpha}$ and it should equal $\frac{1}{2}$. This is very precisely what we observe (Fig. 5, right).

<figure><img src="image_18.png" /><figcaption>Fig. 5 Error of the parallel transport of $v$ along the geodesic with initial velocity $w$ where $v$ and $w$ are orthonormal. Accordingly with our main result, any speed of convergence $1 \le \alpha \le 2$ can be reached. Left: absolute error w.r.t. $\frac{1}{n}$. Right: longitudinal error w.r.t. $\frac{1}{n^\alpha}$: on $S^2$ this is a straight line with slope $\frac{1}{2}$.</figcaption></figure>

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 6 Error of the parallel transport of v along the geodesic with initial velocity w where v and w are orthonormal. Accordingly with our main result, a quadratic speed of convergence $\alpha = 2$ can be reached. Left: absolute error w.r.t. $\frac{1}{n^2}$. Right: longitudinal error w.r.t. $\frac{1}{n^2}$: with the AI metric, this is a straight line with negative slope</figcaption></figure>

## 2.4.2 SPD

We now consider SPD(3) endowed with the affine-invariant metric [30]. Again, the formulas for the exp, log and parallel transport maps are available in closed form (Appendix D). Because the sectional curvature is always non-positive, the slope of the longitudinal error is negative, and this is observed in Fig. 6.

## 3 Variations of SL

We now turn to variations of SL. First, we revisit the Fanning scheme using the neighboring log and show how close to SL it in fact is. Then, we examine the pole ladder and introduce the averaged SL. Furthermore, we introduce infinitesimal schemes, where the exp and log maps are replaced by one step of a numerical integration scheme. We give convergence results and simulations for the PL in this setting. This allows to probe the convergence bounds in settings where the underlying space is not symmetric, and therefore observe the impact of a nonzero $\nabla R$.

### 3.1 The Fanning Scheme, Revisited

The Fanning scheme [22] leverages an identity given in [34] between Jacobi fields and parallel transport. As proved in [22], it converges linearly when dividing the main geodesic into $n$ segments and using second-order integration schemes of the Hamiltonian equations. We recall the algorithm here with our notations and compare it to SL. For $v, w \in T_x \mathcal{M}$, define the Jacobi field (JF) along $\gamma: t \mapsto \exp_x(tw)$ at time $h$ small enough by:

$$J_{\gamma(t)}^{v}(h) = \left. \frac{\partial}{\partial \epsilon} \right|_{\epsilon=0} \exp_{\gamma(t)} (h(\dot{\gamma}(t) + \epsilon v)).$$

<figure><img src="image_21.png" /></figure>
**Fig. 7** Elementary construction of the Fanning scheme. For $\epsilon = \frac{1}{n}$, it is similar to that of SL (with $\alpha = 2$) except that the midpoint is not computed but $w + \epsilon v$ is used instead of $\log_x(m)$

<figure><img src="image_22.png" /></figure>

One can then show that $J_{\gamma(t)}^v(h) = h\Pi_{\gamma,t}^{t+h}v + O(h^2)$ [34]. The Fanning scheme consists in computing perturbed geodesics, i.e., with initial velocity $w + \epsilon v$ for some $\epsilon > 0$ and to approximate the associated JF by

$$J_{\gamma(t)}^{v}(h) = \frac{\exp_{\gamma(t)}(h(\dot{\gamma}(t) + \epsilon v)) - \exp_{\gamma(t)}(h\dot{\gamma}(t))}{\epsilon}, \quad (14)$$

and the parallel transport of $v$ along $\gamma$ between 0 and $h$ by

$$u^w = \frac{\exp_x(h(w + \epsilon v)) - \exp_x(hw)}{h\epsilon}.$$

This defines an elementary construction that is iterated $n$ times with time-steps $h = \frac{1}{n}$, and [22] showed that using $\epsilon = h$ yields the desired convergence. Let us consider that the finite difference approximation of (14) is in fact a first-order approximation of the log of $\exp_x(h(w + \epsilon v))$ from $\exp_x(hw)$. Now this scheme is similar to SL, except that it uses $2a = w + \epsilon v$ instead of (4), i.e., the correction term according to curvature is not accounted for (compare Figs. 3 and 7). When using $h = \epsilon = \frac{1}{n}$, this corresponds to $\alpha = 2$ in our analysis of the sequence defined by 6. Similarly, we can compute a third order Taylor approximation of the error made at each step of the FS. To do so, introduce as in Sect. 2.2 $u = \Pi_{xhw}^x u^w$. Then,

$$\begin{aligned} u &= l_x(hw, h(w + \epsilon v)) \\ &= h\epsilon v + \frac{h^2}{6}R(w, w + \epsilon v)(2h(w + \epsilon v) - hw) + O(h^4\epsilon) \\ &= h\epsilon v + \frac{h^3\epsilon}{6}R(w, v)(w + 2\epsilon v) + O(h^4\epsilon). \end{aligned}$$

But as $\epsilon = h$, the last term including $v$ disappears in the $O(h^4\epsilon)$ and we obtain:

$$\frac{u}{h\epsilon} = v + \frac{h^2}{6}R(w, v)w + O(h^3). \quad (15)$$

We notice that unlike the expression given by Theorem 1, this approximation of FS contains a term linear in $v$ (i.e., where $v$ appears only once), and bilinear in $w$, so that when applied to $h = \frac{1}{n}$, summing error terms yields a global error of order $\frac{1}{n}$. This shows that a better speed cannot be reached. Moreover, we recognize the coefficient $\frac{1}{6n}$ computed explicitly on the sphere in [22, sec. 2.3]. As in the previous section, define $u^w = \text{fs}(x, w, v)$ the result of the FS construction, and the sequence

<figure><img src="image_23.png" /></figure>
$$v_0 = v,$$
$$v_{i+1} = n^2 \cdot \text{fs} \left( x_i, \frac{w_i}{n}, \frac{v_i}{n^2} \right), \qquad (16)$$

where $x_i = \gamma(\frac{i}{n}) = \exp_x(\frac{i}{n}w)$, $w_i = n \log_{x_i}(x_{i+1}) = \Pi_x^{x_i} w$. Then, it is straightforward to reproduce the proof of thm. 2 with (15) instead of (5) to show that $\exists \beta > 0, \exists N \in \mathbb{N}, \forall n > N$,

$$\|v_n - \Pi_x^{x_n} v\| \le \frac{\beta}{n}.$$

This corroborates the result of [22], although at this point geodesics are assumed to be available exactly. An improvement of the FS is to use a second-order approximation of the Jacobi field by computing two perturbed geodesics, but this does not change the precision of the approximation of the parallel transport. We implemented this scheme with this improvement, and the result is compared with ladder methods in paragraph 3.4.3.

## 3.2 Averaged SL

Focusing on Eq. (5) of Theorem 1, we notice that the error term is even with respect to $v$. We therefore propose to average the results of a step of SL applied to $v$ and to $-v$ to obtain $u^{(+)}$ and $u^{(-)}$. Using (32) of Appendix A, we have:

$$u^{(+)} = v + \frac{1}{2}R(w, v)v + \frac{7}{48}\left(\left(\nabla_v R\right)(w, v)v + \left(\nabla_w R\right)(w, v)\left(\frac{8}{7}v + w\right)\right) + O(5), \qquad (17)$$

$$u^{(-)} = -v + \frac{1}{2}R(w, v)v - \frac{7}{48}\left(\left(\nabla_v R\right)(w, v)v + \left(\nabla_w R\right)(w, v)\left(w - \frac{8}{7}v\right)\right) + O(5). \qquad (18)$$

Hence

$$\frac{u^{(+)} - u^{(-)}}{2} = v + \frac{7}{48} \left( (\nabla_v R)(w, v)v + (\nabla_w R)(w, v)w \right) + O(5). \qquad (19)$$

This scheme thus allows to cancel out the third-order term in an elementary construction. However, when iterating this scheme, the dominant term will be $(\nabla_w R)(w, v)w$, so that the speed of convergence remains quadratic anyway. Being also heavier to compute, this scheme has thus very little practical advantage to offer, while the pole ladder, detailed in the next section, present both advantages of being exact at the third order, and being much cheaper to compute.

<figure><img src="image_24.png" /></figure>
## 3.3 Pole Ladder

The pole ladder was introduced by [20] to reduce the computational burden of SL. They showed that at the first order, the constructions were equivalent. [28] latter gave a Taylor approximation of the elementary construction showing that each rung of pole ladder is a third order approximation of parallel transport, and additionally showed that PL is exact in affine (hence in Riemannian) locally symmetric spaces (the proof is reproduced in [11]). We first present the scheme and derive the Taylor approximation using the neighboring log.

### 3.3.1 The Elementary Construction

The construction to parallel transport $v \in T_x M$ along the geodesic $\gamma$ with $\gamma(0) = x$ and $\dot{\gamma}(0) = w \in T_x M$ (such that $(v, w) \in U_x$) is given by the following steps (see Fig. 8):

1. Compute the geodesics from $x$ with initial velocities $v$ and $w$ until time $s = t = 1$ to obtain $x_v$ and $x_w$, and $t = \frac{1}{2}$ to obtain the midpoint $m$. The main geodesic $\gamma$ is a diagonal of the parallelogram.

2. Compute the geodesic between $x_v$ and $m$, let $a \in T_m M$ be its initial velocity. Extend it beyond $m$ for the same length as between $x_v$ and $m$ to obtain $z$, i.e.,

$$a = \log_m(x_v); \quad z = \exp_m(-a) = m_{-a}.$$

This is the second diagonal of the parallelogram.

3. Compute the geodesic between $x_w$ and $z$. The opposite initial velocity $u^w$ is an approximation of the parallel transport of $v$ along the geodesic from $x$ to $x_w$, i.e.,

$$u^w = -\log_{x_w}(m_{-a}).$$

By assuming that there exists a convex neighborhood that contains the entire parallelogram, all the above operations are well defined.

### 3.3.2 Taylor Approximation

We now introduce new notations and center the construction at the midpoint $m$ (see Fig. 8, right). Both $v$, $w$ are now tangent vectors at $m$, and $v$ is the parallel transport

<figure><img src="image_25.png" /><figcaption>Fig. 8 Elementary construction of the pole ladder with the previous notations (left), and new notations (right), in a normal coordinate system at $m$</figcaption></figure>

<figure><img src="image_26.png" /></figure>
of the vectors $v^{-w}$ that we wish to transport between $m_{-w}$ and $m_w$. With these new notations,

$$u = \Pi_{m_w}^m u^w = -l_m(w, -a)$$

where $a = h_m(-w, v)$. Therefore, using (1) and (2), we obtain the following (see Appendix B for the computations)

**Theorem 3** Let $(\mathcal{M}, g)$ be a finite dimensional Riemannian manifold. Let $m \in \mathcal{M}$ and $v, w \in T_m\mathcal{M}$ sufficiently small. Then, the output $u$ of one step of the pole ladder parallel transported back to $m$ is given by

$$u = v + \frac{1}{12} \left( (\nabla_w R)(w, v)(5v - w) + (\nabla_v R)(w, v)(2v - w) \right) + O(5). \quad (20)$$

We notice that an elementary construction of pole ladder is more precise than that of Schild's ladder in the sense that there is no third order term $R(\cdot, \cdot)$. Using $n^{-1}$ to scale $w$ and $n^{-\alpha}$ to scale $v$ as in SL, the term $(\nabla_w R)(w, v)w$ is the dominant term when iterating this construction. Indeed, as it is linear in $v$, the scaling does not influence the convergence speed as one needs to multiply the result by $n^{\alpha}$ at the final rung of the ladder to recover the parallel transport of $v$. The other terms are multilinear in $v$, so if $\alpha > 1$, they are small compared to $\frac{1}{n^3}(\nabla_w R)(w, v)w$. Summing $n$ terms of the form $\frac{1}{n^3}(\nabla_w R)(w, v)w$ we thus obtain a quadratic speed of convergence, as for Schild's ladder, for any $\alpha \ge 1$. Note that using $1 < \alpha \le 2$ will reduce the impact of the other terms in the error without changing the overall order of convergence. We henceforth use $\alpha = 1$.

As in Sect. 3.2, one could think of applying pole ladder to linear combinations of $v, w$ to achieve a better performance, but this strategy is doomed as no linear combination can cancel out the dominant term $(\nabla_w R)(w, v)w$.

As in the previous section, define $u^w = \text{pole}(m, w, v)$ the result of the PL construction (with the new notations), and the sequence

$$\begin{aligned} v_0 &= v \\ v_{i+1} &= n \cdot \text{pole}(m_i, \frac{w_i}{n}, \frac{v_i}{n}), \end{aligned} \quad (21)$$

where $m_i = \gamma(\frac{2i+1}{2n}) = \exp_m(\frac{2i}{n}w)$, $w_i = \dot{\gamma}(\frac{2i+1}{2n}) = \frac{n}{2}\log_{m_i}(m_{i+1}) = \Pi_m^{m_i}w$. Then it is straightforward to reproduce the proof of thm. 2 with (20) instead of (5) to show that $\exists \beta > 0, \exists N \in \mathbb{N}, \forall n > N$,

$$\|v_n - \Pi_m^{m_n} v\| \le \frac{\beta}{n^2}, \quad (22)$$

and $\beta$ can be bounded by the covariant derivative of the curvature.

This proves the convergence of the scheme with the same speed as SL. Note however that when iterating the scheme, unlike SL, the $v_i$'s and the geodesics that correspond to the rungs of the ladder need not being computed explicitly, except at the last step. This is shown in Fig. 9, bottom row. This greatly reduces the computational burden of the PL over SL.

<figure><img src="image_27.png" /></figure>
<figure><img src="image_28.png" /><figcaption>**Fig. 9** The pole ladder, consisting in iterations of elementary constructions (top). It can be simplified to lighten the computational burden of computing geodesics (bottom), only one log and one exp need to be computed from $m_i$ at each iteration</figcaption></figure>

## 3.4 Infinitesimal Schemes with Geodesic Approximations

When geodesics are not available in closed form, we replace the exp map by a fourth-order numerical scheme (e.g., Runge-Kutta (RK)), and the log is obtained by gradient descent over the initial condition of exp. It turns out that only one step of the numerical scheme is sufficient to ensure convergence and keeps the computational complexity reasonable. As only one step of the integration schemes is performed, we are no longer computing geodesic parallelograms, but infinitesimal ones, and thus refer to this variant as *infinitesimal scheme*. We here detail the proof for the PL, but it could be applied to SL as well as the other variants.

**Remark 3** While in this section we use the geodesic equation written with Christoffel symbols, any formulation of the geodesic equation (such as the Euler-Poincaré equation on Lie groups, or the Hamilton Jacobi equations) can be used. In practice, as the number of Christoffel symbols is $n^3$ for an $n$-dimensional manifold, these are rarely computed. This prohibits the direct integration of the parallel transport equation written with Christoffel symbols. On the other hand, as there is a wide literature on variational problems and symplectic integrators, these methods are often preferred to compute geodesics. The infinitesimal ladder schemes thus allow to leverage state-of-the art geodesic solvers for the purpose of computing parallel transport.

However, in some cases the parallel transport equation can be written with known quantities, and a direct integration with high-order integration scheme will then provide faster solutions to compute parallel transport. This is explored, e.g., on Lie groups in [13] and on Kendall shape spaces in [12].

Furthermore, infinitesimal schemes are used in computational anatomy, on approximations of the infinite-dimensional group of diffeomorphisms [20]. This section thus provides a proof that this approach is valid in finite dimension. Finally, infinitesimal ladder schemes may also be useful in the context of discrete geodesics [1].
### 3.4.1 Derivation and Proof

More precisely, consider the geodesic equation written as a first order equation of two variables in a global chart $\Phi$ of $M$, that defines for any $p \in M$ a basis of $T_pM$, written

$$B_x^\Phi = \left. \frac{\partial}{\partial x^i} \right|_x, \quad i = 1, \dots, d$$

$$\left\{ \begin{aligned} \dot{x}^k(t) &= v^k(t), \\ \dot{v}^k(t) &= -\Gamma_{ij}^k v^i(t) v^j(t), \end{aligned} \right. \qquad (23)$$

where $\Gamma_{ij}^k$ are the Christoffel symbols and we use Einstein summation convention. Let $rk: TM \times \mathbb{R}_+ \to TM$ be the map that performs one step of a fourth-order numerical scheme (e.g., RK4), i.e., it takes as input $(x(t), v(t)), h)$ and returns an approximation of $(x(t+h), v(t+h))$ when $(x, v)$ is solution of the system (23). In our case, the step-size $h = \frac{1}{n}$ is used. By fourth order, we mean that we have the following local truncation error relative to the 2-norm of the global coordinate chart $\Phi$:

$$\|x(t+h) - rk_1(x(t), v(t), h)\|_2 \le \tau_1 h^5,$$

and global accumulated error after $n$ steps with step size $h = \frac{1}{n}$:

$$\|x(1) - \tilde{x}_n\|_2 \le \frac{\tau_1}{n^4},$$

where by $rk_1$ we mean the projection on the first variable, and $\tau_1 > 0$ is a constant.

This means that any point on the geodesic is approximated with error $O(\frac{1}{n^4})$. As we are working in a compact domain, and all the norms are equivalent, the previous bounds can be expressed in Riemannian distance $d$:

$$d(x(t+h), rk_1(x(t), v(t), h)) \le \tau_2 h^5 \qquad (24)$$

$$d(x(1), \tilde{x}_n) \le \frac{\tau_2}{n^4}. \qquad (25)$$

Furthermore, the inverse of $v \mapsto rk_1(x, v, h)$ can be computed for any $x$ by gradient descent, and we assume<fnref n="3" /> that any desired accuracy can be reached. We write $\tilde{v} = rk_x^{-1}(y)$ for the optimal $v$ such that $d(rk(x, v, \frac{1}{n}), y) = o(\frac{1}{n^5})$ and assume that $\|rk_x^{-1}(y) - \log_x(y)\| \le \frac{\tau_3}{n^4}$. Note that the step size does not appear explicitly in the notation $rk^{-1}$, but $h = \frac{1}{n}$ is always used. As in practice, $x \in M$, $v, w \in T_x M$ are given, the initial $w_0$ is tangent at $x$, and then for $i > 0$, $w_i$ will be tangent at $m_i$ as in (21), but $w_i$ maps $m_i$ to $m_{i+1}$ and this must correspond to one step of size $\frac{1}{n}$ in the discrete scheme, so there is a factor two that differs from the definition of $w_i$ in (21). Now define also $\tilde{m}_i$, $\tilde{w}_i$, $\tilde{v}_i$, using $rk$ instead of exp and log, that is (see Fig. 10):

<figure><img src="image_29.png" /></figure>

<footnote n="3">³ see Remark 4 thereafter about the validity of this hypothesis</footnote>
<figure><img src="image_30.png" /><figcaption>**Fig. 10** Representation of one iteration of the infinitesimal PL scheme. The exact geodesics and logs are plain lines, while their approximations with a numerical scheme are dashed</figcaption></figure>

$$ 
\begin{aligned}
m_0 &= \exp_x\left(\frac{w}{2n}\right) & z_0 &= \exp_x\left(\frac{v}{n}\right) \\
m_{i+1} &= \exp_{m_i}\left(\frac{w_i}{n}\right) & z_{i+1} &= \exp_{m_i}\left(-\log_{m_i}(z_i)\right) \\
\tilde{m}_0, \frac{\tilde{w}_0}{2} &= rk\left(x, \frac{w}{2}, \frac{1}{n}\right) & \tilde{z}_0 &= rk_1\left(x, v, \frac{1}{n}\right) \\
\tilde{m}_{i+1}, \tilde{w}_{i+1} &= rk\left(\tilde{m}_i, \tilde{w}_i, \frac{1}{n}\right) & \tilde{z}_{i+1} &= rk_1\left(\tilde{m}_i, -rk_{\tilde{m}_i}^{-1}(\tilde{z}_i), \frac{1}{n}\right).
\end{aligned}
$$

Finally, let $x_n = \exp(w)$, $v_n = (-1)^n n \log_{x_n}(z_n)$ and their approximations $\tilde{x}_n = rk(\tilde{m}_n, \frac{\tilde{w}_n}{2}, \frac{1}{n})$ and $\tilde{v}_n = n \cdot rk_{\tilde{x}_n}^{-1}(\tilde{z}_n)$.

We will prove that this approximate sequence converges to the true parallel transport of $v$:

**Theorem 4** Let $(\tilde{v}_n)_n$ be the sequence defined above, corresponding to the result of the pole ladder with approximate geodesics computed by a fourth-order method in a global chart. Then,

$$ \| \Pi_x^{x_n} v - \tilde{v}_n \| = O\left(\frac{1}{n^2}\right). $$

**Proof** It suffices to show that there exists $\beta' > 0$ such that for $n$ large enough $\|\tilde{v}_n - v_n\| \le \frac{\beta'}{n^2}$. Indeed, by (22), for $n$ large enough:

$$ \| \Pi_x^{x_n} v - \tilde{v}_n \| \le \| \Pi_x^{x_n} v - v_n \| + \| \tilde{v}_n - v_n \| \le \frac{\beta + \beta'}{n^2}. $$

The approximations made when computing the geodesics with a numerical scheme accumulate at three steps: (1) the RK scheme compared to the true geodesic, this is controlled with the above hypotheses, (2) the distance between the results of the exp map when both the footpoint and the input vector vary a little, this is handled in Lemma 1, and (3) the difference between the results of the log map when both the foot-point and the input vary a little, this is similar to (2) and is handled in Lemma 2. See Fig. 11 for a visual intuition of those lemmas.

**Lemma 1** $\forall x, \tilde{x} \in M$ such that $x$ and $\tilde{x}$ are close enough, for all $v \in T_x M, \tilde{v} \in T_{\tilde{v}} M$ such that both $v$ and $\delta v = \Pi_{\tilde{x}}^x \tilde{v} - v$ are small enough, we have:

$$ d(\exp_x(v), \exp_{\tilde{x}}(\tilde{v})) \le d(x, \tilde{x}) + \|\Pi_{\tilde{x}}^x \tilde{v} - v\|. $$

<figure><img src="image_31.png" /></figure>
**Fig. 11** Visualization for the two lemmas: Lemma 1 seeks to bound the norm of $d$ given $\delta x$ and $\delta v$, while Lemma 2 bounds the norm of $\delta v$ given $d$ and $\delta x$. Here, $h$ is the shorthand for $h_x(\delta x, v + \delta v)$ and $d$ for $\log_{x_v}(\tilde{x}_{\tilde{v}})$

<figure><img src="image_32.png" /></figure>

**Proof** Let $\delta x = \log_x(\tilde{x})$, $x_v = \exp_x(v)$ and $\tilde{x}_{\tilde{v}} = \exp_{\tilde{x}}(\tilde{v})$. By the definition of the double exp and $\delta v$, $\tilde{x}_{\tilde{v}} = \exp_x(h_x(\delta x, v + \delta v))$. Then, by the definition of the neighboring log, we have

$$ \log_{x_v}(\tilde{x}_{\tilde{v}}) = \log_{x_v}(\exp_x(h_x(\delta x, v + \delta v))) = \Pi_x^{x_v} l_x(v, h_x(\delta_x, v + \delta v)). \quad (26) $$

Using the Taylor approximations (1) and (2) truncated at the order of $\|v\|$, we obtain

$$ \Pi_{x_v}^x \log_{x_v}(\tilde{x}_{\tilde{v}}) = l_x(v, \delta x + v + \delta v) = \delta x + \delta v. $$

So that the Riemannian distance $d$ between $x_v$ and $\tilde{x}_{\tilde{v}}$ is

$$ d = \| \log_{x_v}(\tilde{x}_{\tilde{v}}) \| \le \| \delta x \| + \| \delta v \|, $$

and the result follows.

□

**Lemma 2** $\forall x, \tilde{x}, z \in M$ close enough to one another, we have in the metric norm

$$ \| \log_x(z) - \Pi_x^x \log_{\tilde{x}}(\tilde{z}) \| \le d(x, \tilde{x}) + d(z, \tilde{z}). \quad (27) $$

**Proof** The proof is similar to that of Lemma 1 except that this time it is the norm of $\|\delta v\|$ that needs to be bounded by $d(x_v, \tilde{x}_{\tilde{v}})$. □

Now, we first show that the sequence $(\delta_i = d(\tilde{z}_i, z_i))_i$ verifies an inductive relation, such that it is bounded by $\frac{1}{n^3}$, and then we will use Lemma 2 to conclude. We first write

$$ \| \log_{m_i}(z_i) - r k_{\tilde{m}_i}^{-1}(\tilde{z}_i) \| \le \| \log_{m_i}(z_i) - \log_{\tilde{m}_i}(\tilde{z}_i) \| + \| \log_{\tilde{m}_i}(\tilde{z}_i) - r k_{\tilde{m}_i}^{-1}(\tilde{z}_i) \|. $$

The second term on the r.h.s. corresponds to the approximation of the log by gradient descent and is bounded by hypothesis. For the first term, $d(m_i, \tilde{m}_i) \le \frac{\tau_2}{n^4}$ by hypothesis on the scheme (25). Suppose for a proof by induction on $i$ that $\delta_i = d(\tilde{z}_i, z_i) \le \frac{\tau_2}{n}$ and $d(m_i, z_i) \le \frac{\|w\|+2\|v\|}{n}$. This is verified for $i = 0$ and allows to apply Lemma 2 for $n$ large enough, so

$$ \| \log_{m_i}(z_i) - r k_{\tilde{m}_i}^{-1}(\tilde{z}_i) \| \le d(\tilde{m}_i, m_i) + \delta_i + \frac{\tau_3}{n^4}. \quad (28) $$

<figure><img src="image_33.png" /></figure>
Furthermore, by Lemma 1 applied to $v = -\log_{m_i}(z_i)$ and $\tilde{v} = -rk_{\tilde{m}_i}^{-1}(\tilde{z}_i)$, which are sufficiently small by the induction hypothesis,

$$
\begin{aligned}
\delta_{i+1} &= d(z_{i+1}, \tilde{z}_{i+1}) = d(rk_1(\tilde{m}_i, -rk_{\tilde{m}_i}^{-1}(\tilde{z}_i), \frac{1}{n}), \exp_{m_i}(-\log_{m_i}(z_i))) \\
&\le d(rk_1(\tilde{m}_i, -rk_{\tilde{m}_i}^{-1}(\tilde{z}_i), \frac{1}{n}), \exp_{\tilde{m}_i}(-rk_{\tilde{m}_i}^{-1}(\tilde{z}_i))) \\
&\quad + d(\exp_{\tilde{m}_i}(-rk_{\tilde{m}_i}^{-1}(\tilde{z}_i)), \exp_{m_i}(-\log_{m_i}(z_i))) \\
&\le \frac{\tau_2}{n^5} + d(m_i, \tilde{m}_i) + \|\log_{m_i}(z_i) - rk_{\tilde{m}_i}^{-1}(\tilde{z}_i)\|.
\end{aligned}
$$

And combining the two results, we obtain $\delta_{i+1} \le \frac{(2n+1)\tau_2+n\tau_3}{n^5} + \delta_i$, which completes the induction for $\delta_i$. For $d(z_i, m_i)$ we have:

$$d(z_{i+1}, m_{i+1}) = d(z_i, m_{i+1}) \le d(x_i, z_i) + d(x_i, m_i) = \|\frac{v_i}{n}\| + \|\frac{w_i}{n}\|. \quad (29)$$

In the section on SL, we proved that $\|v_i\| \le 2\|v\|$ for $n$ large enough. This applies here as well, and the fact that $w_i$ is the parallel transport of $w$ completes the proof by induction.

By summing the terms for $i = 0, \dots, n-1$, and using $\delta_0 \le \frac{\tau_2}{n^5}$ by (24), we thus have $\delta_n \le \frac{\tau_3+2\tau_2}{n^3} + \frac{n+1}{n^5}\tau_2$. We finally apply Lemma 2 to the scaled $v_n$, $\tilde{v}_n$, as $x_n$ and $\tilde{x}_n$ are close enough:

$$
\begin{aligned}
\frac{1}{n} \|v_n - \tilde{v}_n\| &= \| \log_{x_n}(z_n) - rk_{\tilde{x}_n}^{-1}(\tilde{z}_n) \| \\
&\le \| \log_{x_n}(z_n) - \log_{\tilde{x}_n}(\tilde{z}_n) \| + \| \log_{\tilde{x}_n}(\tilde{z}_n) - rk_{\tilde{x}_n}^{-1}(\tilde{z}_n) \| \\
&\le d(x_n, \tilde{x}_n) + \delta_n + \frac{\tau_2}{n^5}.
\end{aligned}
$$

So that for $n$ large enough $\|\tilde{v}_n - v_n\| \le \frac{\beta'}{n^2}$ for some $\beta' > 0$. $\square$

**Remark 4** To justify the hypothesis on $rk^{-1}$, namely $\|rk_{x_n}^{-1}(y) - \log_x(y)\|_2 \le \frac{\tau}{n^5}$, we consider the problem (P) which corresponds to an energy minimization. Working in a convex neighborhood, it admits a unique minimizer $v^*$:

$$\min \frac{1}{2} \|v\|_2^2 \quad \text{s.t.} \quad \exp_x(v) = y. \qquad (P)$$

The constraint can be written with Lagrange multipliers,

$$\min \frac{1}{2} \|v\|_2^2 + \lambda \| \exp_x(v) - y \|_2^2. \qquad (P')$$

<figure><img src="image_34.png" /></figure>
Now as the exp map is approximated by $rk_1$ at order 5 locally, writing for any $v$ and $h$ small enough, $\| \exp_x(hv) - y \|_2 \le \| rk_1(x, v, h) - y \|_2 + \tau_1 h^5$, ($P'$) is equivalent to

$$ \min \|rk_1(x, v, h) - y\|_2^2 + \frac{1}{2\lambda} \|v\|_2^2, \qquad (Q) $$

which is solved by gradient descent (GD) until a convergence tolerance $\epsilon \le h^5$ is reached.

### 3.4.2 Numerical Simulations

It is not relevant to compare the PL with SL on spheres or SPD matrices as in the previous section, as these spaces are symmetric and thus the PL is exact [11]. We therefore focus on the Lie group of isometries of $\mathbb{R}^3$, endowed with a left-invariant metric $g$ with the diagonal matrix at identity:

$$ G = \operatorname{diag}(1, 1, 1, \beta, 1, 1), \qquad (30) $$

where the first three coordinates correspond to the basis of the Lie algebra of the group of rotations, while the last three correspond to the translation part, and $\beta > 0$ is a coefficient of anisotropy. These metrics were considered in [35] in relation with kinematics, and visualization of the geodesics was provided, but no result on the curvature was given. Following [23], we compute explicitly the covariant derivative of the curvature and deduce (proof given in Appendix E)

**Lemma 3** ($SE(3), g$) is locally symmetric, i.e., $\nabla R = 0$, if and only if $\beta = 1$.

This allows to observe the impact of $\nabla R$ on the convergence of the PL. In the case where $\beta = 1$, the Riemannian manifold ($SE(3), g$) corresponds to the direct product $SO(3) \times \mathbb{R}^3$ with the left-invariant product metric formed by the canonical bi-invariant metric on the group of rotations $SO(3)$ and the canonical inner product of $\mathbb{R}^3$. Therefore, the geodesics can be computed in closed form. Note that the left-invariance refers to the group law which encompasses a semi-direct action of the rotations on the translations. When $\beta \ne 1$ however, geodesics are no longer available in closed form and the infinitesimal scheme is used, with the geodesic equations computed numerically (detailed in Appendix E). In this case, we compute the pole ladder with an increasing number of steps and then use the most accurate computation as reference to measure the empirical error.

Results are displayed in Fig. 12. Firstly, we observe very precisely the quadratic convergence of the infinitesimal scheme, as straight lines are obtained when plotting the error against $\frac{1}{n^2}$. Secondly, we see how the slope varies with $\beta$: accordingly with our results, it cancels for $\beta = 1$ and increases as $\beta$ grows.

Finally, for completeness, we compare Schild's ladder and the pole ladder in this context. We cannot choose two basis vectors $v, w$ (that don't change when $\beta$ changes) such that $R(w, v)v \ne 0$ when $\beta = 1$ and $\nabla_w R(w, v)w \ne 0$ when $\beta \ne 1$. Indeed, the former condition implies that $v, w$ are infinitesimal rotations, which implies

<figure><img src="image_35.png" /></figure>
<figure><img src="image_36.png" /><figcaption>**Fig. 12** Error of the parallel transport of $v$ along the geodesic with initial velocity $w$ where $v$ and $w$ are orthonormal basis vectors (they stay the same for all values of $\beta$). In accordance with our main result, a quadratic speed of convergence is reached, and the speed depends on the asymmetry of the space, itself induced by the anisotropy $\beta$ of the metric. Left: absolute error w.r.t. $\frac{1}{n^2}$; these are straight lines with slopes growing with $\beta$. Right: absolute error w.r.t. $\frac{1}{n}$ showing the quadratic convergence</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>**Fig. 13** Left: comparison of the speed of convergence of SL and PL for two orthonormal basis vectors of $T_x SE(3)$ depending on the coefficient of anisotropy $\beta$ of the metric. PL converges faster and is much more stable. Right: comparison of SL, the PL and the FS in the same setting with fixed $\beta$. (Note that absolute and relative errors are the same here as $v$ has unit norm)</figcaption></figure>

$\forall \beta, \nabla_w R(w, v)w = 0$. Therefore, we choose $v, w$ such that the latter condition is verified, so that the SL error is of order four in our example (but it is not exact), and it cannot be distinguished from PL, but this is only a particular case. The results are shown in Fig. 13. Note that due to the larger number of operations required for SL and as smaller quantities are involved as $\alpha = 2$, our implementation is less stable and diverges when $n$ grows too much (~ 50). As expected when $\beta > 1$, the speeds of convergence are of the same order for PL and SL, but the multiplicative constant is smaller for PL.

We also compare ladder methods to the Fanning scheme, and as expected the quadratic speed of convergence reached by ladder methods yields a far better accuracy even for small $n$. We now compare the infinitesimal SL, PL and the FS in terms of computational cost.
### 3.4.3 Remarks on Complexity

At initialization, ladder methods require to compute $x_v$, with one call to the numerical scheme $rk$ (e.g., Runge-Kutta). Then, the main geodesic needs to be computed, for SL and FS it requires $n$ calls to $rk$. For PL, we only take a half step at initialization, to compute the first midpoint, then compute all the $m_i$s so that one final call to $rk$ is necessary to compute the final point of the geodesic $x_n$, totaling $n + 1$ calls. Then, at every iteration, considering given $m_i, z_i$, one needs to compute an inversion of $rk$, and then shoot with twice (or minus for PL) the result. In practice, the inversion of $rk$ by gradient descent (i.e., shooting) converges in about five or six iterations, each requiring one call to $rk$. This operation thus requires less than ten calls to $rk$. Additionally, for SL only, the midpoint needs to be computed by one inversion and one call to $rk$, thus adding ten calls to the total. At the final step, another inversion needs to be computed, adding less than ten calls. Therefore, SL requires $21n+11$ calls to $rk$, the PL $11(n+1)$. In contrast, the FS only requires to compute one (or two) perturbed geodesic and finite differences instead of the approximation of the log at every step. It therefore require $3n$ calls to $rk$.

Moreover, as the FS is intrinsically a first-order scheme, a second-order $rk$ step is sufficient to guarantee the convergence, thus requiring only two calls to the geodesic equation—the most significantly expensive computation at each iteration—while ladder methods require a fourth-order scheme, which is twice as expensive. However, this additional cost is quickly balanced as only $O(\sqrt{n})$ steps are needed to reach a given accuracy, where $O(n)$ are required for the FS. Comparing the values of $n \mapsto 4 * \frac{11}{\sqrt{n+1}}$ and $n \mapsto 2 * \frac{3}{n}$, we conclude that PL is the cheapest option as soon as we want to achieve a relative accuracy better than $2 \cdot 10^{-2}$. Note that this does not take into account the constants $\beta$, that may differ, but the previous numerical simulations show that these estimates are valid: a regression on the cost for the PL gives $y = 41n + 45$. Finally, in the experiment of Fig. 13, the PL with $n = 6$ yields a precision of $8 \cdot 10^{-4}$ for a cost of 304 calls to the equations, while $n = 250$ is required to reach that precision with the FS, for a cost of 1500 calls. For a similar computational budget, the PL reaches a precision $\sim 2 \cdot 10^{-5}$ with $n = 37$.

## 4 Conclusion

In this paper, we jointly analyzed the behavior of ladder methods to compute parallel transport. We first gave a Taylor expansion of one step of Schild's ladder. Then, we showed that when scaling the vector to transport by $\frac{1}{n^2}$, a quadratic speed of convergence is reached. Our numerical experiments illustrate that this bound is indeed reached. In the same framework, we bridged the gap between the Fanning scheme and Schild's ladder, shedding light on how SL could be turned into a second-order method, while the FS cannot.

For manifolds with no closed-form geodesics, we introduced the infinitesimal ladder schemes and showed that the PL converges with the same order as its counterpart with exact exp and log maps. The same exercise can be realized with SL. Numerical

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>
experiments were performed on SE(3) with anisotropic metric and allowed to observe the role of $\nabla R$ in the convergence of the PL in a non-symmetric space. This result corroborates our theoretical developments and shows that the bounds on the speed of convergence are reached.

Our last comparison of SL and PL shows that, although more popular, SL is far less appealing that the PL as (i) it is more expensive to compute, (ii) it converges slower, (iii) it is less stable when using approximate geodesics, and (iv) it is not exact in symmetric spaces. Pole ladder is also more appealing than the FS because of its quadratic speed of convergence, which allows to reach mild convergence tolerance at a much lower overall cost despite the higher complexity of each step.

The ladder methods are restricted to transporting along geodesics, but this not a major drawback as this is common to other methods, and one may approximate any curve by a piecewise geodesic curve.

In our work on infinitesimal schemes, we only used basic integration of ODEs in charts, while there is a wide literature on numerical methods on manifolds. In future work, it would be very interesting to consider specifically adapted RK schemes on Lie groups and homogeneous spaces [26,27], or symplectic integrators [4,14] in order to reduce the computational burden and to improve the stability of the scheme. More precisely, the computations of the PL may be hindered by the approximation of the log, while in fact, only a geodesic symmetry $y \mapsto \exp_x(-\log_x(y))$ is necessary and may be computed more accurately with a symmetric scheme.

**Acknowledgements** The authors have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement G-Statistics No 786854). It was also supported by the French government through the 3IA Côte d'Azur Investments ANR-19-P3IA-0002 managed by the National Research Agency. The authors warmly thank Yann Thanwerdas and Paul Balondrade for insightful discussions and careful proofreading of this manuscript.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the 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://creativecommons.org/licenses/by/4.0/.

## A Computation of the Expansion of Schild's Ladder

The details of the Taylor approximation for SL are given below at the fourth order, and a lemma to bound the fourth-order terms is proved. First we combine (1) and (2) to compute an approximation of $a = \log_x(m)$ where $m$ is the midpoint between $x_v$ and $x_w$. That is, $a = h_x(v, \frac{1}{2}b)$ where $b = l_x(v, w)$:

$$ 
\begin{aligned}
2a &= 2v + b + \frac{1}{6}R(b, v)(v + b) + \frac{1}{24}((\nabla_v R)(b, v)(\frac{5}{2}b + 2v) + (\nabla_b R)(\frac{1}{2}b, v)(v + b)) + O(5) \\
&= 2v + w - v + \frac{1}{6}R(v, w)(2w - v) + \frac{1}{24}((\nabla_v R)(v, w)(3w - 2v)
\end{aligned}
$$

<figure><img src="image_40.png" /></figure>
$$
\begin{align*}
& + (\nabla_w R)(v, w)(2w - v) + \frac{1}{6} R(\frac{1}{2}(w - v), v)(v + w - v) \\
& + \frac{1}{24} \left( (\nabla_v R)(w - v, v)\left(\frac{5}{2}(w - v) + 2v\right) + (\nabla_{w-v} R)\left(\frac{1}{2}(w - v), v\right)(v + w - v) \right) + O(5) \\
&= w + v + \frac{1}{6} R(v, w)(2w - v - w) + \frac{1}{24} \left( (\nabla_v R)(v, w)(3w - 2v - 2w - \frac{5}{2}(w - v)) \right. \\
& \qquad \left. + (\nabla_w R)(v, w)(2w - v) + (\nabla_{w-v} R)(v, w)(-\frac{1}{2}w) \right) + O(5).
\end{align*}
$$

Therefore,

$$
\begin{equation}
\begin{split}
2a ={}& w + v + \frac{1}{6} R(v, w)(w - v) \\
      & + \frac{1}{24} \left( (\nabla_v R)(v, w)(w - \frac{3}{2}v) + (\nabla_w R)(w, v)(v - \frac{3}{2}w) \right) + O(5).
\end{split}
\tag{31}
\end{equation}
$$

Now we compute $u = l_x(w, 2a)$:

$$
\begin{align*}
u &= 2a - w + \frac{1}{6} R(2a, w)(w - 4a) \\
  &\quad + \frac{1}{24} \left( (\nabla_w R)(w, 2a)(3 * 2a - 2w) + (\nabla_{2a} R)(w, 2a)(2 * 2a - w) \right) + O(5) \\
  &= w + v + \frac{1}{6} R(v, w)(w - v) \\
  &\quad + \frac{1}{24} \left( (\nabla_v R)(v, w)(w - \frac{3}{2}v) + (\nabla_w R)(w, v)(v - \frac{3}{2}w) \right) \\
  &\quad - w + \frac{1}{6} R(w + v, w)(w - 2(w + v)) + \frac{1}{24} \left( (\nabla_w R)(w, w + v)(3w + 3v - 2w) \right. \\
  &\qquad \left. + (\nabla_{w+v} R)(w, w + v)(2w + 2v - w) \right) + O(5) \\
  &= v + \frac{1}{6} R(v, w)(w - v - w - 2v) \\
  &\quad + \frac{1}{24} \left( (\nabla_v R)(v, w)(w - \frac{3}{2}v - 2v - w) + (\nabla_w R)(v, w)(v - \frac{3}{2}w - w - 3v - 2v - w) \right).
\end{align*}
$$

Thus

$$
\begin{equation}
\begin{split}
u = v & - \frac{1}{2} R(v, w)v \\
      & + \frac{1}{24} \left( (\nabla_v R)(v, w)\left(-\frac{7}{2}v\right) + (\nabla_w R)(v, w)\left(-4v - \frac{7}{2}w\right) \right) + O(5).
\end{split}
\tag{32}
\end{equation}
$$

We deduce the following

**Lemma 4** With the previous notations, at $x$ in a compact set $K$, $u$ can be written
$u = v + \frac{1}{2}R(w, v)v + r_4(v, w)$ and there exists $C > 0$ such that $\exists \eta > 0, \forall v, w \in$
$B_0(\delta) \subset \bar{T}_x M, \forall |s|, |t| \le \eta, r_4$ verifies:

$$
\|r_4(sv, tw)\| \le C(s^3 t + s^2 t^2 + st^3)\delta^4. \quad (33)
$$

If furthermore $|s| \le |t|$, then this reduces to $\|r_4(sv, tw)\| \le Cst^3\delta^4$. Moreover, $C$ can be bounded by bounds on the covariant derivatives of the curvature tensor.

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>
**Proof** Let $t' = \min(|t|, |s|)$. By (32) we get

$$r_4(sv, tw) = \frac{st}{24} \left( (\nabla_v R)(w, sv)(\frac{7}{2}sv) + (\nabla_w R)(tw, v)(4sv + \frac{7t}{2}w) \right) + O(5).$$

Each term of the form $\nabla.R(tw, sv)$ can be bounded, for example:

$$\|st(\nabla_v R)(w, sv)(\frac{7}{2}sv)\| \leq s^3 t \|\nabla R\|_{\infty} \|w\| \|v\|^3, \quad (34)$$

where the infimum norm on $\nabla R$ is taken on the compact set $K$ (thus it exists and it is finite). Similarly, as $\|w\| \leq \delta$ and $\|v\| \leq \delta$

$$\|s(\nabla_v R)(tw, sv)(\frac{7}{2}sv)\| \leq s^3 t \|\nabla R\|_{\infty} \delta^3. \quad (35)$$

By dealing in a similar fashion with the two other terms, we obtain:

$$\|r_4(sv, tw)\| \leq \|\nabla R\|_{\infty}(s^3t + s^2t^2 + st^3)\delta^4 + O(5), \quad (36)$$

where $O(5)$ is a combination of terms of the form $q_5(s\alpha, t\beta)$, which are homogeneous polynomials of degree at least five and variables in the ball of radius $\delta$; hence, they can be bounded above by $C_1(s^3t + s^2t^2 + st^3)\delta^4$ for some $C_1$. The result follows with $C = \|\nabla R\|_{\infty} + C_1$. □

## B Computation of the Expansion of the Pole Ladder

As in the previous appendix, we give here the computations for the fourth-order Taylor approximation of the PL construction. Recall (Fig. 8) that $v, w \in T_m M$ where $m = \gamma(\frac{1}{2})$ is the midpoint between $x = \gamma(0)$ and $\gamma(1)$. The result of the construction parallel transported back to $m$ is $u$ such that:

$$\begin{aligned} -u &= l_m(w, -h_m(-w, v)) \\ &= w + h_m(-w, v) + \frac{1}{6}R(w, -h_m(-w, v))(-2h_m(-w, v) - w) \\ &\quad + \frac{1}{24}((\nabla_w R)(w, -h_m(-w, v))(-3h_m(-w, v) - 2w) \\ &\quad + (\nabla_{-h_m(-w, v)} R)(w, -h_m(-w, v))(-2h_m(-w, v) - w)) + O(5). \end{aligned}$$

Now we plug (1), but it reduces to $-h_m(-w, v) = w - v + O(3)$ when it appears in a term including curvature (because we restrict to fourth-order terms overall). Hence,

$$\begin{aligned} -u &= -(-w + v + \frac{1}{6}R(v, -w)(-w + 2v) + \frac{1}{24}((\nabla_{-w} R)(v, -w)(5v - 2w) \\ &\quad + (\nabla_v R)(v, -w)(-w + 2v)) + q_5(v, v)) + \frac{1}{6}R(w, w - v)(2(w - v) - w) \end{aligned}$$

<figure><img src="image_43.png" /></figure>
$$
\begin{align*}
& + \frac{1}{24}((\nabla_w R)(w, w-v)(3(w-v)-2w) + (\nabla_{w-v} R)(w, w-v)(2(w-v)-w)) + O(5) \\
&= -v + \frac{1}{6}R(v, w)(-w+2v) + \frac{1}{6}R(w, -v)(w-2v) \\
& + \frac{1}{24}((\nabla_w R)(v, -w)(5v-2w) + (\nabla_v R)(v, w)(-w+2v)) \\
& + \frac{1}{24}((\nabla_w R)(v, w)(w-3v) + (\nabla_{w-v} R)(v, w)(w-2v)) + O(5) \\
&= -v + \frac{1}{12}((\nabla_w R)(v, w)(5v-w) + (\nabla_v R)(v, w)(2v-w)) + O(5).
\end{align*}
$$

C Geometry of the Sphere

The unit sphere is defined as $S^2 = \{x \in \mathbb{R}^3, \|x\| = 1\}$. With the canonical scalar product $\langle \cdot, \cdot \rangle$ of $\mathbb{R}^3$, it is a Riemannian manifold of constant curvature whose geodesics are great circles. The tangent space of any $x \in S^2$ is the set of vectors orthogonal to $x$: $T_x S^2 = \{v \in \mathbb{R}^3, \langle x, v \rangle = 0\}$. We have $\forall x, y \in S^2, \forall v, w \in T_x S^2$:

$$
\begin{align*}
\exp_x(w) &= \cos(\|w\|)x + \sin(\|w\|) \frac{w}{\|w\|}, \\
\log_x(y) &= \arccos(\langle y, x \rangle) \frac{y - \langle y, x \rangle x}{\|y - \langle y, x \rangle x\|}, \\
\Pi_x^{xw} v &= \langle v, \frac{w}{\|w\|} \rangle \left( -\sin(\|w\|)x + \cos(\|w\|) \frac{w}{\|w\|} \right) + \left( v - \langle v, \frac{w}{\|w\|} \rangle \frac{w}{\|w\|} \right).
\end{align*}
$$

The expression for parallel transport of $v$ means that the orthogonal projection of $v$ on $\{w\}^{\perp}$ is preserved, while the orthogonal projection on $\mathbb{R}w$ is rotated by an angle $\|w\|$ in the $(x, w)$-plane.

D Affine-Invariant Geometry of SPD Matrices

The cone of 3 × 3 symmetric positive-definite matrices is defined as

$$
SPD(3) = \{ \Sigma \in \mathbb{R}^{3 \times 3}, \Sigma^T = \Sigma, \forall x \in \mathbb{R}^3 \setminus \{0\}, x^T \Sigma x > 0 \}.
$$

The tangent space of any $\Sigma \in SPD(3)$ is the set symmetric matrices: $T_\Sigma SPD(3) = \text{Sym}(3)$. The affine-invariant (AI) metric is defined at any $\Sigma \in SPD(3)$, for all $V, W \in \text{Sym}(3)$ using the matrix trace $\text{tr}$ by

$$
g_{\Sigma}(V, W) = \text{tr}(\Sigma^{-1}V\Sigma^{-1}W).
$$

We have $\forall \Sigma, \Sigma_1, \Sigma_2 \in SPD(3), \forall W \in \text{Sym}(3)$

$$
\begin{align*}
\exp_{\Sigma}(W) &= \Sigma^{\frac{1}{2}} \exp(\Sigma^{-\frac{1}{2}} W \Sigma^{-\frac{1}{2}}) \Sigma^{\frac{1}{2}}, \\
\log_{\Sigma_1}(\Sigma_2) &= \Sigma_1^{\frac{1}{2}} \log(\Sigma_1^{-\frac{1}{2}} \Sigma_2 \Sigma_1^{-\frac{1}{2}}) \Sigma_1^{\frac{1}{2}},
\end{align*}
$$

<figure><img src="image_44.png" /></figure>
where when not indexed, exp and log refer to the matrix operators. Finally, let

$$P_t = \Sigma^{\frac{1}{2}} \exp(\frac{t}{2} \Sigma^{-\frac{1}{2}} W \Sigma^{-\frac{1}{2}}) \Sigma^{-\frac{1}{2}}$$

The parallel transport from $\Sigma$ along the geodesic with initial velocity $W \in Sym(3)$ of $V \in Sym(3)$ a time $t$ is ([33])

$$\Pi_{0,W}^t V = P_t V P_t^T.$$

## E Left-Invariant Metric on SE(3)

In this appendix, we describe the geometry of the Lie group of isometries of $\mathbb{R}^3$, $SE(3)$, endowed with a left-invariant metric $g$. We prove Lemma 3 that gives a necessary and sufficient condition for $(SE(3), g)$ to be a Riemannian locally symmetric space. We also give the details for the computation of the geodesics with numerical integration schemes. Those results and computations are valid in any dimension $d \ge 2$, but we only detail them for $d = 3$ to keep them tractable. For the details on this section, we refer the reader to [7,18,23].

$SE(3)$ is the semi-direct product of the group of three-dimensional rotations $SO(3)$ with $\mathbb{R}^3$, i.e., the group multiplicative law for $R, R' \in SO(3), t, t' \in \mathbb{R}^3$ is given by

$$(R, t) \cdot (R', t') = (RR', t + Rt').$$

It can be seen as a subgroup of $GL(4)$ and represented by homogeneous coordinates:

$$(R, t) = \begin{pmatrix} R & t \\ 0 & 1 \end{pmatrix},$$

and all group operations then correspond to the matrix operations. Let the metric matrix at the identity be diagonal: $G = \text{diag}(1, 1, 1, \beta, 1, 1)$ for some $\beta > 0$, the anisotropy parameter. We write $\langle \cdot, \cdot \rangle$ for the associated inner-product at the identity. An orthonormal basis of the Lie algebra $se(3)$ is

$$e_1 = \frac{1}{\sqrt{2}} \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$$

$$e_2 = \frac{1}{\sqrt{2}} \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$$

$$e_3 = \frac{1}{\sqrt{2}} \begin{pmatrix} 0 & -1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$$

$$e_4 = \frac{1}{\sqrt{\beta}} \begin{pmatrix} 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$$

$$e_5 = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}$$

$$e_6 = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{pmatrix}.$$

<figure><img src="image_45.png" /></figure>
Define the corresponding structure constants $C_{ij}^k =< [e_i, e_j], e_k >$, where the Lie bracket $[\cdot, \cdot]$ is the usual matrix commutator. It is straightforward to compute

$$C_{ij}^k = \frac{1}{\sqrt{2}} \text{ if } ijk \text{ is a direct cycle of } \{1, 2, 3\}; \qquad (37)$$

$$C_{15}^6 = -C_{16}^5 = -\sqrt{\beta}C_{24}^6 = \frac{1}{\sqrt{\beta}}C_{26}^4 = \sqrt{\beta}C_{34}^5 = -\frac{1}{\sqrt{\beta}}C_{35}^4 = \frac{1}{\sqrt{2}}. \qquad (38)$$

and all others that cannot be deduced by skew symmetry of the bracket are equal to 0. Extend the inner product defined in the Lie algebra by $G$ to a left-invariant metric $g$ on $SE(3)$. Let $\nabla$ be its associated Levi-Civita connection. It is also left-invariant, and it is sufficient to know its values on left-invariant vector fields at the identity (identified with tangent vectors at the identity). These are linked to the structure constants by

$$\nabla_{e_i} e_j = \frac{1}{2} \sum_k (C_{ij}^k - C_{jk}^i + C_{ki}^j) e_k,$$

and can thus be computed explicitly thanks to (37). Let $\Gamma_{ij}^k =< \nabla_{e_i} e_j, e_k >$ be the associated Christoffel symbols. Let $\tau = (\sqrt{\beta} + \frac{1}{\sqrt{\beta}})$. We obtain

$$\Gamma_{ij}^k = \frac{1}{2\sqrt{2}} \text{ if } ijk \text{ is a cycle of } [1,2,3], \qquad (39)$$

$$\Gamma_{15}^6 = -\Gamma_{16}^5 = -\frac{2}{\tau}\Gamma_{24}^6 = \frac{2}{\tau}\Gamma_{26}^4 = \frac{2}{\tau}\Gamma_{34}^5 = -\frac{2}{\tau}\Gamma_{35}^4 = \frac{1}{\sqrt{2}}, \qquad (40)$$

and all other are null. Finally, recall that the curvature tensor and its covariant derivative at the identity can be defined $\forall u, v, w, z \in \mathfrak{se}(3)$ by

$$R(u, v) = \nabla_u \nabla_v - \nabla_v \nabla_u - \nabla_{[u, v]} \\ (\nabla_u R)(v, w)_z = \nabla_u(R(v, w)_z) - R(\nabla_u v, w)_z - R(v, \nabla_u w)_z - R(v, w)\nabla_u z.$$

# Geodesics

## General Case

For the computation of the geodesics, let $(x, v) \in TSE(3)$ and $\gamma$ be the geodesic such that $\gamma(0) = x, \dot{\gamma}(0) = v$. At all times $t$, define the Eulerian velocity $X(t) = (L_{\gamma^{-1}(t)})_* \dot{\gamma}(t) \in \mathfrak{se}(3)$ where $*$ refers to the push-forward and $L$ to the left translation $L_g : h \mapsto gh$ of $SE(3)$. Define the co-adjoint action $ad^*$ associated with the metric:

$$\forall a, b, c \in \mathfrak{se}(3), \quad <[a, b], c> = <ad_a^* c, b>. \qquad (41)$$

<figure><img src="image_46.png" /></figure>
It can be computed explicitly using the structure constants. Then, we have the evolution equations:

$$ \dot{\gamma}(t) = (L_{\gamma(t)})_* X(t), \qquad (42) $$

$$ \dot{X}(t) = ad_{X(t)}^* X(t), \qquad (43) $$

with initial conditions $\gamma(0) = (R_0, v_0)$ and $\dot{\gamma}(0) = (Q_0, u_0)$, so that $X(0) = (L_{x^{-1}})_*(Q, u_0) = (R_0^T Q, R_0^T u_0) = (\Omega_0, d_0)$.

### Case $\beta = 1$

Focusing on the case $\beta = 1$ allows to validate our implementation by testing that we recover the direct product exponential map. This is proved in, e.g., [35]. It is straightforward to see that, in this case, $g$ coincides with the (direct) product of the canonical (bi-invariant) metrics $g_{rot}$, $g_{trans}$ of $SO(d)$ and $\mathbb{R}^d$. So we have the general formula for geodesics from $(R, v) \in SE(3)$ with initial velocity $(\Omega, u) \in T_{(R,v)}SE(3): \gamma(t) = (R \exp(tR^T\Omega), tu + v)$, and

$$ \exp_{(R,v)}(\Omega, u) = (R \exp(R^T \Omega), u + v), \qquad (44) $$

$$ \log_{(R,v)}(Q, u) = (R \log(R^T Q), u - v). \qquad (45) $$

These are in fact valid in $SE(d)$ for any $d \ge 2$. These geodesics are represented in Fig. 14 for different values of $\beta$ in $SE(2)$ (where the metric matrix at identity is $G = \text{diag}(1, \beta, 1)$).

<figure><img src="image_47.png" /><figcaption>Fig. 14 Visualization of geodesics of $SE(2)$ with the same initial conditions for different values of the anisotropy parameter $\beta$. The translation part are black dots, and the rotation part is applied to the orthonormal canonical frame of the 2d-plane</figcaption></figure>

<figure><img src="image_48.png" /></figure>
## Proof of Lemma 3

We now prove Lemma 3, formulated as: $(SE(3), g)$ is locally symmetric, i.e., $\nabla R = 0$, if and only if $\beta = 1$. This is valid for any dimension $d \ge 2$ provided that the metric matrix $G$ is diagonal, of size $d(d+1)/2$, with ones everywhere except one coefficient of the translation part.

**Proof** For $\beta = 1$, $(SE(d), g)$ is isometric to $(SO(d) \times \mathbb{R}^d, g_{rot} \oplus g_{trans})$. As the product of two symmetric spaces is again symmetric, $(SE(d), g)$ is symmetric.

We prove the contraposition of the necessary condition. Let $\beta \ne 1$. We give $i, j, k, l$ s.t. $(\nabla_{e_i} R)(e_j, e_k)e_l \ne 0$:

$$ (\nabla_{e_3} R)(e_3, e_2)e_4 = \nabla_{e_3}(R(e_3, e_2)e_4) - R(e_3, \nabla_{e_3} e_2)e_4 - R(e_3, e_2)\nabla_{e_3} e_4 \\ = \nabla_{e_3}(R(e_3, e_2)e_4) + \frac{1}{\sqrt{2}}R(e_3, e_1)e_4 - \frac{\tau}{2\sqrt{2}}R(e_3, e_2)e_5. $$

And from the above

$$ \begin{aligned} R(e_3, e_2)e_4 &= \nabla_{e_3}\nabla_{e_2}e_4 - \nabla_{e_2}\nabla_{e_3}e_4 - \nabla_{[e_3, e_2]}e_4 \\ &= -\frac{\tau}{2\sqrt{2}}\nabla_{e_3}e_6 - \nabla_{e_2}e_5 + \frac{1}{\sqrt{2}}\nabla_{e_1}e_4 \\ &= 0. \end{aligned} $$

$$ \begin{aligned} R(e_3, e_1)e_4 &= \nabla_{e_3}\nabla_{e_1}e_4 - \nabla_{e_1}\nabla_{e_3}e_4 - \nabla_{[e_3, e_1]}e_4 \\ &= -\frac{\tau}{2\sqrt{2}}\nabla_{e_1}e_5 - \frac{1}{\sqrt{2}}\nabla_{e_2}e_4 \\ &= -\frac{\tau}{4}e_6 + \frac{\tau}{4}e_6 = 0. \end{aligned} $$

$$ \begin{aligned} R(e_3, e_2)e_5 &= \nabla_{e_3}\nabla_{e_2}e_5 - \nabla_{e_2}\nabla_{e_3}e_5 - \nabla_{[e_3, e_2]}e_5 \\ &= \frac{\tau}{2\sqrt{2}}\nabla_{e_2}e_4 + \frac{1}{\sqrt{2}}\nabla_{e_1}e_5 \\ &= -\frac{\tau^2}{8}e_6 + \frac{1}{2}e_6 = \frac{1}{2}(1 - \frac{\tau^2}{4})e_6. \end{aligned} $$

And therefore

$$ \beta \ne 1 \implies \tau = (\sqrt{\beta} + \frac{1}{\sqrt{\beta}}) \ne 2 \implies (\nabla_{e_3} R)(e_3, e_1)e_4 = -\frac{\tau}{4\sqrt{2}}(1 - \frac{\tau^2}{4})e_6 \ne 0, $$

which proves Lemma 3. □

<figure><img src="image_49.png" /></figure>
## F Implementation

Our implementation of the ladder methods and the Fanning scheme are made available online at https://github.com/nguigs/ladder-methods. The repository also contains a notebook to reproduce all the experiments of this paper. It relies on the open-source Python package `geomstats`, with its default `numpy` back-end. The automatic differentiation of the `autograd` package is used to compute the logs of the infinitesimal schemes, with `scipy's` L-BFGS-B solver for the gradient descent.

## References

1. B. Berkels, A. Effland, M. Rumpf, Time discrete geodesic paths in the space of images. SIAM J. Imaging Science (SIIMS) 8(3), 1457-1488 (2015). https://doi.org/10.1137/140970719

2. D. Brooks, O. Schwander, F. Barbaresco, J.Y. Schneider, M. Cord, Riemannian batch normalization for SPD neural networks. In: H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché Buc, E. Fox, R. Garnett (eds.) Advances in Neural Information Processing Systems (NeurIPS) 32, pp. 15489-15500. Curran Associates Inc. (2019)

3. C. Cury, M. Lorenzi, D. Cash, et al., Spatio-temporal shape analysis of cross-sectional data for detection of early changes in neurodegenerative disease. In: SeSAMI 2016 - First International Workshop Spectral and Shape Analysis in Medical Imaging, LNCS, vol. 10126, pp. 63-75. Springer (2016). https://doi.org/10.1007/978-3-319-51237-2_6

4. J.P. Dedieu, D. Nowicki, Symplectic methods for the approximation of the exponential map and the Newton iteration on Riemannian submanifolds. Journal of Complexity 21(4), 487-501 (2005). https://doi.org/10.1016/j.jco.2004.09.010

5. J. Ehlers, F.A.E. Pirani, A. Schild, The geometry of free fall and light propagation. In: L. O'Raifeartaigh (ed.) General relativity: papers in honour of J. L. Synge, pp. 63-84. Oxford: Clarendon Press (1972)

6. O. Freifeld, S. Hauberg, M.J. Black, Model transport: towards scalable transfer learning on manifolds. In: Proceedings IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), pp. 1378-1385 (2014)

7. J. Gallier, J. Quaintance, Differential geometry and Lie groups: a computational perspective. Geometry and Computing. Springer International Publishing (2020). https://doi.org/10.1007/978-3-030-46040-2

8. A.V. Gavrilov, Algebraic properties of covariant derivative and composition of exponential maps. Siberian Adv. Math. 16(3), 54-70 (2006)

9. A.V. Gavrilov, The double exponential map and covariant derivation. Sib Math J 48(1), 56-61 (2007). https://doi.org/10.1007/s11202-007-0006-4

10. A.V. Gavrilov, The affine connection in the normal coordinates. Sib. Adv. Math. 23(1), 1-19 (2013). https://doi.org/10.3103/S105513441301001X

11. N. Guigui, S. Jia, M. Sermesant, X. Pennec, Symmetric algorithmic components for shape analysis with diffeomorphisms. In: Proc. of GSI 2019 - 4th conference on Geometric Science of Information, LNCS, vol. 11712, pp. 759-768. Springer (2019). https://doi.org/10.1007/978-3-030-26980-7_79

12. N. Guigui, E. Maignant, A. Trouvé, X. Pennec, Parallel transport on Kendall shape spaces (2021). ArXiv:2103.04611

13. N. Guigui, X. Pennec, A reduced parallel transport equation on Lie groups with a left-invariant metric (2021). ArXiv:2103.01585

14. E. Hairer, G. Wanner, C. Lubich, Symplectic Integration of Hamiltonian Systems, pp. 167-208. Springer Series in Computational Mathematics. Springer, Berlin, Heidelberg (2002). https://doi.org/10.1007/978-3-662-05018-7_6

15. S. Hauberg, F. Lauze, K.S. Pedersen, Unscented Kalman filtering on Riemannian manifolds. J. Math Imaging Vis 46(1), 103-120 (2013). https://doi.org/10.1007/s10851-012-0372-9

16. A. Kheyfets, W.A. Miller, G.A. Newton, Schild's ladder parallel transport procedure for an arbitrary connection. Int. J. of Theoretical Physics 39(12), 2891-2898 (2000). https://doi.org/10.1023/A:1026473418439
17. K.R. Kim, I.L. Dryden, H. Le, K.E. Severn, Smoothing splines on Riemannian manifolds, with applications to 3D shape space. Journal of the Royal Statistical Society: Series B (Statistical Methodology) **83**(1), 108–132 (2020). https://doi.org/10.1111/rssb.12402

18. B. Kolev, Lie groups and mechanics: an introduction. Journal of Nonlinear Mathematical Physics **11**(4), 480–498 (2004). https://doi.org/10.2991/jnmp.2004.11.4.5. ArXiv:math-ph/0402052

19. M. Lorenzi, N. Ayache, X. Pennec, Schild's ladder for the parallel transport of deformations in time series of images. In: G. Székely, H.K. Hahn (eds.) Information Processing in Medical Imaging, LNCS, vol. 6801, pp. 463–474. Springer, Berlin, Heidelberg (2011). https://doi.org/10.1007/978-3-642-22092-0_38

20. M. Lorenzi, X. Pennec, Efficient parallel transport of deformations in time series of images: from Schild to pole ladder. J Math Imaging Vis **50**(1), 5–17 (2014). https://doi.org/10.1007/s10851-013-0470-3

21. M. Louis, Computational and statistical methods for trajectory analysis in a Riemannian geometry setting. Ph.D. thesis, Sorbonne université (2019)

22. M. Louis, B. Charlier, P. Jusselin, S. Pal, S. Durrleman, A fanning scheme for the parallel transport along geodesics on Riemannian manifolds. SIAM Journal on Numerical Analysis **56**(4), 2563–2584 (2018). https://doi.org/10.1137/17M1130617

23. J. Milnor, Curvatures of left invariant metrics on Lie groups. Advances in Mathematics **21**(3), 293–329 (1976). https://doi.org/10.1016/S0001-8708(76)80002-3

24. N. Miolane, N. Guigui, A.L. Brigant, et al., Geomstats: a Python package for Riemannian geometry in machine learning. Journal of Machine Learning Research **21**(223), 1–9 (2020)

25. C.W. Misner, K.S. Thorne, J.A. Wheeler, Gravitation. Princeton University Press (1973)

26. H. Munthe-Kaas, O. Verdier, Integrators on homogeneous spaces: isotropy choice and connections. Found Comput Math **16**(4), 899–939 (2016). https://doi.org/10.1007/s10208-015-9267-7

27. H. Munthe-Kaas, A. Zanna, Numerical integration of differential equations on homogeneous manifolds. In: F. Cucker, M. Shub (eds.) Foundations of Computational Mathematics, pp. 305–315. Springer, Berlin, Heidelberg (1997). https://doi.org/10.1007/978-3-642-60539-0_24

28. X. Pennec, Parallel transport with pole ladder: a third order scheme in affine connection spaces which is exact in affine symmetric spaces (2018). ArXiv:1805.11436

29. X. Pennec, Curvature effects on the empirical mean in Riemannian and affine Manifolds: a non-asymptotic high concentration expansion in the small-sample regime (2019). ArXiv:1906.07418

30. X. Pennec, Manifold-valued image processing with SPD matrices. In: X. Pennec, S. Sommer, T. Fletcher (eds.) Riemannian geometric statistics in medical image analysis, pp. 75–134. Academic Press (2020). https://doi.org/10.1016/B978-0-12-814725-2.00010-8

31. X. Pennec, S. Sommer, T. Fletcher (eds.), Riemannian geometric statistics in medical image analysis, *The Elsevier and MICCAI Society book series*, vol. 3. Elsevier (2020). https://doi.org/10.1016/C2017-0-01561-6

32. J.B. Schiratti, S. Allassonnière, O. Colliot, S. Durrleman, A Bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations. J. Machine Learning Research (JMLR) **18**(133), 33 (2017)

33. O. Yair, M. Ben-Chen, R. Talmon, Parallel transport on the cone manifold of SPD matrices for domain adaptation. In: IEEE Transactions on Signal Processing, vol. 67, pp. 1797–1811 (2019). https://doi.org/10.1109/TSP.2019.2894801

34. L. Younes, Jacobi fields in groups of diffeomorphisms and applications. Quarterly of Applied Mathematics **65**(1), 113–134 (2007). https://doi.org/10.1090/S0033-569X-07-01027-5

35. M. Zefran, V. Kumar, C. Croke, On the generation of smooth three-dimensional rigid body motions. IEEE Transactions on Robotics and Automation **14**(4), 576–589 (1998). https://doi.org/10.1109/70.704225

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>