Numerical Accuracy of Ladder Schemes for Parallel Transport on Manifolds

Parallel transport is a fundamental tool to perform statistics on Rie-mannian manifolds. Since closed formulae don't 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.

to perform statistics and machine learning on manifolds [27]. 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, [1] use it on the manifold of symmetric positive definite (SPD) matrices to centralize batches when training a neural network. In [29], parallel transport is used again on SPD matrices for domain adaptation. In [13], it is a key ingredient to spline-fitting in a Kendall shape space. In computational anatomy, it allows to compare longitudinal, intrasubject evolution across populations [2,15,16,28]. It is also used in computer vision in e.g. [5,11].
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 [18,30]. For instance, [13] derived a homogeneous first-order differential equation stemming from the structure of quotient space that defines Kendall shape spaces. [18] 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 [15,16], 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 1 . Its first appearance in the literature is in [21] 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 [12]. 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 [12], 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 [7,8], 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 1 n α can be obtained where 1 ≤ α ≤ 2 and n is the number of parallelograms, or rungs of the ladder. This improves over the assumption made in [18] that only a first-order convergence speed can be reached with SL. These results are observed with a high accuracy in numerical experiments performed using the open-source Python package geomstats 2 [20].
A slight modification of this scheme, Pole Ladder (PL), was introduced in [16]. 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 [24]. 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) [18], 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 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 [31]. 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 [19] 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 the experiments of the paper are available online at 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. Fig. 2 The double exponential and its inverse (left) and the neighboring logarithm (right) represented in a normal coordinate system centred at x. Geodesics are in black and tangent vectors are coloured

Notations and Assumptions
We consider a complete Riemannian manifold (M, g) of finite dimension d ∈ N. The associated Levi-Civita connection defines a covariant derivative ∇ and the parallel transport map. Denote exp the Riemannian exponential map, and log its inverse, defined locally. For x ∈ M, let T x M be the tangent space at x. The map exp x sends T x M to (a subset of) M, and we will often write of v along γ is defined as the unique solution at time t ≤ 1 to the ODE ∇γ (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 γ is a geodesic starting at x, and let w be its initial velocity, i.e. γ(t) = exp x (tw). Thus the dependence on γ in the notation Π γ will be omitted, and we instead write Π y x 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 · the norm defined on each tangent space by the metric g and R the Riemann curvature tensor that maps for any x ∈ M, u, v, w ∈ T x M to R(u, v)w ∈ T x M. Throughout this paper, we consider γ to be contained in a compact set K ⊂ M of diameter δ > 0, and thus R and all its covariant derivatives ∇ n R can be uniformly bounded.

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) [7,8]. It is defined for x ∈ M, (v, w) ∈ (T x 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 : 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 : As it is a smooth map, we can apply the Taylor theorem to (s, t) → h x (sv, tw). As explained in [7,8], its derivative are invariants of the connection that can be expressed in terms of the curvature and its covariant derivatives.
In this work we will require an approximation at order four, for s, t ∈ R small enough: 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 ∼ 1 n for some n ∈ N large enough, q 5 (sv, tw) = O( 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, [25] introduced the neighboring log and computed its Taylor approximation. It is defined by applying the exp x map to small enough v, w ∈ T x M to obtain x x , 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 : which relates to the double exponential by solving Using (1), one can solve for the first terms of a Taylor expansion and obtains for (v, w) ∈ U x ⊂ T x M × T x M [25]: 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.

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

Elementary Construction
The construction to parallel transport v ∈ T x M along the geodesic γ with γ(0) = x andγ(0) = w ∈ T x M (such that (v, w) ∈ 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.
This is the first diagonal of the parallelogram. 3. Compute the geodesic between x and m, let a ∈ 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.
This is the second diagonal of the parallelogram. 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.
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 γ without further precision.

Taylor expansion
We can now reformulate this elementary construction in terms of successive applications of the double exponential and the neighboring logarithm maps.
be the initial velocity of the geodesic computed in step 2 of the above construction, and b its parallel transport to x. Therefore, Thus step 3. is equivalent to a = h x (v, 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 ∇ · R(·, ·)· are summarized in the term O(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 = Π x w x v −u w , we parallel transport it back to x: . 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 ∈ M and v, w ∈ T x M sufficiently small. Then the output u of one step of Schild's ladder parallel transported back to x is given by 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 transport. Furthermore, this shows that splitting the main geodesic into segments of length 1 n and simply iterating this construction n times will in fact sum n error terms of the form R( w i n , v i )v i , hence by linearity of R, the error won't necessarily vanish as n → ∞. 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.

Numerical Scheme and proof of convergence
With the notations of the previous subsection, let us define schild(x, w, v) = u w ∈ T x w M . We now divide the geodesic γ into segments of length 1 n for some n ∈ N * large enough so that the previous Taylor expansions can be applied to w n and v n . As mentioned before, v needs to be scaled as w. In fact let α ≥ 1 and consider the sequence defined by (see Fig. 4) where We now establish the following result, which ensures convergence of Schild's ladder to the parallel transport of v along γ at order at most two.
Moreover, τ is bounded by a bound on the sectional curvature, and β 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 rungs, parallel transported to x (isometrically): By (5) and lemma 4, the error at each step can be written where r 4 contains the fourth order residual terms, and is given in appendix A. We start as in [17] by assuming that v i doesn't grow too much, i.e. there exists s ∈ {1, . . . , n} such that ∀i < s, v i ≤ 2 v . We will then show that the control obtained on v n is tight enough so that when n is large enough, ∀k ≤ n, v k ≤ 2 v . With this assumption each term of the right-hand side (r.h.s.) can be bounded. As α ≥ 1, we apply lemma 4 (in appendix A): ∃β > 0 such that for n large enough: Similarly, as ∀i ≤ n, w i = w and by assumption v i ≤ 2 v , we have As the r.h.s. does not depend on n, it can be plugged into (7) to obtain by summing for i = 0, . . . , s − 1 ≤ n: Now suppose for contradiction that for arbitrarily large n, there exists k ≤ n such that v k > 2 v and choose u(n) ∈ 1, . . . , 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: But the r.h.s. goes to 0 as n goes to infinity, leading to a contradiction. Therefore, for n large enough, ∀i ≤ n, v i ≤ 2 v , and the previous control on v − Π x x s v s given by (10) is valid for s = n. The result follows by parallel transporting this error to x n , which doesn't change its norm.

Remark 1
The proof allows to grasp the origin of the bounds in Thm. 1 (5), as it is bilinear in v, so that the division by n α is squared in the error, while we only need to multiply by n α 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 doesn't appear in the error.
On the other hand, the O(n −2 ) comes from ∇ w R(w, v)w (from lemma 4). Indeed, this term is invariant by the division/multiplication by n α , 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 α = 2 in practice.

Remark 2
The bound on the curvature that we use can be related to a bound on the sectional curvature κ as follows. Recall Suppose it is bounded on the compact set K: We then could have used τ = δ 3 κ ∞ 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 α ≥ 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.

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×3 symmetric positive-definite matrices SP D (3).
The Sphere We consider the sphere S 2 as the subset of unit vectors of R 3 , endowed with the canonical ambient metric < ·, · >. 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 κ is constant, the fourth-order terms vanish and theorem 2 becomes v − Π x x n v n ≤ τ n α for α ≥ 1. This is precisely observed on figure 5.
Moreover, using Eq. (8), we have at each rung: Therefore, by summing as in the proof of theorem 2 we obtain: 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 −α and it should equal 1 2 . This is very precisely what we observe (Fig.5, right).
SPD We now consider SP D(3) endowed with the affine-invariant metric [26]. 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 on Fig. 6.

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 non-zero ∇R.

The Fanning Scheme, revisited
The Fanning Scheme [18] leverages an identity given in [30] between Jacobi fields and parallel transport. As proved in [18], 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 ∈ T x M, define the Jacobi field (JF) along γ : t → exp x (tw) at time h small enough by: One can then show that J v γ(t) (h) = hΠ t+h γ,t v + O(h 2 ) [30]. The Fanning Scheme consists in computing perturbed geodesics, i.e. with initial velocity w + v for some > 0 and to approximate the associated JF by and the parallel transport of v along γ between 0 and h by Fig. 7 Elementary construction of the Fanning Scheme. For = 1 n , it is similar to that of SL (with α = 2) except that the midpoint is not computed but w + v is used instead of log x (m).
This defines an elementary construction that is iterated n times with time-steps h = 1 n , and [18] showed that using = h yields the desired convergence. Let's consider that the finite difference approximation of (14) is in fact a first-order approximation of the log of exp x h(w + v) from exp x (hw). Now this scheme is similar to SL, except that it uses 2a = w + v instead of (4), i.e. the correction term according to curvature is not accounted for (compare figures 3 and 7). When using h = = 1 n , this corresponds to α = 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 sec.
But as = h, the last term including v disappears in the O(h 4 ) and we obtain: 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 = 1 n , summing error terms yields a global error of order where Then it is straightforward to reproduce the proof of thm. 2 with (15) instead of (5) to show that This corroborates the result of [18], although at this points 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 doesn't 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.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: Hence 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 (∇ 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.

Pole Ladder
The pole ladder was introduced by [16] to reduce the computational burden of SL. They showed that at the first order, the constructions where equivalent.
[24] 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 [9]). We first present the scheme and derive the Taylor approximation using the neighboring log.

The elementary construction
The construction to parallel transport v ∈ T x M along the geodesic γ with γ(0) = x andγ(0) = w ∈ T x M (such that (v, w) ∈ 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 = 1 2 to obtain the midpoint m. The main geodesic γ is a diagonal of the parallelogram. 2. Compute the geodesic between x v and m, let a ∈ 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.
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.

Taylor approximation
We now introduce new notations, and centre 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 of the vectors v −w that we wish to transport between m −w and m w . With these new notations, . Therefore, using (1) and (2), we obtain the following (see appendix B for the computations) Theorem 3 Let (M, g) be a finite dimensional Riemannian manifold. Let m ∈ M and v, w ∈ T m M sufficiently small. Then the output u of one step of the pole ladder parallel transported back to m is given by 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(·, ·)·. Using n −1 to scale w and n −α to scale v as in SL, the term (∇ 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 α at the final rung of the ladder to recover the parallel transport of v. The other terms are multilinear in v, so if α > 1, there are small compared to 1 n 3 (∇ w R)(w, v)w. Summing n terms of the form 1 n 3 (∇ w R)(w, v)w we thus obtain a quadratic speed of convergence, as for Schild's ladder, for any α ≥ 1. We henceforth use α = 1.
As in Sec. 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 (∇ w R)(w, v)w. As in the previous section, define u w = pole(m, w, v) the result of the PL construction (with the new notations), and the sequence where Then it is straightforward to reproduce the proof of thm. 2 with (20) instead of (5) to show that ∃β > 0, ∃N ∈ N, ∀n > N , and β 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 Fig.9, bottom row. This greatly reduces the computational burden of the PL over SL.

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.
More precisely, consider the geodesic equation written as a first order equation of two variables in a global chart Φ of M , that defines for any p ∈ M a basis of where Γ k ij are the Christoffel symbols, and we use Einstein summation convention. Let rk : T M × R + → T M 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 = 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 Φ: and global accumulated error after n steps with step size h = 1 n : where by rk 1 we mean the projection on the first variable, and τ 1 > 0 is a constant. This means that any point on the geodesic is approximated with error O( 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: Furthermore, the inverse of v → rk 1 (x, v, h) can be computed for any x by gradient descent, and we assume 3 that any desired accuracy can be reached. We writeṽ = rk −1 x (y) for the optimal v such that d(rk(x, v, 1 n ), y) = o( 1 n 5 ), and assume that rk −1 x (y)−log x (y) ≤ τ 3 n 4 . Note that the step size does not appear explicitly in the notation rk −1 , but h = 1 n is always used. As in practice, x ∈ M, v, w ∈ 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 1 n in the discrete scheme, so there is a factor two that differs from the definition of w i in (21). Now define alsom i ,w i ,ṽ i , using rk instead of exp and log, that is (see figure 10): Finally let x n = exp(w), v n = (−1) n n log x n (z n ) and their approximationsx n = rk(m n ,w n 2 , 1 n ) andṽ n = n · rk −1 x n (z n ). We will prove that this approximate sequence converges to the true parallel transport of v: Theorem 4 Let (ṽ 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 Proof It suffices to show that there exists β > 0 such that for n large enough ṽ n − v n ≤ β n 2 . Indeed, by (22), for n large enough: 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 below, 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 figure 11 for a visual intuition of those lemmas.
Lemma 1 ∀x,x ∈ M such that x andx are close enough, for all v ∈ T x M,ṽ ∈ TṽM such that both v and δv = Π x xṽ − v are small enough, we have: . By the definition of the double exp and δv,xṽ = exp x (h x (δx, v + δv)). Then by the definition of the neighboring log, we have Using the Taylor approximations (1) and (2) truncated at the order of v , we obtain and the result follows.
Lemma 2 ∀x,x, z ∈ M close enough to one another, we have in the metric norm Proof The proof is similar to that of lemma 1 except that this time it is the norm of δv that needs to be bounded by d(x v ,xṽ). Fig. 11 Visualization for the two lemmas: lemma 1 seeks to bound the norm of d given δx and δv while lemma 2 bounds the norm of δv given d and δx. Here h is the shorthand for hx(δx, v + δv) and d for log xv (xṽ). Now, we first show that the sequence δ i = d(z i , z i ) i verifies an inductive relation, such that it is bounded by 1 n 3 , and then we will use lemma 2 to conclude. We first write 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 ,m i ) ≤ τ 2 n 4 by hypothesis on the scheme (25). Suppose for a proof by induction on i that . This is verified for i = 0 and allows to apply lemma 2 for n large enough, so Furthermore, by lemma 1 applied to v = − log m i (z i ) andṽ = −rk −1 m i (z i ), which are sufficiently small by the induction hypothesis, And combining the two results, we obtain δ i+1 ≤ (2n+1)τ 2 +nτ 3 n 5 + δ i , which completes the induction for δ i . For d(z i , m i ) we have: In the section on SL, we proved that v i ≤ 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, . . . , n − 1, and using δ 0 ≤ τ 2 n 5 by (24), we thus have δ n ≤ τ 3 +2τ 2 n 3 + n+1 n 5 τ 2 . We finally apply lemma 2 to the scaled v n ,ṽ n , as x n andx n are close enough: So that for n large enough ṽ n − v n ≤ β n 2 for some β > 0.

Remark 3
To justify the hypothesis on rk −1 , namely, rk −1 x (y) − log x (y) 2 ≤ τ n 5 , we consider the problem (P) which corresponds to an energy minimization. Working in a convex neighborhood, it admits a unique minimizer v * : The constraint can be written with Lagrange multipliers, 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 ≤ rk 1 (x, v, h)−y 2 +τ 1 h 5 , (P ) is equivalent to which is solved by gradient descent (GD) until a convergence tolerance ≤ h 5 is reached.

Numerical Simulations
It is not relevant to compare the PL with SL on spheres or SPD matrices as in the previous section, as theses spaces are symmetric and thus the PL is exact [9]. We therefore focus on the Lie group of isometries of R 3 , endowed with a left-invariant metric g with the diagonal matrix at identity: 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 β > 0 is a coefficient of anisotropy. These metrics were considered in [31] in relation with kinematics, and visualisation of the geodesics was provided, but no result on the curvature was given. Following [19], we compute explicitly the covariant derivative of the curvature and deduce (proof given in appendix E)  This allows to observe the impact of ∇R on the convergence of the PL. In the case where β = 1, the Riemannian manifold (SE(3), g) corresponds to the direct product SO(3) × 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 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 β = 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 on Fig. 12. Firstly, we observe very precisely the quadratic convergence of the infinitesimal scheme, as straight lines are obtained when plotting the error against 1 n 2 . Secondly, we see how the slope varies with β: accordingly with our results, it cancels for β = 1 and increases as β 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 β changes) such that R(w, v)v = 0 when β = 1 and ∇ w R(w, v)w = 0 when β = 1. Indeed, the former condition implies that v, w are infinitesimal rotations, which implies ∀β, ∇ 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 on Fig. 13. Note that due to the larger number of operations required for SL, our implementation is less stable and diverges when n grows too much (∼ 50). As expected when β > 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.

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 , totalling 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( √ n) steps are needed to reach a given accuracy, where O(n) are required for the FS. Comparing the values of n → 4 * 11 √ n+1 and n → 2 * 3 n , we conclude that PL is the cheapest option as soon as we want to achieve a relative accuracy better than 2.10 −2 . Note that this doesn't take into account the constants β, 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.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 ∼ 2.10 −5 with n = 37.

Conclusion
In this paper, we jointly analysed the behaviour 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 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 experiments were performed on SE(3) with anisotropic metric, and allowed to observe the role of ∇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 [22,23], or symplectic integrators [3,10] 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 → exp x (− log x (y)) is necessary and may be computed more accurately with a symmetric scheme.

C Geometry of the sphere
The unit-sphere is defined as S 2 = {x ∈ R 3 , x = 1}. With the canonical scalar product < ·, · > of R 3 , it is a Riemannian manifold of constant curvature whose geodesics are great circles. The tangent space of any x ∈ S 2 is the set of vectors orthogonal to x: TxS 2 = {v ∈ R 3 , < x, v >= 0}. We have ∀x, y ∈ S 2 , ∀v, w ∈ TxS 2 : exp x (w) = cos( w )x + sin( w ) w w , log x (y) = arccos(< y, x >) y− < y, x > x y− < y, x > x , The expression for parallel transport of v means that the orthogonal projection of v on {w} ⊥ is preserved, while the orthogonal projection on Rw 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 The tangent space of any Σ ∈ SP D(3) is the set symmetric matrices: The affine-invariant (AI) metric is defined at any Σ ∈ SP D (3), for all V, W ∈ Sym(3) using the matrix trace tr by where when not indexed, exp and log refer to the matrix operators. Finally, let The parallel transport from Σ along the geodesic with initial velocity

E Left-invariant metric on SE(3)
In this appendix, we describe the geometry of the Lie group of isometries of 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 ≥ 2, but we only detail them for d = 3 to keep them tractable. For the details on this section, we refer the reader to [6,14,19]. SE (3), is the semi-direct product of the group of three-dimensional rotations SO(3) with R 3 , i.e. the group multiplicative law for R, R ∈ SO(3), t, t ∈ R 3 is given by (R, t) · (R , t ) = (RR , t + Rt ).
It can be seen as a subgroup of GL(4) and represented by homogeneous coordinates: (R, t) = R t 0 1 , and all group operations then correspond to the matrix operations. Let the metric matrix at the identity be diagonal: G = diag(1, 1, 1, β, 1, 1) for some β > 0, the anisotropy parameter. We write < ·, · > for the associated inner-product at the identity. An orthonormal basis of the Lie algebra se (3) is and all others that cannot be deduced by skew-symmetry of the bracket are equal to 0. Extend the inner-produt defined in the Lie algebra by G to a left-invariant metric g on SE(3). Let ∇ 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 and can thus be computed explicitly thanks to (37). Let Γ k ij =< ∇e i e j , e k > be the associated Christoffel symbols. Let τ = ( √ β + 1 √ β and all other are null. Finally, recall that the curvature tensor and its covariant derivative at the identity can be defined ∀u, v, w, z ∈ se(3) by R(u, v) = ∇u∇v − ∇v∇u − ∇ [u,v] (∇uR)(v, w)z = ∇u R(v, w)z − R ∇uv, w z − R v, ∇uw z − R v, w ∇uz.

E.1.2 Case β = 1
Focusing on the case β = 1 allows to validate our implementation by testing that we recover the direct product exponential map. This is proved in e.g. [31]. It is straightforward to see that, in this case, g coincides with the (direct) product of the canonical (bi-invariant) metrics grot, gtrans of SO(d) and R d . So we have the general formula for geodesics from (R, v) ∈ SE(3) with initial velocity (Ω, u) ∈ T (R,v) SE(3): γ(t) = (R exp(tR T Ω), tu + v), and exp (R,v) (Ω, u) = (R exp(R T Ω), u + v), These are in fact valid in SE(d) for any d ≥ 2. These geodesics are represented on Fig. 14 for different values of β in SE(2) (where the metric matrix at identity is G = diag(1, β, 1)).