# Quartic $L^{p}$-convergence of cubic Riemannian splines

Hanne Hardering  Benedikt Wirth

###### Abstract

We prove quartic convergence of cubic spline interpolation for curves into
Riemannian manifolds as the grid size of the interpolation grid tends to zero.
In contrast to cubic spline interpolation in Euclidean space, where this
result is classical, the interpolation operator is no longer linear. Still,
concepts from the linear setting may be generalized to the Riemannian case,
where we try to use intrinsic Riemannian formulations and avoid charts as much
as possible.

##  1 Introduction

Let a smooth function $\gamma:[0,1]\to\mathbb{R}$ be given as well as a grid
on $[0,1]$ of knots $t_{i}=ih$, $i=0,\ldots,N$, where $h=\frac{1}{N}$ denotes
the knot distance or grid size. A standard approximation of $\gamma$ is
obtained by evaluating $\gamma$ at the knots and then computing the cubic
spline interpolant $\gamma_{h}:[0,1]\to\mathbb{R}$ of the resulting data
$(t_{i},\gamma(t_{i}))$, $i=0,\ldots,N$, for instance with Hermite boundary
conditions $\dot{\gamma}_{h}(0)=\dot{\gamma}(0)$,
$\dot{\gamma}_{h}(1)=\dot{\gamma}(1)$. It is a classic result that this
approximation is of quartic order, $\|\gamma_{h}-\gamma\|_{L^{\infty}}\leq
ch^{4}$, proved for instance by de Boor [2] (analogous convergence results
hold for spline interpolation of different order).

The well-known minimum curvature property of cubic splines states that cubic
spline interpolants uniquely minimize the accumulated squared acceleration
$\int_{0}^{1}|\ddot{\gamma}_{h}(t)|^{2}\,\mathrm{d}t$ among all functions
interpolating $(t_{i},\gamma(t_{i}))$, $i=0,\ldots,N$, which allows to
generalize the notion of a cubic spline interpolant to functions or curves
$\gamma:[0,1]\to{\mathcal{M}}$ with values in a smooth Riemannian manifold
${\mathcal{M}}$.

######  Definition 1 (Riemannian cubic spline interpolation).

Let $\gamma:[0,1]\to{\mathcal{M}}$ be a continuous curve and set $t_{i}=ih$
for $i=0,\ldots,N$ and $h=\frac{1}{N}$. A corresponding _cubic spline
interpolation_ of $\gamma$ is defined as

$\gamma_{h}\in{\mathrm{argmin}}\left\\{\int_{0}^{1}|D_{t}^{2}\gamma_{h}|^{2}\,\mathrm{d}t\,\middle|\,\gamma_{h}(t_{i})=\gamma(t_{i})\text{
for }i=0,\ldots,N\text{ and }\dot{\gamma}_{h}(t)=\dot{\gamma}(t)\text{ for
}t=0,1\right\\}.$

Above, $D_{t}^{2}\gamma_{h}$ denotes the intrinsic second derivative, the
covariant derivative of $\dot{\gamma}_{h}$ along the curve $\gamma_{h}$, and
$|\cdot|$ represents the Riemannian norm on the tangent space bundle
$T{\mathcal{M}}$. A natural question is whether the approximation properties
of cubic splines transfer from the setting of real-valued functions to the
manifold-valued generalization. Our main result is that they do.

######  Theorem 2 (Quartic convergence of Riemannian cubic spline
interpolation).

Let ${\mathcal{M}}$ be complete with bounded Riemann curvature and its
derivatives. Let $\gamma:[0,1]\to{\mathcal{M}}$ be four times differentiable.
For $h$ small enough (depending on ${\mathcal{M}}$ and $\gamma$)

  1. 1.

the cubic spline interpolation $\gamma_{h}$ exists and is unique,

  2. 2.

and it satisfies

$\|d(\gamma_{h},\gamma)\|_{L^{p}}\leq ch^{4},$

where the constant $c>0$ only depends on ${\mathcal{M}}$, $\gamma$, and
$p\in[1,\infty]$.

Above, $d(\cdot,\cdot)$ denotes the Riemannian distance and
$\|d(\gamma_{h},\gamma)\|_{L^{p}}$ the $L^{p}$-norm of the function $t\mapsto
d(\gamma_{h}(t),\gamma(t))$ for $t\in[0,1]$. Note that the conditions for
existence and uniqueness can be a little relaxed, see proposition 9. One could
envision two strategies to arrive at this approximation result. One strategy
would be to reduce the problem to the real-valued case by reformulating the
Riemannian cubic spline interpolation as an interpolation problem in a linear
space (using charts or a tangent space to the manifold) which can be viewed as
a perturbation of Euclidean or real-valued cubic spline interpolation. The
complementary strategy consists in the generalization of all ingredients to
the Euclidean analysis to the Riemannian setting, using a mostly coordinate-
free intrinsic Riemannian formulation. We will follow the latter strategy (but
also discuss the former along the way). In particular, the Euclidean proof in
[2] makes use of the fact that the second derivative of the cubic spline
interpolant $\gamma_{h}$ is the $L^{2}$-best approximation to the the second
derivative of $\gamma$, and for our proof in the nonlinear manifold setting we
will quantify the deviation from a Riemannian generalization of that property.
This will result in an $L^{\infty}$-bound for the second derivative
$D_{t}^{2}\gamma_{h}$, which is crucial for the convergence estimate.

For the sake of completeness, along the way we will deal with linear spline
interpolation since it is a widespread alternative for applications in
Riemannian manifolds and since its convergence properties follow without
additional effort.

######  Definition 3 (Riemannian linear spline interpolation).

Let $\gamma:[0,1]\to{\mathcal{M}}$ be a continuous curve and set $t_{i}=ih$
for $i=0,\ldots,N$ and $h=\frac{1}{N}$. A corresponding _linear spline
interpolation_ of $\gamma$ is defined as

$\gamma_{h}\in{\mathrm{argmin}}\left\\{\int_{0}^{1}|\dot{\gamma}_{h}|^{2}\,\mathrm{d}t\,\middle|\,\gamma_{h}(t_{i})=\gamma(t_{i})\text{
for }i=0,\ldots,N\right\\}.$

######  Theorem 4 (Quadratic convergence of Riemannian linear spline
interpolation).

Let ${\mathcal{M}}$ be complete with bounded Riemann curvature and its
derivatives. Let $\gamma:[0,1]\to{\mathcal{M}}$ be twice differentiable. For
$h$ small enough (depending on ${\mathcal{M}}$ and $\gamma$)

  1. 1.

the linear spline interpolation $\gamma_{h}$ exists and is unique,

  2. 2.

and it satisfies

$\|d(\gamma_{h},\gamma)\|_{L^{p}}\leq ch^{2},$

where the constant $c>0$ only depends on ${\mathcal{M}}$, $\gamma$, and
$p\in[1,\infty]$.

Higher order spline interpolation of odd degree could in principle be tackled
in a similar way as cubic spline interpolation, but of course requires more
technical estimates.

In section 2 we will summarize the required basic Riemannian notions and our
notation after which we prove the well-posedness of (linear and) cubic spline
interpolation in section 3 and derive the corresponding optimality conditions.
We then present two convergence proofs for the classical Euclidean setting in
section 4 before adapting these to the Riemannian setting in section 5.

##  2 Required Riemannian notions and preliminaries

Throughout we will assume ${\mathcal{M}}$ to be a smooth, complete, connected,
finite-dimensional Riemannian manifold with Riemannian metric $g$ such that
the Riemann curvature tensor and all its derivatives are bounded. The tangent
space at a point $p\in{\mathcal{M}}$ will be denoted $T_{p}{\mathcal{M}}$, and
for $v,w\in T_{p}{\mathcal{M}}$ we will write

$(v,w)=g_{p}(v,w)\qquad\text{as well as}\qquad|v|=\sqrt{(v,v)},$

where the explicit dependence on the base point $p$ is suppressed for better
readability. The induced Riemannian distance is denoted by
$d:{\mathcal{M}}\times{\mathcal{M}}\to[0,\infty)$, the Levi-Civita connection
on ${\mathcal{M}}$ is denoted by $\nabla$, and the Riemann curvature tensor is
denoted by ${\mathrm{R}}$. For an absolutely continuous curve
$\gamma:[0,1]\to{\mathcal{M}}$ we will denote the parameterization variable by
$t$ and its velocity by $\partial_{t}\gamma(t)=\dot{\gamma}(t)\in
T_{\gamma(t)}{\mathcal{M}}$ (time derivatives of functions into a vector space
will also be denoted with powers of $\partial_{t}$ or dots). A lift of the
curve $\gamma$ is a mapping $v$ from $[0,1]$ into the tangent bundle
$T{\mathcal{M}}$ such that $v(t)\in T_{\gamma(t)}{\mathcal{M}}$ for all
$t\in[0,1]$; the vector space of all liftings is abbreviated
$\gamma^{\\#}T{\mathcal{M}}$. For a vector field
$v\in\gamma^{\\#}T{\mathcal{M}}$ along the curve $\gamma$ and $k\in\mathbb{N}$
we will write

$D_{t}^{k}v=\nabla_{\dot{\gamma}}^{k}v,$

where with a slight misuse of notation the covariant derivative of $v$ along
$\gamma$ shall be

$\nabla_{\dot{\gamma}}v(t)=\lim_{\epsilon\to
0}\tfrac{1}{\epsilon}(\pi_{\gamma}v(t+\epsilon)-v(t))$

with $\pi_{\gamma}$ being the parallel transport along $\gamma$ from
$\gamma(t+\epsilon)$ to $\gamma(t)$ (start and end point of the parallel
transport along a curve will usually be clear from the context). We will also
write

$D_{t}^{k}\gamma=\nabla_{\dot{\gamma}}^{k-1}\dot{\gamma}$

for $k\in\mathbb{N}$. A family of curves
$\gamma(s,\cdot):[0,1]\to{\mathcal{M}}$ will be parameterized by the variable
$s\in\mathbb{R}$, and analogously to before, for a family of vector fields
$v(s,\cdot)$ along the curves $\gamma(s,\cdot)$ we introduce the notation

$D_{s}v(s,t)=\nabla_{\partial_{s}\gamma}v(s,t)=\lim_{\epsilon\to
0}\tfrac{1}{\epsilon}(\pi_{\gamma(\cdot,t)}v(s+\epsilon,t)-v(s,t)),$

where $\pi_{\gamma(\cdot,t)}$ is parallel transport along the curve
$s\mapsto\gamma(s,t)$. Below we collect a few classical Riemannian calculus
rules to be used in the sequel.

######  Lemma 5 (Differentiation rules).

  1. 1.

Parallel transport along a curve commutes with differentiation along the curve
in the following sense. Let $\gamma:[0,1]\to{\mathcal{M}}$ be differentiable
and denote by $\pi_{\gamma}:T_{\gamma(t)}{\mathcal{M}}\to
T_{\gamma(0)}{\mathcal{M}}$ the parallel transport along $\gamma$, then for
any differentiable $v\in\gamma^{\\#}T{\mathcal{M}}$ we have

$\tfrac{\mathrm{d}}{\mathrm{d}t}\pi_{\gamma}v=\pi_{\gamma}D_{t}v.$

  2. 2.

Let $\gamma:\mathbb{R}\times[0,1]\to{\mathcal{M}}$ be a differentiable family
of differentiable curves, then
$D_{s}\partial_{t}\gamma=D_{t}\partial_{s}\gamma$.

  3. 3.

Let $\gamma:\mathbb{R}\times[0,1]\to{\mathcal{M}}$ be a differentiable family
of differentiable curves and $v\in\gamma^{\\#}T{\mathcal{M}}$ differentiable,
then
$D_{s}D_{t}v=D_{t}D_{s}v+{\mathrm{R}}(\partial_{s}\gamma,\partial_{t}\gamma)v$.

###### Proof.

  1. 1.

For all $w\in T_{\gamma(0)}{\mathcal{M}}$ we have

$(\tfrac{\mathrm{d}}{\mathrm{d}t}\pi_{\gamma}v,w)=\tfrac{\mathrm{d}}{\mathrm{d}t}(\pi_{\gamma}v,w)=\tfrac{\mathrm{d}}{\mathrm{d}t}(v,\pi_{\gamma}^{-1}w)=(D_{t}v,\pi_{\gamma}^{-1}w)+(v,D_{t}\pi_{\gamma}^{-1}w)=(D_{t}v,\pi_{\gamma}^{-1}w)=(\pi_{\gamma}D_{t}v,w).$

  2. 2.

Since the Levi-Civita connection is torsion-free, we have
$D_{s}\partial_{t}\gamma=D_{t}\partial_{s}\gamma+[\partial_{s}\gamma,\partial_{t}\gamma]$
for the Lie bracket $[\cdot,\cdot]$ of vector fields. Now for any smooth
function $f:{\mathcal{M}}\to\mathbb{R}$ we have

$[\partial_{s}\gamma,\partial_{t}\gamma]f=\partial_{s}\gamma(\partial_{t}\gamma(f))-\partial_{t}\gamma(\partial_{s}\gamma(f))=\partial_{s}\partial_{t}f\circ\gamma-\partial_{t}\partial_{s}f\circ\gamma=0.$

  3. 3.

We have

${\mathrm{R}}(\partial_{s}\gamma,\partial_{t}\gamma)v=D_{s}D_{t}v-D_{t}D_{s}v-\nabla_{[\partial_{s}\gamma,\partial_{t}\gamma]}v,$

where the last term is zero as shown in the previous point. ∎

For a Lebesgue measurable function $v:I\to T{\mathcal{M}}$ on some measurable
domain $I$ we denote its $L^{p}$-norm for $p\in[1,\infty]$ by

$\|v\|_{L^{p}}=\|\,|v|\,\|_{L^{p}(I)}.$

If the domain is not clear from the context, we will write $\|v\|_{L^{p}(I)}$
instead. Similarly, for a function $v\in\gamma^{\\#}T{\mathcal{M}}$ we
introduce its Hölder norm of exponent $\alpha\in(0,1]$ by

$\|v\|_{C^{0,\alpha}}=\|v\|_{L^{\infty}}+\sup_{r,t\in[0,1]}\frac{|\pi_{\gamma}v(t)-\pi_{\gamma}v(r)|}{|t-r|^{\alpha}},$

where $\pi_{\gamma}$ is parallel transport to $\gamma(0)$ along $\gamma$.

Finally, the Riemannian exponential with base point $p\in{\mathcal{M}}$ is
denoted by $\exp_{p}:T_{p}{\mathcal{M}}\to{\mathcal{M}}$. Letting $\rho$ be
the injectivity radius of ${\mathcal{M}}$, the inverse of $\exp_{p}$ can be
defined on the geodesic ball $B_{\rho}(p)\subset{\mathcal{M}}$ of radius
$\rho$ and is called the Riemannian logarithm $\log_{p}:B_{\rho}(p)\to
T_{p}{\mathcal{M}}$. We close by providing some bounds on the Riemannian
logarithm and exponential.

######  Lemma 6 (Bounds on Riemannian logarithm and exponential).

Let ${\mathcal{M}}$ be a $C^{k+1}$-manifold, $k\geq 2$, with injectivity
radius $\rho$, sectional curvature bounded from above by $K$, and bounds
$\|{\mathrm{R}}\|_{\infty}$ on ${\mathrm{R}}$ and
$\|\nabla{\mathrm{R}}\|_{\infty}$ on $\nabla{\mathrm{R}}$. Then the bivariate
mapping $\log(p,q)=\log_{p}q$ is in $C^{k}$ for $d(p,q)<\rho$. If $K>0$,
assume additionally $d(p,q)\leq\frac{\pi}{2\sqrt{K}}$. The operator norms of
the first derivatives of $\log$ satisfy

$\displaystyle\|\partial_{1}\log_{p}q+\mathrm{Id}_{T_{p}M}\|$
$\displaystyle\leq\frac{\|{\mathrm{R}}\|_{L^{\infty}}}{2}d^{2}(p,q),$
$\displaystyle\|\partial_{2}\log_{p}q-\pi_{q\to p}\|$
$\displaystyle\leq\frac{\|{\mathrm{R}}\|_{L^{\infty}}}{2}d^{2}(p,q),$

where $\pi_{q\to p}$ denotes parallel transport along the geodesic from $q$ to
$p$. The covariant second order derivatives, denoted by $D^{2}$, satisfy

$\displaystyle\|D^{2}\log_{p}q\|\leq
C(\|{\mathrm{R}}\|_{\infty},\|\nabla{\mathrm{R}}\|_{\infty})d(p,q)$

for a constant $C(\|{\mathrm{R}}\|_{\infty},\|\nabla{\mathrm{R}}\|_{\infty})$
depending only on $\|{\mathrm{R}}\|_{\infty}$ and
$\|\nabla{\mathrm{R}}\|_{\infty}$. Similarly, the bivariate mapping
$\exp(p,v)=\exp_{p}v$ is in $C^{k}$, and for $|v|$ small enough, depending
only on $\rho$, $K$, and $\|{\mathrm{R}}\|_{\infty}$,
$\|\nabla{\mathrm{R}}\|_{\infty}$, it satisfies

$\displaystyle\|\partial_{2}\exp_{p}v-\pi_{p\to q}\|$ $\displaystyle\leq
C(\|{\mathrm{R}}\|_{\infty},\|\nabla{\mathrm{R}}\|_{\infty})|v|^{2},$
$\displaystyle\|\partial_{2}^{2}\exp_{p}v\|$ $\displaystyle\leq
C(\|{\mathrm{R}}\|_{\infty},\|\nabla{\mathrm{R}}\|_{\infty})|v|.$

The proof for the derivatives of the logarithm follows from Jacobi field
estimates and can be found in [5, Prop. A.1-A.2] (the estimates for the first
order derivatives can also be found in [7]). The estimates for the derivatives
of the exponential are then straightforward applications of the inverse
function theorem.

##  3 Well-posedness and Euler–Lagrange equations of Riemannian spline
interpolations

All throughout the article, $\gamma:[0,1]\to{\mathcal{M}}$ will be a given
curve, $k$ times differentiable ($k=2$ for linear and $k=4$ for cubic spline
interpolation) in the sense that $\gamma$ is continuous and $D_{t}^{l}\gamma$
is well-defined and bounded for $l=1,\ldots,k$. Furthermore,
$\gamma_{h}:[0,1]\to{\mathcal{M}}$ will denote the spline interpolation
according to definition 1 or definition 3 at points $t_{i}=ih$,
$i=0,\ldots,N$, where $h=\frac{1}{N}$.

######  Remark 7 (Boundary conditions for Riemannian cubic spline
interpolation).

In definition 1 we chose to impose the so-called Hermite boundary conditions
$\dot{\gamma}_{h}(t)=\dot{\gamma}(t)\text{ for }t=0,1$. An alternative would
be to solely require $\gamma_{h}=\gamma$ at all interpolation points. This is
known to result in so-called natural boundary conditions, essentially a
vanishing acceleration of $\gamma_{h}$ at $t=0$ and $t=1$. Unfortunately,
though, existence of cubic spline interpolations with natural boundary
conditions cannot be guaranteed due to lacking control of the curve velocity
(a simple example is provided in [6]).

We begin with straightforward a priori bounds that follow from the coercivity
of the spline interpolation energy, which is nothing else but
$\|\dot{\gamma}_{h}\|_{L^{2}}^{2}$ in the case of linear splines and
$\|D_{t}^{2}\gamma_{h}\|_{L^{2}}^{2}$ in case of cubic splines. Boundedness of
higher derivatives for the cubic spline interpolation will be proven later.

######  Proposition 8 (A priori bounds).

Let $\gamma:[0,1]\to{\mathcal{M}}$ be an absolutely continuous curve.

  1. 1.

If $\|\dot{\gamma}\|_{L^{2}}^{2}$ is bounded, then $\gamma$ is Hölder
continuous with exponent $\frac{1}{2}$, and for any $r\in[0,1]$ we have

$\|d(\gamma,\gamma(r))\|_{C^{0,\frac{1}{2}}}\leq 2\|\dot{\gamma}\|_{L^{2}}.$

  2. 2.

If $\|D_{t}^{2}\gamma\|_{L^{2}}^{2}$ is bounded, then $\gamma$ is even Hölder
continuously differentiable with exponent $\frac{1}{2}$ and satisfies

$\displaystyle\|\dot{\gamma}\|_{C^{0,\frac{1}{2}}}$
$\displaystyle\leq|\dot{\gamma}(0)|+2\|D_{t}^{2}\gamma\|_{L^{2}},$
$\displaystyle d(\gamma(t),\gamma(r))$
$\displaystyle\leq|t-r|\|\dot{\gamma}\|_{L^{\infty}}\quad\text{for all
}r,t\in[0,1].$

###### Proof.

  1. 1.

For any $t,r\in[0,1]$, by Hölder’s inequality we have

$d(\gamma(t),\gamma(r))\leq\int_{r}^{t}|\dot{\gamma}(\tau)|\,\mathrm{d}\tau\leq\||\dot{\gamma}|\|_{L^{2}}\sqrt{|t-r|}=\|\dot{\gamma}\|_{L^{2}}\sqrt{|t-r|}$

so that $\gamma$ is indeed Hölder continuous with exponent $2$. Furthermore,
this implies

$\sup_{t\in[0,1]}d(\gamma(t),\gamma(r))+\sup_{t,\tilde{t}\in[0,1]}\frac{|d(\gamma(t),\gamma(r))-d(\gamma(\tilde{t}),\gamma(r))|}{\sqrt{|t-\tilde{t}|}}\leq\|\dot{\gamma}\|_{L^{2}}+\sup_{t,\tilde{t}\in[0,1]}\frac{d(\gamma(t),\gamma(\tilde{t}))}{\sqrt{|t-\tilde{t}|}}\leq
2\|\dot{\gamma}\|_{L^{2}},$

from which the estimate follows.

  2. 2.

Letting $\pi_{\gamma}$ denote parallel transport along $\gamma$ to
$T_{\gamma(r)}{\mathcal{M}}$ we can analogously estimate

$|\pi_{\gamma}\dot{\gamma}(t)-\dot{\gamma}(r)|=\left|\int_{r}^{t}\frac{\mathrm{d}}{\mathrm{d}\tau}\pi_{\gamma}\dot{\gamma}(\tau)\,\mathrm{d}\tau\right|=\left|\int_{r}^{t}\pi_{\gamma}D_{\tau}\dot{\gamma}(\tau)\,\mathrm{d}\tau\right|\leq\sqrt{|t-r|}\|\pi_{\gamma}D_{t}\dot{\gamma}\|_{L^{2}}=\sqrt{|t-r|}\|D_{t}^{2}\gamma\|_{L^{2}}$

by lemma 5(1) and Hölder’s inequality. Thus

$\sup_{t\in[0,1]}|\dot{\gamma}(t)|+\sup_{t,r\in[0,1]}\frac{|\pi_{\gamma}\dot{\gamma}(t)-\dot{\gamma}(r)|}{\sqrt{|t-r|}}\leq|\dot{\gamma}(0)|+\sqrt{t}\|D_{t}^{2}\gamma\|_{L^{2}}+\|D_{t}^{2}\gamma\|_{L^{2}},$

which is the first norm bound. The second bound trivially follows from
$d(\gamma(t),\gamma(r))\leq\int_{r}^{t}|\dot{\gamma}(\tau)|\,\mathrm{d}\tau$.
∎

Obviously, an immediate consequence is that the linear spline interpolation
$\gamma_{h}$ of $\gamma$, if it exists, is Hölder continuous with

$d(\gamma_{h}(t),\gamma_{h}(t_{i}))\leq
2\sqrt{|t-t_{i}|}\|\dot{\gamma}\|_{L^{2}}\quad\text{for all
}t\in[0,1],i\in\\{1,\ldots,N\\}$ (1)

and that the cubic spline interpolation, if it exists, is Hölder continuously
differentiable with

$\displaystyle\|\dot{\gamma}_{h}\|_{C^{0,\frac{1}{2}}}$
$\displaystyle\leq|\dot{\gamma}(0)|+2\|D_{t}^{2}\gamma\|_{L^{2}},$ (2)
$\displaystyle d(\gamma_{h}(t),\gamma(t_{i}))$
$\displaystyle\leq|t-t_{i}|\|\dot{\gamma}\|_{L^{\infty}}\quad\text{for all
}t\in[0,1],i\in\\{1,\ldots,N\\}.$ (3)

We proceed by proving the well-posedness of spline interpolation, theorem 4(1)
and theorem 2(1).

######  Proposition 9 (Well-posedness of Riemannian spline interpolation).

Linear and cubic spline interpolations $\gamma_{h}$ of
$\gamma:[0,1]\to{\mathcal{M}}$ exist. Furthermore, they are unique if $\gamma$
is continuous (for linear spline interpolation) or has finite
$\|D_{t}^{2}\gamma\|_{L^{2}}$ (for cubic spline interpolation) and $h$ is
small enough (depending on $\gamma$ and ${\mathcal{M}}$).

###### Proof.

Obviously, a linear spline interpolation is obtained by minimizing for each
$i=1,\ldots,N$ the energy

$\int_{t_{i-1}}^{t_{i}}|\dot{\gamma}_{h}|^{2}\,\mathrm{d}t$

among all curve segments $\gamma_{h}:[t_{i-1},t_{i}]\to{\mathcal{M}}$ that
satisfy $\gamma_{h}(t_{i-1})=\gamma(t_{i-1})$ and
$\gamma_{h}(t_{i})=\gamma(t_{i})$. Those curve segments are known to be
exactly the constant speed-parameterized geodesics connecting
$\gamma(t_{i-1})$ and $\gamma(t_{i})$, which always exist on a complete
Riemannian manifold and which are unique if the end points are close enough to
each other (depending on ${\mathcal{M}}$). By uniform continuity of $\gamma$
the latter condition will be fulfilled for all $i=1,\ldots,N$ as soon as $h$
is small enough.

As for the well-posedness of the cubic spline interpolation, consider a
minimizing sequence $\gamma_{h}^{n}$, $n=1,2,\ldots$, of curves satisfying the
interpolation conditions such that the cubic spline energy
$E[\gamma_{h}^{n}]=\|D_{t}^{2}\gamma_{h}^{n}\|_{L^{2}}^{2}$ converges
monotonically to the infimum,
$E[\gamma_{h}^{n}]\to\inf\\{E[\gamma_{h}]\,|\,\gamma_{h}(t_{i})=\gamma(t_{i})\text{
for }i=0,\ldots,N\text{ and }\dot{\gamma}_{h}(t)=\dot{\gamma}(t)\text{ for
}t=0,1\\}$. Letting $\pi_{\tilde{\gamma}}$ denote parallel transport to
$T_{\tilde{\gamma}(0)}{\mathcal{M}}$ along a curve $\tilde{\gamma}$, we have

$\|\tfrac{\mathrm{d}}{\mathrm{d}t}\pi_{\gamma_{h}^{n}}\dot{\gamma}_{h}^{n}\|_{L^{2}}^{2}=\|\pi_{\gamma_{h}^{n}}D_{t}^{2}\gamma_{h}^{n}\|_{L^{2}}^{2}=\|D_{t}^{2}\gamma_{h}^{n}\|_{L^{2}}^{2}=E[\gamma_{h}^{n}]$

so that the functions $\pi_{\gamma_{h}^{n}}\dot{\gamma}_{h}^{n}:[0,1]\to
T_{\gamma(0)}{\mathcal{M}}$ have uniformly bounded $H^{1}$-seminorm. Together
with
$\pi_{\gamma_{h}^{n}}\dot{\gamma}_{h}^{n}(0)=\dot{\gamma}_{h}^{n}(0)=\dot{\gamma}(0)$
it follows by Poincaré’s inequality that the functions even have uniformly
bounded $H^{1}$-norm so that a subsequence converges weakly in $H^{1}$ (and
thus strongly in $C^{0}([0,1])$) to some function $\beta:[0,1]\to
T_{\gamma(0)}{\mathcal{M}}$. Now define $\gamma_{h}:[0,1]\to{\mathcal{M}}$ as
the solution to the ordinary differential equation

$\dot{\gamma}_{h}(t)=\pi_{\gamma_{h}}^{-1}\beta(t)$

with initial value $\gamma_{h}(0)=\gamma(0)$, where we write
$\pi_{\gamma_{h}}^{-1}\beta(t)$ for the parallel transport of $\beta(t)$ to
$\gamma_{h}(t)$ along $\gamma_{h}$. The solution
$\gamma_{h}:[0,1]\to{\mathcal{M}}$ exists and is unique by the theorem of
Picard and Lindelöf, and it satisfies the interpolation conditions by
Gronwall’s lemma. Indeed, let us express all quantities in local coordinates
$\mathbb{R}^{d}$ and indicate this by a hat. Further, let
$\hat{\Pi}:[0,1]\to\mathbb{R}^{d\times d}$ denote the matrix representation in
local coordinates of $\pi_{\gamma_{h}}^{-1}$. Then $\hat{\gamma}_{h}$ solves
the initial value problem

$\frac{\mathrm{d}}{\mathrm{d}t}{\hat{\gamma}_{h}\choose\hat{\Pi}}(t)=F(t,\hat{\gamma}_{h}(t),\hat{\Pi}(t)),\;{\hat{\gamma}_{h}\choose\hat{\Pi}}(0)={\hat{\gamma}_{h}(0)\choose\mathrm{Id}}\\\
\quad\text{for }F:\mathbb{R}\times\mathbb{R}^{d}\times\mathbb{R}^{d\times
d}\to\mathbb{R}^{d}\times\mathbb{R}^{d\times
d},\;F(t,\hat{\gamma}_{h},\hat{\Pi})={\hat{\Pi}\hat{\beta}(t)\choose-\Gamma_{\hat{\gamma}_{h}}(\hat{\Pi}\hat{\beta}(t),\hat{\Pi})},$

where
$\Gamma_{\hat{\gamma}_{h}}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{d}$
is the Christoffel operator at coordinates $\hat{\gamma}_{h}$ and the second
entry of $F(t,\hat{\gamma}_{h},\hat{\Pi})$ refers to the matrix defined by
$\mathbb{R}^{d}\ni\xi\mapsto\Gamma_{\hat{\gamma}_{h}}(\hat{\Pi}\hat{\beta}(t),\hat{\Pi}\xi)\in\mathbb{R}^{d}$.
The function $F$ is obviously Lipschitz so that the Picard–Lindelöf theorem
can be applied (it is straightforward to show that $\hat{\gamma}_{h}$ and
$\hat{\Pi}$ stay bounded for $t\in[0,1]$). Since $\hat{\gamma}_{h}^{n}$
satisfies the same initial value problem, only with $\hat{\beta}$ replaced by
$\widehat{\pi_{\gamma_{h}^{n}}\dot{\gamma}_{h}^{n}}$, the curves
$\hat{\gamma}_{h}^{n}$ converge uniformly on $[0,1]$ to $\hat{\gamma}_{h}$ by
Gronwall’s lemma, as do the corresponding $\hat{\alpha}$. Consequently,
$\hat{\gamma}_{h}$ satisfies the interpolation conditions. Finally note that

$E[\gamma_{h}]=\|\pi_{\gamma_{h}}D_{t}^{2}\gamma_{h}\|_{L^{2}}^{2}=\|\tfrac{\mathrm{d}}{\mathrm{d}t}\pi_{\gamma_{h}}\dot{\gamma}_{h}\|_{L^{2}}^{2}=\|\dot{\beta}\|_{L^{2}}^{2}\leq\liminf_{n\to\infty}\|\tfrac{\mathrm{d}}{\mathrm{d}t}\pi_{\gamma_{h}^{n}}\dot{\gamma}_{h}^{n}\|_{L^{2}}^{2}=\liminf_{n\to\infty}E[\gamma_{h}^{n}]$

due to the sequential weak lower semi-continuity of the $L^{2}$-norm so that
$\gamma_{h}$ minimizes $E$ under the interpolation constraints.

It remains to prove uniqueness of $\gamma_{h}$. To this end we will show local
convexity of the problem for $h$ small enough. Indeed, for $h$ small enough
the interpolation problem can equivalently be formulated in terms of functions

$v_{i}:[t_{i-1},t_{i}]\to T_{\gamma_{i}}{\mathcal{M}},\quad
v_{i}(t)=\log_{\gamma_{i}}\gamma_{h}(t),\;i=1,\ldots,N,$

where we abbreviated $\gamma_{i}=\gamma(t_{i-1})$ (equivalently,
$\gamma_{h}(t)=\exp_{\gamma_{i}}v_{i}(t)$ for $t\in[t_{i-1},t_{i}]$). This is
possible since $E[\gamma_{h}]\leq E[\gamma]$ and thus by proposition 8(2),

$d(\gamma_{h}(t),\gamma_{i})\leq h\|\dot{\gamma}\|_{L^{\infty}}\leq
h\left(|\dot{\gamma}(0)|+2\|D_{t}^{2}\gamma\|_{L^{2}}\right)\quad\text{for
}t\in[t_{i-1},t_{i}]$

so that for $h$ small enough the Riemannian logarithm in the definition of
$v_{i}$ is uniquely defined and smooth. Furthermore, again by proposition
8(2), the $v_{i}$ satisfy

$\displaystyle\|[\partial_{2}\exp_{\gamma_{i}}v_{i}]\dot{v}_{i}\|_{C^{0,\frac{1}{2}}}$
$\displaystyle\leq C\quad\text{for
}C=|\dot{\gamma}(0)|+2\|D_{t}^{2}\gamma\|_{L^{2}},$
$\displaystyle\|v_{i}\|_{L^{\infty}}$ $\displaystyle\leq Ch,$

for $i=1,\ldots,N$, which by lemma 6 implies $\|\dot{v}_{i}\|_{L^{\infty}}\leq
2C$ for $h$ small enough. Now introduce

$J[v_{1},\ldots,v_{N}]=\sum_{i=1}^{N}\int_{t_{i-1}}^{t_{i}}\left|[\partial_{2}\exp_{\gamma_{i}}v_{i}(t)]\ddot{v}_{i}(t)+[\partial_{2}^{2}\exp_{\gamma_{i}}v_{i}(t)](\dot{v}_{i}(t),\dot{v}_{i}(t))\right|^{2}\mathrm{d}t,$

which equals $E[\gamma_{h}]$. Due to the above bounds on
$\|v_{i}\|_{L^{\infty}}$ and $\|\dot{v}_{i}\|_{L^{\infty}}$ we can apply lemma
6 to obtain $E[\gamma]\geq
E[\gamma_{h}]=J[v_{1},\ldots,v_{N}]\geq\frac{1}{2}\|\ddot{v}_{i}\|_{L^{2}}^{2}-1$
for $h$ small enough. Thus, the cubic spline interpolation problem in terms of
the $v_{i}$ can thus be written as

$\displaystyle\\!\min\left\\{J[v_{1},\ldots,v_{N}]\,\middle|\,(v_{1},\ldots,v_{N})\in\mathcal{A}_{h}\right\\}\quad\text{with
}$ $\displaystyle\mathcal{A}_{h}=\big{\\{}(v_{1},\ldots,v_{N})\in
L^{\infty}((t_{0},t_{1}))\times\ldots\times
L^{\infty}((t_{N-1},t_{N}))\,\big{|}\,\|v_{i}\|_{L^{\infty}}\leq
Ch,\,\|\dot{v}_{i}\|_{L^{\infty}}\leq 2C,\,\|\ddot{v}_{i}\|_{L^{2}}^{2}\leq
2E[\gamma]+2,$ $\displaystyle\hskip
245.4159ptv_{i}(t_{i-1})=0,\,v_{i}(t_{i})=\log_{\gamma_{i}}\gamma(t_{i})\text{
for }i=1,\ldots,N$ $\displaystyle\hskip 223.88818pt\text{and
}\dot{v}_{1}(0)=\dot{\gamma}(0),\,[\partial_{2}\exp_{\gamma(t_{N-1})}v_{N}(1)]\dot{v}_{N}(1)=\dot{\gamma}(1)\big{\\}}.$

It is obvious that $\mathcal{A}_{h}$ is convex, so it remains to show
convexity of $J$ on $\mathcal{A}_{h}$. Abbreviating
$\gamma_{v}(t)=\exp_{\gamma_{i}}v_{i}(t)$ for $t\in[t_{i-1},t_{i}]$, the
Gâteaux derivatives of $J$ in a direction
$\phi\equiv(\varphi_{1},\ldots,\varphi_{N})$ with
$\varphi_{i}(t_{i-1})=\varphi(t_{i})=0$ and
$\dot{\varphi}_{1}(0)=\dot{\varphi}_{N}(1)=0$ are

$\displaystyle\partial J[v_{1},\ldots,v_{N}](\phi)$
$\displaystyle=2\sum_{i=1}^{N}\int_{t_{i-1}}^{t_{i}}(D_{t}^{2}\gamma_{v},\underbrace{[\partial_{2}\exp_{\gamma_{i}}v_{i}]\ddot{\varphi}_{i}+[\partial_{2}^{2}\exp_{\gamma_{i}}v_{i}](\ddot{v}_{i},\varphi_{i})+[\partial_{2}^{3}\exp_{\gamma_{i}}v_{i}](\dot{v}_{i},\dot{v}_{i},\varphi_{i})+2[\partial_{2}^{2}\exp_{\gamma_{i}}v_{i}](\dot{v}_{i},\dot{\varphi}_{i})}_{U_{i}})\,\mathrm{d}t,$
$\displaystyle\partial^{2}J[v_{1},\ldots,v_{N}](\phi,\phi)$
$\displaystyle=2\sum_{i=1}^{N}\int_{t_{i-1}}^{t_{i}}(U_{i},U_{i})+(D_{t}^{2}\gamma_{v},\underbrace{2[\partial_{2}^{2}\exp_{\gamma_{i}}v_{i}](\varphi_{i},\ddot{\varphi}_{i})+[\partial_{2}^{3}\exp_{\gamma_{i}}v_{i}](\varphi_{i},\varphi_{i},\ddot{v}_{i})+4[\partial_{2}^{3}\exp_{\gamma_{i}}v_{i}](\dot{\varphi}_{i},\dot{v}_{i},\varphi_{i})}_{V_{i}^{1}}$
$\displaystyle\hskip
236.80481pt+\underbrace{[\partial_{2}^{4}\exp_{\gamma_{i}}v_{i}](\varphi_{i},\varphi_{i},\dot{v}_{i},\dot{v}_{i})+2[\partial_{2}^{2}\exp_{\gamma_{i}}v_{i}](\dot{\varphi}_{i},\dot{\varphi}_{i})}_{V_{i}^{2}})\,\mathrm{d}t.$

Now by lemma 6,

$\displaystyle\|U_{i}-\ddot{\varphi}_{i}\|_{L^{2}}$ $\displaystyle\leq
K\|v_{i}\|_{L^{\infty}}\left(\|v_{i}\|_{L^{\infty}}\|\ddot{\varphi}_{i}\|_{L^{2}}+\|\ddot{v}_{i}\|_{L^{2}}\|\varphi_{i}\|_{L^{\infty}}+\|\dot{v}_{i}\|_{L^{\infty}}\|\dot{\varphi}_{i}\|_{L^{2}}\right)+K\|\dot{v}_{i}\|_{L^{\infty}}^{2}\|\dot{\varphi}_{i}\|_{L^{2}},$
$\displaystyle\|V_{i}^{1}+V_{i}^{2}\|_{L^{2}}$ $\displaystyle\leq
K\|v_{i}\|_{L^{\infty}}\left(\|\varphi_{i}\|_{L^{\infty}}\|\ddot{\varphi}_{i}\|_{L^{2}}+\|\dot{\varphi}_{i}\|_{L^{\infty}}\|\dot{\varphi}_{i}\|_{L^{2}}\right)$
$\displaystyle\quad+K\left(\|\varphi_{i}\|_{L^{\infty}}^{2}\|\ddot{v}_{i}\|_{L^{2}}+\|\varphi_{i}\|_{L^{\infty}}\|\dot{v}_{i}\|_{L^{\infty}}\|\ddot{\varphi}_{i}\|_{L^{2}}+\|\dot{v}_{i}\|_{L^{\infty}}^{2}\|\varphi_{i}\|_{L^{\infty}}\|\varphi_{i}\|_{L^{2}}\right)$

for a constant $K>0$ depending only on ${\mathcal{M}}$. Now recall that for
$(v_{1},\ldots,v_{N})\in\mathcal{A}_{h}$ we have

$\|v_{i}\|_{L^{\infty}}\leq Ch,\quad\|\dot{v}_{i}\|_{L^{\infty}}\leq
2C,\quad\|\ddot{v}_{i}\|_{L^{2}}^{2}\leq 2\|D_{t}^{2}\gamma\|_{L^{2}}^{2}+2.$

Furthermore, due to the homogeneous boundary conditions on each $\varphi_{i}$
there exists $c>0$ independent of $h$ or $\varphi_{i}$ such that

$\|\varphi_{i}\|_{L^{\infty}}\leq
ch^{3/2}\|\ddot{\varphi}_{i}\|_{L^{2}},\quad\|\dot{\varphi}_{i}\|_{L^{\infty}}\leq
ch^{1/2}\|\ddot{\varphi}_{i}\|_{L^{2}},\quad\|\dot{\varphi}_{i}\|_{L^{2}}\leq
ch\|\ddot{\varphi}_{i}\|_{L^{2}}$

(as follows from straightforward scaling arguments). Summarizing, there is a
constant $\kappa>0$ depending only on ${\mathcal{M}}$ and $\gamma$ such that

$\|U_{i}-\ddot{\varphi}_{i}\|_{L^{2}}\leq\kappa
h\|\ddot{\varphi}_{i}\|_{L^{2}},\quad\|V_{i}^{1}+V_{i}^{2}\|_{L^{2}}\leq\kappa
h^{3/2}\|\ddot{\varphi}_{i}\|_{L^{2}}^{2}.$

Furthermore, for $h$ small enough we have

$\|D_{t}^{2}\gamma_{v}\|_{L^{2}((t_{i-1},t_{i}))}^{2}=\|[\partial_{2}\exp_{\gamma_{i}}v_{i}]\ddot{v}_{i}+[\partial_{2}^{2}\exp_{\gamma_{i}}v_{i}](\dot{v}_{i},\dot{v}_{i})\|_{L^{2}((t_{i-1},t_{i}))}^{2}\\\
\leq
2\|[\partial_{2}\exp_{\gamma_{i}}v_{i}]\ddot{v}_{i}\|_{L^{2}((t_{i-1},t_{i}))}^{2}+2\|[\partial_{2}^{2}\exp_{\gamma_{i}}v_{i}](\dot{v}_{i},\dot{v}_{i})\|_{L^{2}((t_{i-1},t_{i}))}^{2}\leq
4\|\ddot{v}_{i}\|_{L^{2}((t_{i-1},t_{i}))}^{2}+1\leq
8\|D_{t}^{2}\gamma\|_{L^{2}}^{2}+9$

by lemma 6 and the previous estimates on $v_{i}$ so that

$\partial^{2}J[v_{1},\ldots,v_{N}](\phi,\phi)\geq
2\sum_{i=1}^{N}\left(\frac{1}{2}\|\ddot{\varphi}_{i}\|_{L^{2}}^{2}-\kappa^{2}h^{2}\|\ddot{\varphi}_{i}\|_{L^{2}}^{2}-\|D_{t}^{2}\gamma_{v}\|_{L^{2}((t_{i-1},t_{i}))}\kappa
h^{\frac{3}{2}}\|\ddot{\varphi}_{i}\|_{L^{2}}^{2}\right)\\\
\geq(1-2\kappa^{2}h^{2}-\kappa
h^{\frac{3}{2}}(8\|D_{t}^{2}\gamma\|_{L^{2}}+9))\sum_{i=1}^{N}\|\ddot{\varphi}_{i}\|_{L^{2}}^{2}>0$

for $h$ small enough (depending just on ${\mathcal{M}}$ and $\gamma$).
Therefore, $J$ is strictly convex on $\mathcal{A}_{h}$. ∎

We close the section by deriving the Euler–Lagrange equations satisfied by the
linear and cubic spline interpolations.

######  Proposition 10 (Optimality conditions).

  1. 1.

A linear spline interpolation $\gamma_{h}$ satisfies the optimality conditions

$\displaystyle 0=D_{t}^{2}\gamma_{h}\text{ on
}(t_{k-1},t_{k}),\,k=1,\ldots,N,$ $\displaystyle\gamma_{h}=\gamma\text{ at
}t_{0},\ldots,t_{N}.$

  2. 2.

A cubic spline interpolation $\gamma_{h}$ satisfies the optimality conditions

$\displaystyle
0=D_{t}^{4}\gamma_{h}+{\mathrm{R}}(D_{t}^{2}\gamma_{h},\dot{\gamma}_{h})\dot{\gamma}_{h}\text{
on }(t_{k-1},t_{k}),\,k=1,\ldots,N,$ $\displaystyle D_{t}\gamma_{h}\text{ and
}D_{t}^{2}\gamma_{h}\text{ exist and }\gamma_{h}=\gamma\text{ at
}t_{0},\ldots,t_{N},$ $\displaystyle D_{t}\gamma_{h}=D_{t}\gamma\text{ at
}t_{0}\text{ and }t_{N}.$

###### Proof.

  1. 1.

Extend $\gamma_{h}$ to a family
$\gamma_{h}:\mathbb{R}\times[0,1]\to{\mathcal{M}}$ of competing curves with
$\gamma(0,\cdot)$ being the minimizer. Then at $s=0$ we necessarily have

$0=\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}s}\int_{0}^{1}|\partial_{t}\gamma_{h}|^{2}\,\mathrm{d}t=\int_{0}^{1}(\partial_{t}\gamma_{h},D_{s}\partial_{t}\gamma_{h})\,\mathrm{d}t=\int_{0}^{1}(\partial_{t}\gamma_{h},D_{t}\partial_{s}\gamma_{h})\,\mathrm{d}t\\\
=\sum_{k=1}^{N}\left[(\partial_{t}\gamma_{h},\partial_{s}\gamma_{h})\right]_{t_{k-1}}^{t_{k}}-\int_{t_{k-1}}^{t_{k}}(D_{t}^{2}\gamma_{h},\partial_{s}\gamma_{h})\,\mathrm{d}t=-\sum_{k=1}^{N}\int_{t_{k-1}}^{t_{k}}(D_{t}^{2}\gamma_{h},\partial_{s}\gamma_{h})\,\mathrm{d}t$

where we used lemma 5, performed integration by parts using
$\frac{\mathrm{d}}{\mathrm{d}t}(v,w)=(D_{t}v,w)+(v,D_{t}w)$ and where we
exploited $\partial_{s}\gamma_{h}(t_{k})=0$ due to the admissibility of the
curve family $\gamma_{h}$. Since $\gamma_{h}$ was extended arbitrarily, the
result now follows from the fundamental lemma of the calculus of variations.

  2. 2.

Extend $\gamma_{h}$ to a family
$\gamma_{h}:\mathbb{R}\times[0,1]\to{\mathcal{M}}$ of competing curves with
$\gamma(0,\cdot)$ being the minimizer. Then at $s=0$ we necessarily have

$\displaystyle 0$
$\displaystyle=\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}s}\int_{0}^{1}|D_{t}^{2}\gamma_{h}|^{2}\,\mathrm{d}t$
$\displaystyle=\int_{0}^{1}(D_{t}^{2}\gamma_{h},D_{s}D_{t}^{2}\gamma_{h})\,\mathrm{d}t$
$\displaystyle=\int_{0}^{1}(D_{t}^{2}\gamma_{h},D_{t}D_{s}\partial_{t}\gamma_{h}+{\mathrm{R}}(\partial_{s}\gamma_{h},\partial_{t}\gamma_{h})\partial_{t}\gamma_{h})\,\mathrm{d}t$
$\displaystyle=\sum_{k=1}^{N}\left[(D_{t}^{2}\gamma_{h},D_{s}\partial_{t}\gamma_{h})\right]_{t_{k-1}}^{t_{k}}+\int_{t_{k-1}}^{t_{k}}(D_{t}^{2}\gamma_{h},{\mathrm{R}}(\partial_{s}\gamma_{h},\partial_{t}\gamma_{h})\partial_{t}\gamma_{h})-(D_{t}^{3}\gamma_{h},D_{s}\partial_{t}\gamma_{h})\,\mathrm{d}t$
$\displaystyle=\sum_{k=1}^{N}\left[(D_{t}^{2}\gamma_{h},D_{t}\partial_{s}\gamma_{h})\right]_{t_{k-1}}^{t_{k}}+\int_{t_{k-1}}^{t_{k}}(\partial_{s}\gamma_{h},{\mathrm{R}}(D_{t}^{2}\gamma_{h},\partial_{t}\gamma_{h})\partial_{t}\gamma_{h})-(D_{t}^{3}\gamma_{h},D_{t}\partial_{s}\gamma_{h})\,\mathrm{d}t$
$\displaystyle=\sum_{k=1}^{N}\left[(D_{t}^{2}\gamma_{h},D_{t}\partial_{s}\gamma_{h})\right]_{t_{k-1}}^{t_{k}}-\left[(D_{t}^{3}\gamma_{h},\partial_{s}\gamma_{h})\right]_{t_{k-1}}^{t_{k}}+\int_{t_{k-1}}^{t_{k}}(D_{t}^{4}\gamma_{h}+{\mathrm{R}}(D_{t}^{2}\gamma_{h},\partial_{t}\gamma_{h})\partial_{t}\gamma_{h},\partial_{s}\gamma_{h})\,\mathrm{d}t$
$\displaystyle=\sum_{k=1}^{N-1}((D_{t}^{-})^{2}\gamma_{h}(t_{k})-(D_{t}^{+})^{2}\gamma_{h}(t_{k}),D_{t}\partial_{s}\gamma_{h}(t_{k}))+\sum_{k=1}^{N}\int_{t_{k-1}}^{t_{k}}(D_{t}^{4}\gamma_{h}+{\mathrm{R}}(D_{t}^{2}\gamma_{h},\partial_{t}\gamma_{h})\partial_{t}\gamma_{h},\partial_{s}\gamma_{h})\,\mathrm{d}t,$

where $D_{t}^{+}$ and $D_{t}^{-}$ denote the left and right derivative. Above,
we again repeatedly exploited lemma 5 and performed integration by parts.
Furthermore we exploited $\partial_{s}\gamma_{h}(t_{k})=0$ and
$D_{t}\partial_{s}\gamma_{h}(t)=D_{s}\partial_{t}\gamma_{h}(t)=0$ at
$t=t_{0},t_{N}$ due to the admissibility of the curve family $\gamma_{h}$.
Since $\gamma_{h}$ was extended arbitrarily, the result now follows from the
fundamental lemma of the calculus of variations. ∎

##  4 Two short convergence proofs in the Euclidean case

Linear and cubic spline interpolation of curves in Euclidean space converge
quadratically and quartically, repectively. As our aim is to transfer this
result to the Riemannian setting, we briefly provide the corresponding
Euclidean proof here. We first give a simple proof exploiting the piecewise
polynomial structure of splines, in particular that the second derivative of
the cubic spline interpolation $\gamma_{h}$ is the $L^{2}$-best piecewise
linear approximation of $\ddot{\gamma}$ (in the remainder of the article we
will try to mimick this idea). Afterwards we give a second, more direct but
less elegant proof for cubic spline interpolation which essentially just
checks the stability of interpolation as well as its consistency (by comparing
the Taylor expansions of spline interpolant and original curve). In principle,
an analogous approach should also be possible on a manifold, but the
corresponding stability and Taylor estimates would best be obtained in a chart
and promise to be quite cumbersome. Below we will denote the Sobolev space of
integrable functions with essentially bounded derivatives up to fourth order
by $W^{4,\infty}$.

######  Theorem 11 (Spline convergence).

Let $\gamma\in W^{4,\infty}((0,1))$. Its linear spline interpolation
$\gamma_{h}$ satisfies

$\|\gamma_{h}-\gamma\|_{L^{p}}\leq ch^{2}$

for some $c>0$ depending only on $\gamma$ and $p$. The corresponding cubic
spline interpolation $\gamma_{h}$ satisfies

$\|\gamma_{h}-\gamma\|_{L^{p}}\leq ch^{4}$

for some $c>0$ depending only on $\gamma$ and $p$.

###### Proof.

For a function $f:[0,1]\to\mathbb{R}^{d}$ with zeros at $t_{0},\ldots,t_{N}$
one has $\|f\|_{L^{p}}\leq Ch^{2}\|\ddot{f}\|_{L^{p}}$ for some $C>0$
depending on $p$. Picking $f=\gamma_{h}-\gamma$ yields for the linear spline
interpolation (which satisfies $(\ddot{\gamma}_{h}=0$)

$\|\gamma_{h}-\gamma\|_{L^{p}}\leq Ch^{2}\|\ddot{\gamma}\|_{L^{p}}.$

For the cubic spline interpolation one obtains

$\|\gamma_{h}-\gamma\|_{L^{p}}\leq
Ch^{2}\|\ddot{\gamma}_{h}-\ddot{\gamma}\|_{L^{p}}.$

Now it turns out that $\ddot{\gamma}_{h}$ is the piecewise linear $L^{2}$-best
approximation of $\ddot{\gamma}$. Indeed, for any piecewise linear
$w:[0,1]\to\mathbb{R}^{d}$ one obtains the Galerkin orthogonality

$\int_{0}^{1}(w,\ddot{\gamma}_{h}-\ddot{\gamma})\,\mathrm{d}t=\sum_{i=1}^{N}\left[(w,\dot{\gamma}_{h}-\dot{\gamma})\right]_{t=t_{i-1}}^{t_{i}}-\int_{t_{i-1}}^{t_{i}}(\dot{w},\dot{\gamma}_{h}-\dot{\gamma})\,\mathrm{d}t=\left[(w,\dot{\gamma}_{h}-\dot{\gamma})-(\dot{w},\gamma_{h}-\gamma)\right]_{t=t_{i-1}}^{t_{i}}=0$

from two integrations by parts, exploiting that $\gamma_{h}-\gamma=0$ at the
interpolation points, that $\dot{\gamma}_{h}$ exists by proposition 10, and
that $\dot{\gamma}_{h}-\dot{\gamma}=0$ at $t=0$ and $t=1$. This immediately
implies (e. g. via Céa’s lemma) that $\ddot{\gamma}_{h}$ is the
$L^{2}$-projection or -best approximation of $\ddot{\gamma}$. Thus, by [4,
Lemma 4.3] (or by the $L^{\infty}$-boundedness of the $L^{2}$-projection onto
splines [3] and the resulting best approximation argument as in [1, Thm.
3.3.7]) we have

$\|\ddot{\gamma}_{h}-\ddot{\gamma}\|_{L^{\infty}}\leq
Ch^{2}\|\gamma\|_{W^{4,\infty}}$

for some $C>0$. ∎

The alternative, more direct proof for cubic spline interpolation just uses
stability and fourth order consistency of cubic spline interpolation via
Taylor expansion.

###### Proof.

The piecewise cubic spline interpolation can be written as a linear
combination of shifted and scaled versions of the cubic B-spline basis
function

$B(t)=\begin{cases}(3|t|^{3}-6t^{2}+4)/6&\text{if }|t|\leq 1,\\\
(2-|t|)^{3}/6&\text{if }|t|\in(1,2],\\\ 0&\text{else.}\end{cases}$

Writing $\gamma_{h}$ as a linear combination

$\gamma_{h}(t)=\sum_{i=-1}^{N+1}x_{i}B(\tfrac{t}{h}-i)$

of these basis functions, the interpolation constraints translate into the
linear system

$Ax=b\quad\text{with }A=\frac{1}{6}\begin{pmatrix}\frac{1}{2}&1\\\ 1&4&1\\\
&1&4&1\\\ &&\ddots&\ddots&\ddots\\\ &&&1&4&1\\\
&&&&1&\frac{1}{2}\end{pmatrix},\,x=\begin{pmatrix}x_{-1}\\\ \vdots\\\
x_{N+1}\end{pmatrix},\,b=\begin{pmatrix}\tfrac{\gamma(t_{0})}{4}-\tfrac{\dot{\gamma}(0)h}{12}\\\
\gamma(t_{0})\\\ \vdots\\\ \gamma(t_{N})\\\
\tfrac{\gamma(t_{N})}{4}+\tfrac{\dot{\gamma}(1)h}{12}\end{pmatrix}.$

Now for a point $t\in[t_{i},t_{i+1}]$ we abbreviate $a_{i}=B(\tfrac{t}{h}-i)$
(thus all $a_{i}$ are zero except for $a_{i-1},\ldots,a_{i+2}$) so that

$\gamma_{h}(t)-\gamma(t)=a^{T}x-\gamma(t)=a^{T}A^{-1}b-\gamma(t)\quad\text{with
}a=(a_{-1},\ldots,a_{N+1})^{T}.$

Now define
$p(r)=\gamma(t)+\dot{\gamma}(t)(r-t)+\frac{1}{2}\ddot{\gamma}(t)(r-t)^{2}+\frac{1}{6}\dddot{\gamma}(t)(r-t)^{3}$
to be the third order Taylor expansion of $\gamma$ and notice

$b=b_{1}+b_{2}=\begin{pmatrix}p(t_{0})/4-\dot{p}(0)h/12\\\ p(t_{0})\\\
\vdots\\\ p(t_{N})\\\
p(t_{N})/4+\dot{p}(1)h/12\end{pmatrix}+\begin{pmatrix}O(|0-t|^{4})\\\
O(|t_{0}-t|^{4})\\\ \vdots\\\ O(|t_{N}-t|^{4})\\\ O(|1-t|^{4})\end{pmatrix}.$

Since $a^{T}A^{-1}b_{1}=p(t)=\gamma(t)$, we arrive at

$|\gamma_{h}(t)-\gamma(t)|=|a^{T}A^{-1}b_{2}|=\left|\sum_{j=i-1}^{i+2}a_{j}A^{-1}_{j,:}b_{2}\right|\leq\sum_{j=i-1}^{i+2}|A^{-1}_{j,:}b_{2}|,$

where $A^{-1}_{j,:}=(q_{-1},q_{0},\ldots,q_{N+1})$ denotes the $j$th row of
$A^{-1}$ and thus necessarily satisfies

$\displaystyle q_{-1}=-2q_{0}\quad\text{and}\quad
q_{k}=-(q_{k-2}+4q_{k-1}),\,k=1,\ldots,j,$ $\displaystyle
q_{N+1}=-2q_{N}\quad\text{and}\quad
q_{k}=-(q_{k+2}+4q_{k+1}),\,k=N-1,\ldots,j.$

By induction this immediately implies

$q_{k}q_{k-1}<0,\,k=0,\ldots,N+1,\qquad|q_{k}|\begin{cases}\geq|q_{k-1}|,\,k=1,\ldots,j,\\\
\leq|q_{k-1}|,\,k=j+1,\ldots,N,\end{cases}\quad\text{ and }|q_{k}|\leq
3^{-|k-j|}|q_{j}|,\,k=0,\ldots,N.$

Furthermore, the equation $1=q_{j-1}+4q_{j}+q_{j+1}$ together with the fact
that $q_{j-1},q_{j+1}$ have a different sign than $q_{j}$ and are at least by
the factor $3$ smaller in absolute value implies

$|q_{j}|\leq 1\quad\text{and thus}\quad|q_{k}|\leq 3^{-|k-j|},\,k=0,\ldots,N.$

Summarizing, we obtain

$|\gamma_{h}(t)-\gamma(t)|\leq\sum_{j=i-1}^{i+2}|A^{-1}_{j,:}b_{2}|\leq\sum_{j=i-1}^{i+2}\left|\sum_{k=-1}^{N+1}3^{-|k-j|}O(|t_{k}-t|^{4})\right|\leq
4h^{4}\sum_{k=-\infty}^{\infty}3^{-|k-j|}(|k-j|+1)^{4}\leq Ch^{4}$

for a constant depending only on the maximum fourth derivative of $\gamma$. ∎

######  Corollary 12 (Perturbed spline).

Let $g:[0,1]\to\mathbb{R}^{d}$ be twice differentiable with
$g(t_{i})=\dot{g}(0)=\dot{g}(1)=0$ for $i=0,\ldots,N$ and $g\in
W^{4,\infty}((t_{i-1},t_{i}))$ with
$\|\partial_{t}^{4}g\|_{L^{\infty}((t_{i-1},t_{i}))}\leq C$ for
$i=1,\ldots,N$. Then $\|g\|_{L^{\infty}}\leq ch^{4}$ with a constant $c>0$
depending only on $C$.

###### Proof.

Define $\gamma:[0,1]\to\mathbb{R}^{d}$ to be the curve in
$W^{4,\infty}((0,1))$ with
$\gamma(0)=\dot{\gamma}(0)=\gamma(1)=\dot{\gamma}(1)=0$ whose fourth weak
derivative is $f\in L^{\infty}((0,1))$ with $f=\partial_{t}^{4}g$ in each
$(t_{i-1},t_{i})$, that is, we set

$\displaystyle\tilde{\gamma}(t)$
$\displaystyle=\int_{0}^{t}\int_{0}^{s}\int_{0}^{r}\int_{0}^{q}f(p)\,\mathrm{d}p\,\mathrm{d}q\,\mathrm{d}r\,\mathrm{d}s,$
$\displaystyle\gamma(t)$
$\displaystyle=\tilde{\gamma}(t)+(\dot{\tilde{\gamma}}(1)-3\tilde{\gamma}(1))t^{2}+(2\tilde{\gamma}(1)-\dot{\tilde{\gamma}}(1))t^{3}.$

Letting now $\gamma_{h}$ be the cubic spline interpolation of $\gamma$, we
obviously have $\gamma-\gamma_{h}=g$ and thus $\|g\|_{L^{\infty}}\leq ch^{4}$,
where the constant only depends on the maximum fourth derivative of $\gamma$
(which is bounded by $C$). ∎

##  5 The approximation properties of Riemannian cubic splines

In this section we transfer the proof of theorem 11 for the Euclidean case to
the Riemannian setting. Essentially, we will proceed in three steps: we will
first derive a generalization of the estimate
$\|\gamma_{h}-\gamma\|_{L^{p}}\leq
Ch^{2}\|\ddot{\gamma}_{h}-\ddot{\gamma}\|_{L^{p}}$, then we will introduce a
linear interpolation operator ${\mathcal{I}^{1}_{0}}$ of vector fields as a
preparation for estimating $\|\ddot{\gamma}_{h}-\ddot{\gamma}\|_{L^{p}}$
(which we think of as some kind of best approximation error) by
$\|\ddot{\gamma}_{h}-{\mathcal{I}^{1}_{0}}\ddot{\gamma}\|_{L^{p}}$, and
finally we will derive boundedness of $\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}$
using the $L^{2}$-best approximation idea from theorem 11 and bootstrap this
$L^{2}$-best approximation argument with the help of the bound on
$\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}$ to obtain a generalization of the
Euclidean estimate $\|\ddot{\gamma}_{h}-\ddot{\gamma}\|_{L^{\infty}}\leq
Ch^{2}\|\gamma\|_{W^{4,\infty}}$.

###  5.1 Quadratic convergence rate exploiting local interpolation

Just like in the proof of theorem 11 for the Euclidean case we first derive a
quadratic dependence of the spline interpolation error on $h$.

######  Lemma 13 (Approximation lemma).

Let $\gamma:[0,1]\to{\mathcal{M}}$ be a curve and
$v\in\gamma^{\\#}T{\mathcal{M}}$ a vector field along the curve with
$v(t_{k})=v(t_{k+1})=\ldots=v(t_{k+l})=0$, then for $p\in[1,\infty]$ and
$I=(t_{k},t_{k+l})$ we have

$\|v\|_{L^{p}(I)}\leq ch^{l+1}\|D_{t}^{l+1}v\|_{L^{p}(I)}$

for a constant $c$ depending only on $l$ and $p$.

###### Proof.

Let $\pi_{\gamma}:T_{\gamma(t)}{\mathcal{M}}\to T_{\gamma(0)}{\mathcal{M}}$
denote the parallel transport along the curve $\gamma$. Then
$\pi_{\gamma}v:[0,1]\to T_{\gamma(0)}{\mathcal{M}}$ satisfies
$\pi_{\gamma}v(t_{k})=\ldots=\pi_{\gamma}v(t_{k+l})=0$ so that the standard
approximation result for curves in Euclidean space [1, Thm. 3.1.5] or a
straightforward scaling argument can be applied, giving

$\|v\|_{L^{p}(I)}=\|\pi_{\gamma}v\|_{L^{p}(I)}\leq
ch^{l+1}\|\partial_{t}^{l+1}(\pi_{\gamma}v)\|_{L^{p}(I)}=ch^{l+1}\|D_{t}^{l+1}v\|_{L^{p}(I)}$

by lemma 5(1). ∎

This immediately implies the generalization of the Euclidean estimate
$\|\gamma_{h}-\gamma\|_{L^{p}}\leq
Ch^{2}\|\ddot{\gamma}_{h}-\ddot{\gamma}\|_{L^{p}}$.

######  Proposition 14 (Quadratic convergence).

For any $p\in[1,\infty]$ the linear and cubic spline interpolation of an
absolutely continuous curve $\gamma:[0,1]\to{\mathcal{M}}$ satisfy

$\|d(\gamma_{h},\gamma)\|_{L^{p}(I)}\leq
ch^{2}\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{p}(I)}$

on each interval $I=(t_{i-1},t_{i})$ for a constant $c>0$ independent of
$\gamma$, $\gamma_{h}$, or $h$, for all $h$ small enough (both depending only
on ${\mathcal{M}}$ and $\gamma$).

###### Proof.

Simply apply lemma 13 with $l=1$ to $v=\log_{\gamma_{h}}\\!\gamma$ and note
$d(\gamma_{h},\gamma)=|\log_{\gamma_{h}}\\!\gamma|_{\gamma_{h}}$ as long as
$h$ is small enough for $v$ to be well-defined (follows from proposition 8 and
continuity of $\gamma$). ∎

The full convergence result of theorems 2 and 4 for linear and cubic spline
interpolation will follow from estimating the right-hand side of the previous
inequality. To this end we expand

$D_{t}^{2}\log_{\gamma_{h}}\\!\gamma=\partial_{2}\log_{\gamma_{h}}\\!\gamma\,(D_{t}^{2}\gamma)+\partial_{1}\log_{\gamma_{h}}\\!\gamma\,(D_{t}^{2}\gamma_{h})+\partial_{2}^{2}\log_{\gamma_{h}}\\!\gamma\,(\dot{\gamma},\dot{\gamma})+\partial_{2}\partial_{1}\log_{\gamma_{h}}\\!\gamma\,(\dot{\gamma},\dot{\gamma}_{h})\\\
+\partial_{1}\partial_{2}\log_{\gamma_{h}}\\!\gamma\,(\dot{\gamma}_{h},\dot{\gamma})+\partial_{1}^{2}\log_{\gamma_{h}}\\!\gamma\,(\dot{\gamma}_{h},\dot{\gamma}_{h}).$

Note that at this point, proposition 14 can already be turned into a quadratic
convergence result (which for linear splines already is the optimal rate),
which is briefly done in the remainder of this paragraph. Indeed, with lemma 6
and the boundedness of $\gamma$ and $\gamma_{h}$ from proposition 8 or more
specifically from (1)-(3) we can estimate

$\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{p}(I)}\leq
C\left(\|D_{t}^{2}\gamma_{h}\|_{L^{p}(I)}+\|D_{t}^{2}\gamma\|_{L^{p}(I)}+\|(|\dot{\gamma}|+|\dot{\gamma}_{h}|)^{2}\|_{L^{p}(I))}\right)$

for some $C>0$ if $h$ is small enough (depending on $\gamma$ and
${\mathcal{M}}$), which leads to the following result.

######  Corollary 15 (Quadratic $L^{1}$\- and $L^{2}$-convergence).

Let $\gamma:[0,1]\to{\mathcal{M}}$ be twice differentiable, then its linear
spline interpolation $\gamma_{h}$ satisfies

$\|d(\gamma_{h},\gamma)\|_{L^{1}}\leq ch^{2}$

for some constant $c>0$ and all $h$ small enough (both depending only on
${\mathcal{M}}$ and $\gamma$). If $\gamma$ is four times differentiable, then
its cubic spline interpolation $\gamma_{h}$ satisfies

$\|d(\gamma_{h},\gamma)\|_{L^{2}}\leq ch^{2}$

for some constant $c>0$ and all $h$ small enough (again both depending only on
${\mathcal{M}}$ and $\gamma$).

###### Proof.

By proposition 14, for the linear spline we have

$\|d(\gamma_{h},\gamma)\|_{L^{1}}\leq
ch^{2}\sum_{i=1}^{N}\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{1}((t_{i-1},t_{i}))}\leq
ch^{2}\sum_{i=1}^{N}C\left(\|D_{t}^{2}\gamma\|_{L^{1}((t_{i-1},t_{i}))}+\||\dot{\gamma}|+|\dot{\gamma}_{h}|\|_{L^{2}((t_{i-1},t_{i}))}^{2}\right)\\\
=cCh^{2}\left(\|D_{t}^{2}\gamma\|_{L^{1}}+\||\dot{\gamma}|+|\dot{\gamma}_{h}|\|_{L^{2}}^{2}\right)\leq
cCh^{2}\left(\|D_{t}^{2}\gamma\|_{L^{1}}+2\|\dot{\gamma}\|_{L^{2}}^{2}+2\|\dot{\gamma}_{h}\|_{L^{2}}^{2}\right)\leq
cCh^{2}\left(\|D_{t}^{2}\gamma\|_{L^{1}}+4\|\dot{\gamma}\|_{L^{2}}^{2}\right),$

where we exploited the optimality condition $D_{t}^{2}\gamma_{h}=0$ as well as
$\|\dot{\gamma}_{h}\|_{L^{2}}^{2}\leq\|\dot{\gamma}\|_{L^{2}}^{2}$. Similarly,
for the cubic spline we have

$\|d(\gamma_{h},\gamma)\|_{L^{2}}\leq
ch^{2}\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{2}}\leq
cCh^{2}\left(\|D_{t}^{2}\gamma\|_{L^{2}}+\|D_{t}^{2}\gamma_{h}\|_{L^{2}}+\|(|\dot{\gamma}|+|\dot{\gamma}_{h}|)^{2}\|_{L^{2}}\right)\\\
\leq
cCh^{2}\left(2\|D_{t}^{2}\gamma\|_{L^{2}}+(\|\dot{\gamma}\|_{L^{\infty}}+\|\dot{\gamma}_{h}\|_{L^{\infty}})^{2}\right)\leq
cCh^{2}\left(2\|D_{t}^{2}\gamma\|_{L^{2}}+(\|\dot{\gamma}\|_{L^{\infty}}+|\dot{\gamma}(0)|+2\|D_{t}^{2}\gamma\|_{L^{2}})^{2}\right),$

where we exploited
$\|D_{t}^{2}\gamma_{h}\|_{L^{2}}\leq\|D_{t}^{2}\gamma\|_{L^{2}}$ as well as
(2). ∎

To improve the previous convergence estimate for linear spline interpolation
to the final estimate of theorem 4 we require the following uniform bound on
$\dot{\gamma}_{h}$, which immediately follows from properties of unit speed
parameterized geodesics.

######  Proposition 16 (Boundedness of first derivative).

The linear spline interpolation $\gamma_{h}$ of a curve $\gamma$ satisfies

$\|\dot{\gamma}_{h}\|_{L^{\infty}}\leq\|\dot{\gamma}\|_{L^{\infty}}.$

###### Proof.

Recall that between $t_{i-1}$ and $t_{i}$ the curve $\gamma_{h}$ is a shortest
unit speed parameterized geodesic between $\gamma(t_{i-1})$ and
$\gamma(t_{i})$, thus on $(t_{i-1},t_{i})$ we have

$|\dot{\gamma}_{h}|=d(\gamma(t_{i-1}),\gamma(t_{i}))/h\leq\int_{t_{i-1}}^{t_{i}}|\dot{\gamma}|\,\mathrm{d}t/h\leq\|\dot{\gamma}\|_{L^{\infty}}.\qed$

Together with proposition 16 for linear spline interpolation and
$\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}\leq
3\|D_{t}^{2}\gamma\|_{L^{\infty}}+Ch$ for cubic spline interpolation, derived
later in corollary 24, which by a bootstrapping argument makes use of
corollary 15, one finally arrives at the following quadratic convergence
estimate.

######  Corollary 17 (Quadratic $L^{\infty}$-convergence).

The linear and cubic spline interpolation satisfy

$\|d(\gamma_{h},\gamma)\|_{L^{\infty}}\leq ch^{2}$

for some constant $c>0$ depending only on $\|\dot{\gamma}\|_{L^{\infty}}$,
$\|D_{t}^{2}\gamma\|_{L^{\infty}}$, and ${\mathcal{M}}$, for $h$ small enough
depending on ${\mathcal{M}}$ and $\gamma$.

###### Proof.

Abbreviating $I_{i}=(t_{i-1},t_{i})$, by proposition 14 and proposition 16 or
corollary 24 we have

$\displaystyle\|d(\gamma_{h},\gamma)\|_{L^{\infty}}$ $\displaystyle\leq
ch^{2}\max_{i=1,\ldots,N}\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{\infty}(I_{i})}$
$\displaystyle\leq
ch^{2}\max_{i=1,\ldots,N}C\left(\|D_{t}^{2}\gamma\|_{L^{\infty}(I_{i})}+\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}(I_{i})}+\||\dot{\gamma}|+|\dot{\gamma}_{h}|\|_{L^{\infty}(I_{i})}^{2}\right)$
$\displaystyle\leq\tilde{C}h^{2}\left(\|D_{t}^{2}\gamma\|_{L^{\infty}}+Ch+\|\dot{\gamma}\|_{L^{\infty}}^{2}\right),$

for some $\tilde{C}>0$, where in the case of linear splines we exploited that
$D_{t}^{2}\gamma_{h}=0$. ∎

###  5.2 A linear interpolation operator of vector fields

As mentioned before, the guiding idea of our approach is to interpret
$D_{t}^{2}\gamma_{h}$ as an approximation to $D_{t}^{2}\gamma$ which is (maybe
up to a constant factor) as good as the $L^{2}$-best approximation or a
piecewise linear interpolation. To formalize this, below we introduce a
piecewise linear interpolation operator of vector fields along curves and
analyse its well-posedness and approximation properties.

######  Definition 18 (Linear interpolation operator).

Let $\gamma:[0,1]\to{\mathcal{M}}$ be a piecewise differentiable and globally
continuous curve and $v\in\gamma^{\\#}T{\mathcal{M}}$ a continuous vector
field along $\gamma$. We define the piecewise linear interpolation
${\mathcal{I}^{1}_{0}}v$ of $v$ as solution $w$ to the linear boundary value
problem

$\displaystyle D_{t}^{2}w+{\mathrm{R}}(w,\dot{\gamma})\dot{\gamma}=0\text{ on
}(t_{k-1},t_{k}),\,k=1,\ldots,N,$ (4) $\displaystyle
w(t_{k})=v(t_{k}),\,k=0,\ldots,N.$

######  Lemma 19 (Well-posedness and boundedness of linear interpolation).

For $h$ small enough depending on ${\mathcal{M}}$ and
$\|\dot{\gamma}\|_{L^{\infty}}$, ${\mathcal{I}^{1}_{0}}$ is well-defined, that
is, the boundary value problem (4) has a unique solution. Furthermore, for
every continuous vector field $v\in\gamma^{\\#}T{\mathcal{M}}$ we have

$\displaystyle\|{\mathcal{I}^{1}_{0}}v\|_{L^{\infty}}$ $\displaystyle\leq
2\max_{k=0,\ldots,N}|v(t_{k})|\leq 2\|v\|_{L^{\infty}}.$

###### Proof.

As the boundary value problem (4) is a linear system of second order ordinary
differential equations, we obtain existence of solutions by standard ordinary
differential equation theory. Uniqueness and the estimate on the norm follow
similarly to Rauch comparison estimates for classical Jacobi fields on each
interval $(t_{k-1},t_{k})$. In detail, assume that the scalar curvature of
${\mathcal{M}}$ is bounded from above by $\mu\geq 0$, and let $h$ be small
enough such that $h\|\dot{\gamma}\|_{L^{\infty}}^{2}\leq\frac{\pi}{2\mu}$. Let
$w$ be defined on $(t_{k-1},t_{k})$ by (4) with $w(t_{k-1})=0$ and
$w(t_{k})=\hat{v}\in T_{\gamma(t_{k})}{\mathcal{M}}$. Set

$\displaystyle
f(t)=\partial_{t}|w|(t_{k-1})\;\tfrac{2h}{\pi}\sin\left(\tfrac{\pi(t-t_{k-1})}{2h}\right)\quad\text{for
}t\in[t_{k-1},t_{k}],$

then

$\displaystyle\partial_{t}^{2}f+\tfrac{\pi}{2h}f=0\text{ on
}(t_{k-1},t_{k})\quad\text{with}$ $\displaystyle
f(t_{k-1})=0,\qquad\partial_{t}f(t_{k-1})=\partial_{t}|w|(t_{k-1}),\quad\text{and
}0<f(t)<f(t_{k})\quad\text{for }t\in(t_{k-1},t_{k}).$

As long as $w(t)\neq 0$, abbreviating by
$K(w,\dot{\gamma})=(w,{\mathrm{R}}(w,\dot{\gamma})\dot{\gamma})/\left(|w|^{2}|\dot{\gamma}|^{2}-(w,\dot{\gamma})^{2}\right)$
the sectional curvature, we have

$\displaystyle\partial_{t}^{2}|w|$
$\displaystyle=-\frac{1}{|w|}K(w,\dot{\gamma})\left(|w|^{2}|\dot{\gamma}|^{2}-(w,\dot{\gamma})^{2}\right)+\frac{1}{|w|^{3}}\left(|w|^{2}|D_{t}w|^{2}-(w,D_{t}w)^{2}\right)\geq-\mu|\dot{\gamma}|^{2}|w|$

and thus

$\displaystyle\partial_{t}\left(f\partial_{t}|w|-|w|\partial_{t}f\right)=f\partial_{t}^{2}|w|-|w|\partial_{t}^{2}f\geq
0.$

As $(f\partial_{t}|w|-|w|\partial_{t}f)(t_{k-1})=0$, this implies

$\displaystyle\partial_{t}\left(\tfrac{|w|}{f}\right)=\frac{1}{f^{2}}\left(f\partial_{t}|w|-|w|\partial_{t}f\right)\geq
0$

as long as $w(t)\neq 0$, and thus for $s>t>t_{k-1}$,

$\displaystyle\frac{|w|(s)}{f(s)}\geq\frac{|w|(t)}{f(t)}\geq\lim_{r\searrow
t_{k-1}}\frac{|w|(r)}{f(r)}=1,$

where the limit follows by de l’Hôpital’s rule. This shows in particular, that
$|w|$ cannot have a zero before $f$ does, and thus the estimate

$\displaystyle|w|(t)\leq\frac{f(t)}{f(t_{k})}|w|(t_{k})\leq|\hat{v}|$

holds for all $t\in[t_{k-1},t_{k}]$. Together with the linearity of (4), this
stability result implies uniqueness of the solution to the boundary value
problem as well as (using that the same estimate is obtained when the roles of
$t_{k-1}$ and $t_{k}$ are swapped) the estimate

$\|{\mathcal{I}^{1}_{0}}v\|_{L^{\infty}}\leq
2\max_{k=0,\ldots,N}|v(t_{k})|\leq 2\|v\|_{L^{\infty}}.\qed$

######  Lemma 20 (Approximation properties of linear interpolation).

Let $h$ be small enough such that lemma 19 applies, and let $p\in[1,\infty]$.
There exists a constant $c>0$ depending on $p$ and ${\mathcal{M}}$ such that
for any vector field $v\in\gamma^{\\#}T{\mathcal{M}}$ the estimate

$\displaystyle\|D_{t}^{2}{\mathcal{I}^{1}_{0}}v\|_{L^{p}(I_{k})}$
$\displaystyle\leq
c\|\dot{\gamma}\|_{L^{\infty}(I_{k})}^{2}\|{\mathcal{I}^{1}_{0}}v\|_{L^{p}(I_{k})}$

holds on each interval $I_{k}=(t_{k},t_{k+1})$. For $p>q\geq 1$ and
$\alpha=(q^{-1}-p^{-1})$ there exist constants $c>0$ depending on $p$, $q$ and
${\mathcal{M}}$ such that the inverse estimates

$\displaystyle\|{\mathcal{I}^{1}_{0}}v\|_{L^{p}}$ $\displaystyle\leq
c\left(\|\dot{\gamma}\|^{\alpha}_{L^{\infty}}+h^{-\alpha}\right)\|{\mathcal{I}^{1}_{0}}v\|_{L^{q}},$
$\displaystyle\|D_{t}{\mathcal{I}^{1}_{0}}v\|_{L^{p}}$ $\displaystyle\leq
c\left(\|\dot{\gamma}\|_{L^{\infty}}^{1+\alpha}+h^{-1-\alpha}\right)\|{\mathcal{I}^{1}_{0}}v\|_{L^{q}}$

hold. If $\|D_{t}^{l}v\|_{L^{p}}$ is bounded for $l=0,1,2$, we further have
the interpolation error estimate

$\displaystyle\|v-{\mathcal{I}^{1}_{0}}v\|_{L^{p}}\leq
c\;h^{2}\left(\|D_{t}^{2}v\|_{L^{p}}+\|\dot{\gamma}\|_{L^{\infty}}^{2}\|v\|_{L^{\infty}}\right).$

###### Proof.

The estimate for $D_{t}^{2}{\mathcal{I}^{1}_{0}}v$ follows on each interval
$I_{k}$ directly from the definition of ${\mathcal{I}^{1}_{0}}$ with
$c=\|{\mathrm{R}}\|_{L^{\infty}}$. The inverse estimates follow from the
scaled Gagliardo–Nirenberg–Sobolev inequality and the estimates on the second
derivatives,

$\displaystyle\|D_{t}^{j}{\mathcal{I}^{1}_{0}}v\|_{L^{p}(I_{k})}$
$\displaystyle\leq
c\|D_{t}^{2}{\mathcal{I}^{1}_{0}}v\|_{L^{q}(I_{k})}^{\beta}\|{\mathcal{I}^{1}_{0}}v\|_{L^{q}(I_{k})}^{1-\beta}+ch^{-j+p^{-1}-q^{-1}}\|{\mathcal{I}^{1}_{0}}v\|_{L^{q}(I_{k})}$
$\displaystyle\leq
c\left(\|\dot{\gamma}\|_{L^{\infty}}^{2\beta}+h^{-j+p^{-1}-q^{-1}}\right)\|{\mathcal{I}^{1}_{0}}v\|_{L^{q}(I_{k})}$

with $p^{-1}=j-2\beta+q^{-1}$ and $j=0,1$. The global estimates are now
obtained by taking the $p$th power, summing over all subintervals $I_{k}$, and
exploiting that $\sum_{i}a_{i}^{p/q}\leq(\sum_{i}a_{i})^{p/q}$ for any vector
of nonnegative values $a_{i}$. The interpolation error follows from lemma 13
with $l=1$, the estimate for $D_{t}^{2}{\mathcal{I}^{1}_{0}}v$ on each
$I_{k}$, and lemma 19. ∎

######  Remark 21 (Vector field interpolation along splines).

We will always apply the vector field interpolation ${\mathcal{I}^{1}_{0}}$
along the cubic spline interpolations $\gamma_{h}$ of some continuous curve
$\gamma$. The dependence on $\gamma_{h}$ in the estimates of lemma 19 and
lemma 20 only comprises bounds on $\|\dot{\gamma}_{h}\|_{L^{\infty}}$ that can
be estimated in terms of $\gamma$ using proposition 8. Thus, by slight abuse
of notation, also for $v\in\gamma^{\\#}T{\mathcal{M}}$ we will take
${\mathcal{I}^{1}_{0}}v\in\gamma_{h}^{\\#}T{\mathcal{M}}$.

######  Remark 22 (Alternative interpolation operators).

Note that the linear vector field interpolation defined by the boundary value
problem (4) is natural along both linear and cubic splines in the following
sense:

  1. 1.

Let $\gamma_{h}(s,t)$ be a family of linear splines, then its variation is
linear in the sense
${\mathcal{I}^{1}_{0}}\partial_{s}\gamma_{h}=\partial_{s}\gamma_{h}$. Indeed,
the optimality conditions in proposition 10 imply
$0=D_{s}(D_{t}^{2}\gamma_{h})=D_{t}D_{s}\partial_{t}\gamma_{h}+{\mathrm{R}}(\partial_{s}\gamma_{h},\partial_{t}\gamma_{h})\partial_{t}\gamma_{h}=D_{t}^{2}\partial_{s}\gamma_{h}+{\mathrm{R}}(\partial_{s}\gamma_{h},\partial_{t}\gamma_{h})\partial_{t}\gamma_{h}$
(where we used lemma 5), the well-known equation of Jacobi fields.

  2. 2.

For cubic splines the second derivatives $D_{t}^{2}\gamma_{h}$ are linear in
the sense that ${\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma_{h}=D_{t}^{2}\gamma_{h}$.

It is obvious that there are multiple possible alternatives to generalize
piecewise linear interpolation from Euclidean space to vector fields on
Riemannian manifolds. In fact replacing (4) by $D_{t}^{2}w=0$ would also work
in the following convergence analysis with just slight modifications to the
proofs, however, it seems a little less natural. Indeed, the motivation for
definition 18 is to view a (piecewise) linear vector field as a second
derivative of a cubic spline (which in Euclidean space is a piecewise cubic
polynomial). Would one consider higher order spline interpolation, that is,
minimization of $\int_{0}^{1}|D_{t}^{k}\gamma_{h}|^{2}\,\mathrm{d}t$ under
interpolation constraints, then again a corresponding natural definition of
${\mathcal{I}^{1}_{0}}$ would have been the $(2k-2)$th derivative of local
minimizers of that energy (and the special case $k=2$ leads to definition 18).

###  5.3 Quartic convergence of cubic splines

We aim to derive a quartic convergence rate by applying corollary 12 to
$g=\pi_{\gamma_{h}}\log_{\gamma_{h}}\\!\gamma$. To this end we need to show
boundedness of $\partial_{t}^{4}g$, which can be reduced to the boundedness of
$D_{t}^{4}\gamma_{h}$ and thus, via proposition 10(2), of
$D_{t}^{2}\gamma_{h}$. The derivation of such a bound occupies the major part
of this paragraph, at the end of which the quartic convergence is deduced.
Exploiting ${\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma_{h}=D_{t}^{2}\gamma_{h}$, we
could bound $D_{t}^{2}\gamma_{h}$ by lemma 19 if we knew its values at the
interpolation times $t_{0},\ldots,t_{N}$, however, those are not available.
This is in contrast to the piecewise linear interpolation
${\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma$ of $D_{t}^{2}\gamma$ (along
$\gamma_{h}$), which thus can readily be bounded in $L^{\infty}((0,1))$.
Therefore we will show boundedness of $D_{t}^{2}\gamma_{h}$ by bounding its
difference to ${\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma$. The proof, in turn,
approximates this difference by $D_{t}^{2}\log_{\gamma_{h}}\\!\gamma$, the
generalization of $\ddot{\gamma}_{h}-\ddot{\gamma}$ from the Euclidean case,
for which there is a Galerkin orthogonality and consequently a best-
approximation result.

######  Proposition 23 (Approximation error for $D_{t}^{2}\gamma_{h}$).

For $\gamma$ four times differentiable and $h$ small enough (depending on
${\mathcal{M}}$ and $\gamma$) the cubic spline interpolation satisfies

$\displaystyle\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma-(D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma)\|_{L^{p}}$
$\displaystyle\leq
c\left(h^{2}+\|d(\gamma_{h},\gamma)\|_{L^{p}}\right)\left(\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}+1\right),$
$\displaystyle\|D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma\|_{L^{2}}$
$\displaystyle\leq
c\left(h^{2}+\|d(\gamma_{h},\gamma)\|_{L^{2}}\right)\left(\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}+1\right)$

for a constant $c$ only depending on ${\mathcal{M}}$ and $\gamma$.

###### Proof.

We will estimate the difference by adding a zero and using the triangle
inequality to split into easier terms. In particular, we will add

$0=-D_{t}^{2}\log_{\gamma_{h}}\\!\gamma+D_{t}(\partial_{2}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma}))+D_{t}(\partial_{1}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma}_{h}))$

to obtain

$\displaystyle D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma$
$\displaystyle=-D_{t}^{2}\log_{\gamma_{h}}\\!\gamma+\left[D_{t}(\dot{\gamma}_{h}+\partial_{1}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma}_{h}))\right]+\left[D_{t}(\partial_{2}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma}))-{\mathcal{I}^{1}_{0}}(\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma))\right]$
$\displaystyle=-D_{t}^{2}\log_{\gamma_{h}}\\!\gamma+\left[D_{t}^{2}\gamma_{h}+\partial_{1}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma_{h})\right]+\left[\partial_{1}^{2}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma}_{h},\dot{\gamma}_{h})+\partial_{2}\partial_{1}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma},\dot{\gamma}_{h})\right]$
$\displaystyle\quad+\left[\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma)-{\mathcal{I}^{1}_{0}}(\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma))\right]+\left[\partial_{1}\partial_{2}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma}_{h},\dot{\gamma})+\partial_{2}^{2}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma},\dot{\gamma})\right],$

where we additionally exploited
${\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma={\mathcal{I}^{1}_{0}}(\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma))$.
We now estimate the bracketed terms. By lemma 20, we have

$\displaystyle\|\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma)-{\mathcal{I}^{1}_{0}}(\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma))\|_{L^{p}}$
$\displaystyle\leq
c\;h^{2}\left(\|D_{t}^{2}\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma)\|_{L^{p}}+\|\dot{\gamma}\|_{L^{\infty}}^{2}\|\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma)\|_{L^{\infty}}\right).$

We have

$D_{t}^{2}\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma)=\partial_{1}^{2}\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma_{h},D_{t}^{2}\gamma)+r,$

where $r$ only contains derivatives up to third order of the logarithm
evaluated in directions of derivatives of $\gamma$ up to fourth order and of
$\gamma_{h}$ up to first order. Thus, $r$ can be bounded by a constant
(depending on ${\mathcal{M}}$ and $\gamma$). Together we obtain

$\|\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma)-{\mathcal{I}^{1}_{0}}(\partial_{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma))\|_{L^{p}}\leq
ch^{2}(\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}+1).$

By lemma 6 we furthermore have

$\displaystyle\left\|\left[\partial_{1}^{2}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma}_{h},\dot{\gamma}_{h})+\partial_{2}\partial_{1}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma},\dot{\gamma}_{h})\right]+\left[\partial_{1}\partial_{2}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma}_{h},\dot{\gamma})+\partial_{2}^{2}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma},\dot{\gamma})\right]\right\|_{L^{p}}$
$\displaystyle\hskip 216.81pt\leq
c\|d(\gamma_{h},\gamma)\|_{L^{p}}(\|\dot{\gamma}\|_{L^{\infty}}+\|\dot{\gamma}_{h}\|_{L^{\infty}})^{2}\leq\tilde{c}\|d(\gamma_{h},\gamma)\|_{L^{p}},$
$\displaystyle\left\|D_{t}^{2}\gamma_{h}+\partial_{1}\log_{\gamma_{h}}\\!\gamma(D_{t}^{2}\gamma_{h})\right\|_{L^{p}}\leq
c\|d(\gamma_{h},\gamma)\|_{L^{\infty}}\|d(\gamma_{h},\gamma)\|_{L^{p}}\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}\leq\tilde{c}\|d(\gamma_{h},\gamma)\|_{L^{p}}\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}$

so that in summary

$\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma-(D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma)\|_{L^{p}}\leq
c\left(h^{2}+\|d(\gamma_{h},\gamma)\|_{L^{p}}\right)\left(\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}+1\right).$

Next let $w\in\gamma_{h}^{\\#}T{\mathcal{M}}$ be piecewise linear in the sense
of definition 18, then applying integration by parts twice we obtain

$\int_{0}^{1}(-D_{t}^{2}\log_{\gamma_{h}}\\!\gamma,w)\,\mathrm{d}t=\sum_{k=1}^{N}\int_{t_{k-1}}^{t_{k}}(\log_{\gamma_{h}}\\!\gamma,{\mathrm{R}}(w,\dot{\gamma}_{h})\dot{\gamma}_{h})\,\mathrm{d}t\leq
c\|d(\gamma_{h},\gamma)\|_{L^{2}}\|w\|_{L^{2}}\|\dot{\gamma}_{h}\|_{L^{\infty}}^{2}\leq\tilde{c}\|d(\gamma_{h},\gamma)\|_{L^{2}}\|w\|_{L^{2}}.$

By the above we have

$\int_{0}^{1}(D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma,w)\,\mathrm{d}t\leq
c\|w\|_{L^{2}}\left(h^{2}+\|d(\gamma_{h},\gamma)\|_{L^{2}}\right)\left(\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}+1\right).$

The choice $w=D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma$ now
finishes the proof. ∎

######  Corollary 24 (Uniform boundedness of $D_{t}^{2}\gamma_{h}$).

For $\gamma$ four times differentiable and $h$ small enough (depending on
${\mathcal{M}}$ and $\gamma$) the cubic spline interpolation satisfies

$\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}\leq
3\|D_{t}^{2}\gamma\|_{L^{\infty}}+ch^{\frac{3}{2}}\qquad\text{and}\qquad\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{2}},\|D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma\|_{L^{2}}\leq
ch^{2}$

for a constant $c$ only depending on ${\mathcal{M}}$ and $\gamma$.

###### Proof.

Using lemma 19 we have

$\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}\leq\|D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma\|_{L^{\infty}}+\|{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma\|_{L^{\infty}}\leq\|D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma\|_{L^{\infty}}+2\|D_{t}^{2}\gamma\|_{L^{\infty}}.$

By lemma 20 and $D_{t}^{2}\gamma_{h}={\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma_{h}$
we have

$\displaystyle\|D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma\|_{L^{\infty}}$
$\displaystyle\leq
c\left(1+h^{-\frac{1}{2}}\right)\|D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma\|_{L^{2}}.$

Using proposition 23 and corollary 15 this turns into

$\displaystyle\|D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma\|_{L^{\infty}}$
$\displaystyle\leq
c\left(h^{2}+h^{\frac{3}{2}}\right)(\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}+1)\;,$

and thus by lemma 19

$\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}\leq
2ch^{\frac{3}{2}}(1+\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}})+2\|D_{t}^{2}\gamma\|_{L^{\infty}}.$

Taking $h$ small enough such that $ch^{\frac{3}{2}}\leq\frac{1}{6}$, we arrive
at the desired bound for $\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}}$. This bound in
turn can be applied in proposition 23 together with corollary 15 to yield the
desired bound on $\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{2}}$ and
$\|D_{t}^{2}\gamma_{h}-{\mathcal{I}^{1}_{0}}D_{t}^{2}\gamma\|_{L^{2}}$. ∎

######  Remark 25 (Boundedness of $\dot{\gamma}_{h}$ for linear splines).

Note that a slightly weaker boundedness result than proposition 16 for linear
splines could have been obtained following exactly the same argument as the
one above for cubic splines. In that case, one would have used a piecewise
constant interpolation operator ${\mathcal{I}^{0}_{0}}$ on vector fields
$v\in\gamma^{\\#}T{\mathcal{M}}$ which would simply parallel transport
$v(t_{k-1})$ along $\gamma((t_{k-1},t_{k}))$ and which would satisfy analogous
properties to those of ${\mathcal{I}^{1}_{0}}$ in lemma 20. In the subsequent
estimate of

$\dot{\gamma}_{h}-{\mathcal{I}^{0}_{0}}\dot{\gamma}=-D_{t}\log_{\gamma_{h}}\\!\gamma+\left[\dot{\gamma}_{h}+\partial_{1}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma}_{h})\right]+\left[\partial_{2}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma})-{\mathcal{I}^{0}_{0}}(\partial_{2}\log_{\gamma_{h}}\\!\gamma(\dot{\gamma}))\right]$

one would this time exploit the orthogonality property
$\int_{0}^{1}(D_{t}\log_{\gamma_{h}}\\!\gamma,w)\,\mathrm{d}t=0$ for any
piecewise constant vector field along $\gamma_{h}$.

With the above preparations we can now turn to the final convergence result.
Corollary 24 and proposition 14 already prove the desired result
$\|d(\gamma_{h},\gamma)\|_{L^{2}}\leq ch^{4}$ in the $L^{2}$-norm; for the
$L^{\infty}$-norm (or any other $L^{p}$-norm) we will proceed by bounding
$\|D_{t}^{4}\log_{\gamma_{h}}\\!\gamma\|_{L^{\infty}}$.

######  Corollary 26 (Quartic convergence).

For $\gamma$ four times differentiable, $p\in[1,\infty]$, and $h$ small enough
(depending on ${\mathcal{M}}$ and $\gamma$) the cubic spline interpolation
satisfies

$\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{p}}\leq
ch^{2},\qquad\|d(\gamma_{h},\gamma)\|_{L^{p}}\leq ch^{4},$

where the constant $c>0$ only depends on ${\mathcal{M}}$, $\gamma$, and $p$.

###### Proof.

We just show the $L^{\infty}$ estimate, which is the strongest among all $p$,
and apply corollary 12 to $g=\pi_{\gamma_{h}}\log_{\gamma_{h}}\\!\gamma$ to
obtain the desired result. To this end it remains to show uniform boundedness
of $|\partial_{t}^{4}g|=|D_{t}^{4}\log_{\gamma_{h}}\\!\gamma|$ on each
interval $I_{k}=(t_{k-1},t_{k})$. Now

$D_{t}^{4}\log_{\gamma_{h}}\\!\gamma=\partial_{1}\log_{\gamma_{h}}\\!\gamma(D_{t}^{4}\gamma_{h})+c_{1}\partial_{1}^{2}\log_{\gamma_{h}}\\!\gamma(D_{t}^{3}\gamma_{h},D_{t}\gamma)+c_{2}\partial_{2}\partial_{1}\log_{\gamma_{h}}\\!\gamma(D_{t}^{3}\gamma_{h},D_{t}\gamma_{h})+r$

with some integers $c_{1},c_{2}$, where $r$ is a term containing derivatives
of the logarithm up to fourth order evaluated in directions of derivatives of
$\gamma$ up to fourth order and of $\gamma_{h}$ up to second order. Thus,
$|r|$ can be bounded by a constant (depending on ${\mathcal{M}}$ and
$\gamma$). Consequently, exploiting lemma 6 it remains to show uniform
boundedness of $|D_{t}^{4}\gamma_{h}|$ and
$d(\gamma_{h},\gamma)|D_{t}^{3}\gamma_{h}|$, where by corollary 17 the latter
is no larger than $h^{2}|D_{t}^{3}\gamma_{h}|$ up to a constant factor (in
fact, below we will even obtain boundedness of $h|D_{t}^{3}\gamma_{h}|$). By
proposition 10(2) as well as proposition 8(2) and corollary 24 we have

$|D_{t}^{4}\gamma_{h}|=|{\mathrm{R}}(D_{t}^{2}\gamma_{h},\dot{\gamma}_{h})\dot{\gamma}_{h}|\leq
C$

for a constant $C>0$ depending on ${\mathcal{M}}$ and $\gamma$. Using the
scaled Gagliardo–Nirenberg–Sobolev inequality or an inverse inequality, we
then obtain

$\|D_{t}^{3}\gamma_{h}\|_{L^{\infty}(I_{k})}\leq
c\|D_{t}^{4}\gamma_{h}\|_{L^{\infty}(I_{k})}^{1/2}\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}(I_{k})}^{1/2}+ch^{-1}\|D_{t}^{2}\gamma_{h}\|_{L^{\infty}(I_{k})}\leq
ch^{-1},$

which finally proves boundedness $|\partial_{t}^{4}g|<C$ for a constant $C$
only depending on ${\mathcal{M}}$ and $\gamma$.

The estimate of $\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{\infty}}$ now
follows again from the the scaled Gagliardo–Nirenberg–Sobolev inequality.
Indeed, there exists $c>0$ such that on each interval $I_{k}=(t_{k-1},t_{k})$

$\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{\infty}(I_{k})}\leq
c\|D_{t}^{4}\log_{\gamma_{h}}\\!\gamma\|_{L^{\infty}(I_{k})}^{1/2}\|\log_{\gamma_{h}}\\!\gamma\|_{L^{\infty}(I_{k})}^{1/2}+ch^{-2}\|\log_{\gamma_{h}}\gamma\|_{L^{\infty}(I_{k})}.$

Using $\|d(\gamma_{h},\gamma)\|_{L^{\infty}}\leq Ch^{4}$ as well as the
uniform boundedness of $D_{t}^{4}\log_{\gamma_{h}}\\!\gamma$, the estimate
turns into

$\|D_{t}^{2}\log_{\gamma_{h}}\\!\gamma\|_{L^{\infty}(I_{k})}\leq
c(\sqrt{C}\|D_{t}^{4}\log_{\gamma_{h}}\\!\gamma\|_{L^{\infty}(I_{k})}^{1/2}+C)h^{2}.\qed$

## Acknowledgements

B. Wirth’s research was supported by the Alfried Krupp Prize for Young
University Teachers awarded by the Alfried Krupp von Bohlen und Halbach-
Stiftung. He also acknowledges support by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany’s Excellence Strategy EXC
2044–390685587, Mathematics Münster: Dynamics–Geometry–Structure, and under
the Collaborative Research Centre 1450, InSight, University of Münster.

## References

  * [1] Philippe G. Ciarlet.  The finite element method for elliptic problems.  North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978.  Studies in Mathematics and its Applications, Vol. 4.
  * [2] Carl de Boor.  Bounding the error in spline interpolation.  SIAM Review, 16(4):531–544, 1974.
  * [3] Carl de Boor.  A bound on the $L_{\infty}$-norm of $L_{2}$-approximation by splines in terms of a global mesh ratio.  Math. Comp., 30(136):765–771, 1976.
  * [4] Jim Douglas, Jr., Todd Dupont, and Lars Wahlbin.  Optimal $L_{\infty}$ error estimates for Galerkin approximations to solutions of two-point boundary value problems.  Math. Comp., 29:475–483, 1975.
  * [5] Hanne Hardering.  Intrinsic Discretization Error Bounds for Geodesic Finite Elements.  PhD thesis, Freie Universität Berlin, 2015.
  * [6] Behrend Heeren, Martin Rumpf, and Benedikt Wirth.  Variational time discretization of Riemannian splines.  IMA J. Numer. Anal., 39(1):61–104, 2019.
  * [7] H. Karcher.  Riemannian center of mass and mollifier smoothing.  Comm. Pure Appl. Math., 30:509–541, 1977.