A Numerical Algorithm for C2-splines on Symmetric Spaces

TL;DR

This paper introduces an efficient iterative numerical algorithm for constructing C2-splines on Riemannian symmetric spaces, addressing a key challenge in manifold interpolation with applications in computer vision and quantum control.

Contribution

It provides the first convergent iterative method for C2-splines on symmetric spaces, surpassing existing approaches in efficiency and parallelizability.

Findings

Algorithm converges linearly.

Applicable to spheres, complex projective spaces, and Grassmannians.

Outperforms Riemannian cubic-based methods in computational efficiency.

Abstract

Cubic spline interpolation on Euclidean space is a standard topic in numerical analysis, with countless applications in science and technology. In several emerging fields, for example computer vision and quantum control, there is a growing need for spline interpolation on curved, non-Euclidean space. The generalization of cubic splines to manifolds is not self-evident, with several distinct approaches. One possibility is to mimic the acceleration minimizing property, which leads to Riemannian cubics. This, however, require the solution of a coupled set of non-linear boundary value problems that cannot be integrated explicitly, even if formulae for geodesics are available. Another possibility is to mimic De~Casteljau's algorithm, which leads to generalized B\'ezier curves. To construct C2-splines from such curves is a complicated non-linear problem, until now lacking numerical methods. Here we provide an iterative algorithm for C2-splines on Riemannian symmetric spaces, and we prove convergence of linear order. In terms of numerical tractability and computational efficiency, the new method surpasses those based on Riemannian cubics. Each iteration is parallel, thus suitable for multi-core implementation. We demonstrate the algorithm for three geometries of interest: the $n$-sphere, complex projective space, and the real Grassmannian.

Figures

Tables

Table 1. Table 1: Some examples of Riemannian symmetric spaces. For the details of the corresponding symmetric space decomposition, see  [ 10 , § 5] . In the right-most column, we give the geometry class to use when interpolating with our Python package; we refer to Appendix D for example usage.

	$M$	$G$	$H$	Implementation
Euclidean space	${ℝ}^n$	$\mathrm{SE}\left(n\right)$	$\mathrm{SO}\left(n\right)$	Flat
$n$-sphere	$S^n$	$O\left(n+1\right)$	$O\left(n\right)$	Sphere
Hyperbolic space	$H^n$	$O(n,1)$	$O\left(n\right)$	Hyperbolic
Real Grassmannian	$\mathrm{Gr}(k,n)$	$O\left(n\right)$	$O\left(k\right)\times O\left(n-k\right)$	Grassmannian
Orthogonal group	$O\left(n\right)$	$O\left(n\right)\times O\left(n\right)$	$O\left(n\right)$
Complex projective space	$ℂP^n$	$U\left(n+1\right)$	$U\left(n\right)\times U\left(1\right)$	Projective
Complex Grassmannian	${\mathrm{Gr}}_{ℂ}(k,n)$	$U\left(n\right)$	$U\left(k\right)\times U\left(n-k\right)$
Unitary group	$U\left(n\right)$	$U\left(n\right)\times U\left(n\right)$	$U\left(n\right)$
Co-variance matrices	$P\left(n\right)$	$\mathrm{GL}\left(n\right)$	$O\left(n\right)$

Equations

$$\underset{\gamma (i)=p_i}{\min }{\int }_0^N|{\nabla }_{\dot{\gamma }(t)}\dot{\gamma }(t){|}^2\mathrm{d}t.$$

$$G\times M\ni (g,p)\mapsto g\cdot p\in M.$$

$$H_p:=\{g\in G|g\cdot p=p\}.$$

$$\mathfrak{g}={\mathfrak{h}}_p\oplus {\mathfrak{m}}_p,$$

$$[{\mathfrak{h}}_p,{\mathfrak{h}}_p]\subset {\mathfrak{h}}_p,[{\mathfrak{h}}_p,{\mathfrak{m}}_p]\subset {\mathfrak{m}}_p,[{\mathfrak{m}}_p,{\mathfrak{m}}_p]\subset {\mathfrak{h}}_p.$$

$${\mathfrak{h}}_{g\cdot p}=g\cdot {\mathfrak{h}}_p,{\mathfrak{m}}_{g\cdot p}=g\cdot {\mathfrak{m}}_p,$$

$$(\xi ,p)\mapsto \xi \cdot p\in T_{p}M.$$

$$\eta (t)=\exp (t\xi )\cdot p_0\eta (1)=p_1\xi \in {\mathfrak{m}}_{p_0}.$$

$$[p_0,p_1{]}_t:=\eta (t).$$

$$[p_0,p_1{]}_0$$

$$\gamma(t)\coloneqq\Big{[}\big{[}[p_{0},q_{0}]_{t},[q_{0},q_{1}]_{t}\big{]}_{t},\big{[}[q_{0},q_{1}]_{t},[q_{0},p_{1}]_{t}\big{]}_{t}\Big{]}_{t}.$$

$$g\cdot x\in M.$$

$$E_p(\xi \cdot p)=\exp (\xi )\xi \in {\mathfrak{m}}_p.$$

$$v\in T_{p}M,E_p(v)\cdot p=q,⟺v={\log }_p(q).$$

$${\log }_{p_i}(q_i^{+})$$

$${\log }_{p_i}(E_{p_i}(-v_i)\cdot q_{i+1}^{-})-{\log }_{p_i}(q_i^{+})$$

$$q_i^{+}-p_i=-(q_i^{-}-p_i)$$

$$q_{i-i}^{+}-2q_i^{-}+p_i=q_{i+1}^{-}-2q_i^{+}+p_i.$$

$$\begin{split}&w_{i}(v_{i-1},v_{i},v_{i+1})\coloneqq\\ &\qquad\log_{p_{i}}\big{\lparen}E_{p_{i}}(-v_{i})E_{p_{i+1}}(-v_{i+1})\cdot p_{i+1}\big{\rparen}\\ &\qquad-\log_{p_{i}}\big{\lparen}E_{p_{i}}(v_{i})E_{p_{i-1}}(v_{i-1})\cdot p_{i-1}\big{\rparen}.\end{split}$$

$$\mathbf{v}\mapsto \frac{1}{4}\mathbf{w}(\mathbf{v})+\frac{1}{2}\mathbf{v},$$

$$\bm{w}(\bm{v})\coloneqq\big{\lparen}w_{1}(v_{0},v_{1},v_{2}),\dotsc,w_{N-1}(v_{N-2},v_{N-1},v_{N})\big{\rparen}.$$

$$v_0(v)=v_{\text{start}},v_N(v)=v_{\text{end}}$$

$$v_0(v)=\frac{1}{2}{\log }_{p_0}(q_1^{-}),v_N(v)=-\frac{1}{2}{\log }_{p_N}(q_{N-1}^{+}).$$

$$D=\underset{0\le i\le N-1}{\max }d(p_i,p_{i+1}).$$

$$K^2=\underset{\xi ,\eta ,\zeta \in {\mathfrak{m}}_p\backslash 0}{\sup }\frac{\Vert [\xi ,[\eta ,\zeta ]]\Vert }{\Vert \xi \Vert \Vert \eta \Vert \Vert \zeta \Vert }.$$

$$q_0^2+q_1^2+q_2^2+q_3^2=1.$$

$${\mathrm{rot}}_q(p)=qp{q}^{-1}$$

$$\{{\mathrm{rot}}_{\gamma (t)}(i),{\mathrm{rot}}_{\gamma (t)}(j),{\mathrm{rot}}_{\gamma (t)}(k)\}.$$

$$\big{\{}\,(x,y,z)\in{\mathbb{R}}^{2}\;\big{|}\;z\leq 0\,\big{\}}.$$

$$(x,y,z)\mapsto\Big{\lparen}\frac{x}{1-z},\frac{y}{1-z}\Big{\rparen}.$$

$${\Psi }^p(g)={\tilde{L}}_g(p)=g\cdot p,$$

Full text

MSC-2010: 41A15, 65D07, 65D05, 53C35, 53B20, 14M17

Keywords: De Casteljau, cubic spline, Riemannian symmetric space, Bézier curve

1. 1 Introduction
2. 2 Generalized Bézier curves on symmetric spaces
3. 3 Numerical algorithm
4. 4 Examples
5. 5 Outlook
6. A Identities and bounds on symmetric spaces
7. B $C^{2}$ continuity
8. C Convergence of the fixed-point method
9. D Python implementation

1 Introduction

We address the following.

Problem 1.1**.**

Given a set of points $\{\,p_{i}\,\}_{i=0}^{N}$ on a connected manifold $M$, construct a $C^{2}$ path $\gamma\colon[0,N]\to M$ such that $\gamma(i)=p_{i}$.

For $M={\mathbb{R}}^{n}$ every textbook in numerical analysis teaches cubic splines, i.e., piecewise polynomials of order 3 with matching first and second derivatives. However, the generalization of cubic splines to curved space is non-trivial, essentially because polynomials are not well-defined on manifolds. Interpolating paths on manifolds are, nevertheless, needed in a growing number of applications. In the realization of quantum computers, for example, quantum control of qubits leads to an interpolation problem in complex projective space ${\mathbb{C}}P^{n}$ (see §​ 4.2). In computer vision, as another example, the recognition of a point cloud up to affine transformations is naturally identified with an element in the Grassmannian $\mathrm{\mathrm{Gr}}(k,n)$ (see §​ 4.3).

If $M$ is a Riemannian manifold, a natural generalization of cubic splines are Riemannian cubics [12, 2]. They are based on minimizing acceleration under interpolation constraints, by solving the problem

[TABLE]

The corresponding Euler–Lagrange equation is a fourth order ODE on $M$ with complicated, coupled boundary conditions [12]. In addition to the Euclidean case, explicit solutions to the initial-value problem are known for so called null Lie quadratics on $\mathrm{SO}(3)$ and $\mathrm{SO}(2,1)$ equipped with the bi-invariant Riemannian metric [11]. For other cases, even when geodesics are explicitly known, one is left with numerical ODE methods, combined, for example, with a shooting method to match the boundary conditions. This approach is computationally costly and the associated convergence analysis is involved. Another possibility is to only approximatelly fulfill the interpolation constraints by incorporating them in the minimization functional and then use a steepest-descent method on the space of curves [15]. This approach is also costly.

If $M$ is a manifold for which the geodesic boundary problem can be solved explicitly, there is a natural concept of generalized Bézier curves (GBC). They are based on a direct generalization of De Casteljau’s algorithm [13]. Splines can then be constructed by gluing piecewise GBC with matching boundary conditions. In this way, a fully explicit algorithm for $C^{1}$-splines on Grassmannian and Stiefel manifolds has been developed [8]. The $C^{2}$ condition, however, first derived by Popiel and Noakes [14], is significantly more complicated than the $C^{1}$ condition. For this reason, there is a lack of numerical algorithms for $C^{2}$-splines. Nevertheless, $C^{2}$ continuity is often required in applications, to avoid discontinuities in accelerations.

In this paper we give the first numerical algorithm for $C^{2}$-splines based on GBC. Instead of treating all Riemannian manifolds, we focus on the subclass of symmetric spaces (see §​ 2.1). Such manifolds are of high interest in applications, and include, e.g., spheres, hyperbolic space, real and complex projective spaces, and Grassmannian manifolds (see Table 1). In many cases, the geodesic boundary value problem on a symmetric space has an explicit solution (a prerequisite for De Casteljau’s algorithm). The original $C^{2}$ condition of Popiel and Noakes is, however, still complicated. We now list the specific contributions of our paper.

Using the special homogeneous space structure of symmetric spaces, we give a significant simplification of the $C^{2}$ condition (Theorem 2.2). 2. 2.

Using the new $C^{2}$ condition, we construct an algorithm based on fixed-point iterations (Algorithm 1). 3. 3.

We prove convergence of the algorithm under conditions on the maximal distance between consecutive interpolation points (Theorem 3.1).

The paper is organized as follows. In §​ 2 we first give a brief description of Riemannian symmetric spaces and thereafter describe the construction of generalized Bézier curves. We also formulate the $C^{2}$ condition. In §​ 3 we give the numerical algorithm and we formulate the convergence result. Numerical examples are given in §​ 4. An outlook towards future work is given in §​ 5. The proofs of the $C^{2}$ condition and the convergence result are long and technical, and therefore given in Appendix A-C. In Appendix D we give a brief demonstration of an easy-to-use open-source Python package implementing our interpolation algorithm for several geometries.

2 Generalized Bézier curves on symmetric spaces

Let us describe the construction of generalized Bézier curves on a symmetric space. The construction relies on a notion of interpolation between two points on the manifold, generalizing the notion of a geodesic on a Riemannian manifold [9].

2.1 Symmetric spaces

In this section we briefly recall symmetric spaces. A list of the most common examples of symmetric spaces is given in Table 1. For details, we refer to [16] or [4]. We assume that a Lie group $G$ acts transitively on $M$, making $M$ a homogeneous space [16, § 4]. We denote the action by

[TABLE]

For every point $p\in M$, we define the isotropy subgroup $H_{p}$ at $p$

[TABLE]

Now, the homogeneous space $M$ is called a symmetric space if the Lie algebra $\mathfrak{g}$ of $G$ has a direct sum decomposition

[TABLE]

where $\mathfrak{h}_{p}$ is the Lie algebra of $H_{p}$, and the summands fulfill the following algebraic conditions

[TABLE]

If these conditions are satisfied at one point $p\in M$, then they are automatically valid at every point with

[TABLE]

where we follow the convention formulated in Appendix A that $g\cdot$ denotes the canonical action on the space at hand (here, the adjoint action of a Lie group on its Lie algebra).

By differentiating the action map $(g,p)\mapsto g\cdot p$ with respect to the arguments, we also obtain, on the one hand, an infinitesimal action of $\mathfrak{g}$ on $M$, and, on the other hand, a prolonged action of $G$ on $TM$. For $\xi\in\mathfrak{g}$ and $p\in M$, we denote the infinitesimal action by

[TABLE]

For a fixed $p$, this map induces an isomorphism between $\mathfrak{m}_{p}$ and $T_{p}M$.

2.2 Interpolating Curves

We now define the interpolating curve $\eta\colon[0,1]\to M$ between two (close enough) points $p_{0}$ and $p_{1}$ as the unique curve satisfying

[TABLE]

Let us introduce the notation

[TABLE]

In particular, notice that

[TABLE]

In many (but not all) cases, a symmetric space $M$ can be equipped with the structure of a Riemannian manifold, given through an $H_{p}$-invariant inner product $\langle\cdot,\cdot\rangle_{p}$ on $\mathfrak{m}_{p}$. In this case, the interpolating curves are geodesics. From here on, we shall assume that $M$ is such a Riemannian symmetric space.

2.3 De Casteljau’s construction

We briefly describe the De Casteljau construction, illustrated on Figure 1. This algorithm constructs Bézier curves from interpolating curves. We focus on cubic Bézier curves. From four points $p_{0}$, $q_{0}$, $q_{1}$ and $p_{1}$, one defines the cubic Bézier curve $\gamma\colon[0,1]\to M$ by

[TABLE]

Note that, by construction, $\gamma(0)=p_{0}$ and $\gamma(1)=p_{1}$.

2.4 Exponential and logarithm on symmetric spaces

For the $C^{2}$ condition in §​ 2.5 and implementation of the algorithm we need the following functions (see Figure 2).

The action of a group element $g\in G$ on a point $x\in M$:

[TABLE] 2. 2.

A movement function $E$, which, to a point $p$ and a velocity $v$ at $p$, assigns an element $E_{p}(v)\in G$. This function should generate interpolating curves (6), i.e., it should be of the form

[TABLE] 3. 3.

The logarithm; given two points $p$ and $q$ in $M$, it is defined by

[TABLE]

Throughout the paper we make the blanket assumption that $p$ and $q$ are never past conjugate points, so that the logarithm is always well-defined.

2.5 Cubic splines and $C^{2}$ condition

We shall now consider curves consisting of piecewise cubic Bézier curves. As we have seen in §​ 2.3, each cubic Bézier curve is determined by two interpolating points and two control points. To construct a curve of piecewise Bézier curves we therefore need to impose conditions on the control points (in addition to the interpolating conditions). On Euclidean space, the standard approach is to use the conditions for $C^{2}$ continuity at the interpolating points, which leads to a linear set of equations. On a general Riemannian manifold, the corresponding $C^{2}$ condition, given by Popiel and Noakes [14], is highly nonlinear, involving the inverse of the derivative of the Riemannian exponential. In the special case of symmetric spaces, however, we now show that the $C^{2}$ condition simplifies significantly, involving only the three operations (9), (10), and (11).

Definition 2.1**.**

Let $\bm{p}=(p_{0},\ldots,p_{N})$ be points on a symmetric space $M$. A composite cubic Bézier curve on $M$ is a curve $\gamma\colon[0,N]\to M$ such that for each integer $0\leq i\leq N$, $\gamma(i)=p_{i}$, and such that for each integer $0\leq i<N$, the restriction $\gamma_{i}\coloneqq\gamma|_{[i,i+1]}$ is a cubic Bézier curve (as defined in §​ 2.3). To each such piece $\gamma_{i}$, we associate two control points, denoted $q_{i}^{+}$ and $q_{i+1}^{-}$. A $C^{2}$-continuous composite cubic Bézier curve is called a cubic spline.

The interpolating conditions ensure that $\gamma$ is continuous (or $C^{0}$). However, for higher degree of continuity ($C^{1}$ or $C^{2}$), we need additional conditions on the control points. We now give such conditions.

Theorem 2.2**.**

Let $M$ be a Riemannian symmetric space and let $\gamma\colon[0,N]\to M$ be a composite cubic Bézier curve according to §​ 2.5. Consider the conditions on the control points given by

[TABLE]

1. (i)

If condition (12) is fulfilled, then $\gamma$ is a $C^{1}$-curve. 2. (ii)

If conditions (12) and (13) are fulfilled, then $\gamma$ is a $C^{2}$-curve, i.e., a cubic spline.

Furthermore, if $\gamma$ fulfills the $C^{2}$ conditions (12) and (13), then it is uniquely determined either by clamped boundary conditions, where $\gamma^{\prime}(0)$ and $\gamma^{\prime}(N)$ are prescribed, or by natural boundary conditions, where $\nabla_{\gamma^{\prime}(0)}\gamma^{\prime}(0)=\nabla_{\gamma^{\prime}(0)}\gamma^{\prime}(0)=0$.

Remark 2.3**.**

Notice how (12) and (13) are generalizations of the first and second order finite differences in Euclidean geometry. Indeed, in this case condition (12) reads

[TABLE]

and condition (13) reads

[TABLE]

See Figure 3 for a geometric illustration of the $C^{2}$ condition (13).

Outline of proof.

The details of the proof are rather technical and therefore given in Appendix B. However, for convenience we also provide an outline here:

1. (i)

The $C^{2}$ condition in the case of Riemannian symmetric spaces is given by Popiel and Noakes [14, Sec. 4], expressed in terms of the derivative of the symmetry function $I_{p}(q)=E_{p}(-\log_{p}(q))$. 2. (ii)

Using the special symmetric space structure, we show that the derivative of $I_{p}$ can be expressed solely in terms of the movement function $E_{p}$ (Appendix A). 3. (iii)

Together with a result on equivariance of the $\log$ function (Appendix A) we can then simplify the original Riemannian symmetric space $C^{2}$ condition by Popiel and Noakes [14] to the one given by (13).

∎

3 Numerical algorithm

We shall now construct a fixed point iteration algorithm for cubic splines on symmetric spaces, based on the $C^{2}$ condition in Theorem 2.2. To this end, we parameterize the control points $q_{i}^{-}$ and $q_{i}^{+}$ by the corresponding tangent vectors $v_{i}$ as defined in Theorem 2.2. Next, let

[TABLE]

The iteration map over $\bm{v}\coloneqq(v_{1},\dotsc,v_{N-1})$ is then given by

[TABLE]

where

[TABLE]

Pseudo code for the complete algorithm is given in Algorithm 1. We consider two options for boundary conditions.

1. (i)

Clamped spline:

[TABLE]

where $v_{\text{start}}$ and $v_{\text{end}}$ are prescribed constants. 2. (ii)

Natural spline:

[TABLE]

3.1 Convergence result

In this section we give a result on convergence of the fixed point method given by Algorithm 1. To do so, we first give some preliminary definition.

The Riemannian distance between $p,q\in M$ is denoted by $d(p,q)$. Our proof of convergence uses that two consecutive points $p_{i},p_{i+1}$ are close. Thus, we define the constant

[TABLE]

The other constant that comes into play is a bound on (essentially) the curvature of $M$. Indeed, if $\lVert\cdot\rVert_{p}$ denotes the $H_{p}$-invariant norm on $\mathfrak{m}_{p}$ associated with the Riemannian structure of $M$, then we define the constant $K\geq 0$ as

[TABLE]

If $M$ is flat, then $K=0$. If $M$ is the $n$-sphere of radius $r$, then $K=1/r$.

Our convergence result is now formulated as follows.

Theorem 3.1**.**

If $KD$ is sufficiently small, then Algorithm 1 with natural boundary conditions converges linearly. If $KD$ and $\max\{\,K\lVert v_{0}\rVert,K\lVert v_{N}\rVert\,\}$ are both sufficiently small, then Algorithm 1 with clamped boundary conditions converges linearly.

Outline of proof.

Again, the proof is long and technical, and therefore given in Appendix C. Here we give a brief outline.

The proof is based on showing that the iteration mapping (17) is a contraction. For convenience, we use the variables $\xi_{i}\in\mathfrak{m}_{p_{i}}$ instead of $v_{i}\in T_{p_{i}}M$. The iteration mapping is denoted $\phi$.

1. (i)

The first step is to give an invariant region of the iteration mapping. This is given in §​ C.1, which proves existence of $V>0$ (depending on both $D$ and $K$), such that $\lVert\bm{\xi}\rVert\leq V$ implies $\lVert\phi(\bm{\xi})\rVert\leq V$. 2. (ii)

The next step is to prove that $\phi$ is a contraction, which is established by a series of estimates on the partial derivatives on $\phi$ (§​ C.2, §​ C.2, and §​ C.2). These estimates depend on $K$. 3. (iii)

The final step, in §​ C.3, is to use the contraction mapping theorem.

∎

4 Examples

4.1 Unit quaternions

Unit quaternions (versors) are used extensively in computer graphics to represent 3-D rotations. They are elements of the form $q=q_{0}+q_{1}i+q_{2}j+q_{3}k$ with

[TABLE]

Thus, we can identify unit quaternions with $S^{3}$. In turn, $S^{3}$ is a double cover of $\mathrm{SO}(3)$. A point on $S^{3}$ therefore gives a rotation matrix, and likewise any rotation matrix corresponds to two antipodal points on $S^{3}$. The rotation of a vector $p=(p_{x},p_{y},p_{z})$ by a unit quaternion $q$ can be compactly written using quaternion multiplication as

[TABLE]

where $p=p_{x}i+p_{y}j+p_{z}k$ is thought of as a pure imaginary quaternion.

Since $S^{3}\simeq\mathrm{O}(4)/\mathrm{O}(3)$ is a Riemannian symmetric space (with respect to the standard Riemannian metric), we can use Algorithm 1 to obtain $C^{2}$ continuous spline interpolation between rotations. Since geodesics on $S^{3}$ are given by great circles, it is straightforward to derive the mappings $E$ and $\log$.

The resulting $C^{2}$-curve interpolating 5 orientations is shown in Figure 5. In the figure, an element $\gamma(t)\in S^{3}$ is represented by the rotated basis vectors

[TABLE]

The actual interpolation points are marked with bolder lines.

Note that $S^{3}$ can be canonically identified with $\mathrm{SU}(2)$. Note also that $\mathrm{SU}(2)$ has a symmetric space structure as the quotient $(\mathrm{SU}(2)\times\mathrm{SU}(2))/\mathrm{SU}(2)$ [10]. Here, however, we considered the symmetric space structure $S^{3}=\mathrm{SO}(4)/\mathrm{SO}(3)$. These two symmetric space structure are locally equivalent. This is because, $\mathrm{SU}(2)\times\mathrm{SU}(2)$ is a double cover of $\mathrm{SO}(4)$ and $S^{3}$ is a double cover of $\mathrm{SO}(3)$. Our algorithm gives the same result regardless of the choice of either of those two symmetric space structures on $S^{3}$. From this point of view, our algorithm could be used for spline interpolation on any compact Lie group $G$, using the Cartan–Schouten symmetric space structure $(G\times G)/G$.

4.2 Quantum states

The control of quantum states is an important subproblem in quantum information and the realization of quantum computers [3]. The objective is to find a time-dependent Hamiltonian $H(t)$ designed to ‘steer’ a given quantum state $|\psi_{0}\rangle$ through a sequence of states $|\psi_{1}\rangle,|\psi_{2}\rangle,\ldots,|\psi_{N}\rangle$ at given times $t_{1},\ldots,t_{N}$. As instantaneous switching of the Hamiltonian is not experimentally feasible [1], the interpolating curve $|\psi(t)\rangle$ should be at least $C^{2}$ continuous.

If the quantum state space is finite dimensional, corresponding to an ensemble of qubits, then, in the geometric description of quantum mechanics [6, 7], the phase space is given by complex projective space ${\mathbb{C}}P^{n}$. The quantum control problem can then be seen as a two-step process:

Find an interpolating curve $t\mapsto|\psi(t)\rangle\in{\mathbb{C}}P^{n}$ such that $|\psi(0)\rangle=|\psi_{0}\rangle$ and $|\psi(t_{k})\rangle=|\psi_{k}\rangle$. 2. 2.

Using the homogeneous structure ${\mathbb{C}}P^{n}\simeq\mathrm{U}(n+1)/\mathrm{U}(n)\times\mathrm{U}(1)$, lift $|\psi(t)\rangle$ to a curve $t\mapsto g(t)\in\mathrm{U}(n+1)$ such that $g(0)=e$ and $\pi(g(t))=|\psi(t)\rangle$.111The lifting is, of course, not unique. It is natural to minimize the change in Hamiltonian $\dot{H}$, as suggested by Brody, Holm, and Meier [1]. The time-dependent Hamiltonian is then given by $H(t)=-i\dot{g}(t)g(t)^{-1}$.

We can use Algorithm 1 for the first step. (The second step is not treated in this paper.)

The simplest case of a single qubit corresponds to the phase space ${\mathbb{C}}P^{1}$, which is isomorphic to a sphere (called the Bloch sphere). Using 6 interpolation points, the resulting interpolating curve is visualized on the Bloch sphere in Figure 6.

4.3 Affine shapes

The space of affine shapes is the quotient set whose elements consist of $k+1$ labelled (i.e. ordered) points in ${\mathbb{R}}^{n}$ defined up to affine transformations. Thus, $(p_{0},\ldots,p_{k})\in({\mathbb{R}}^{n})^{k+1}$ and $(\tilde{p}_{0},\ldots,\tilde{p}_{k})\in({\mathbb{R}}^{n})^{k+1}$ define the same affine shape if there exists an affine transformation $p\mapsto Ap+q$ such that $(\tilde{p}_{0},\ldots,\tilde{p}_{k})=(Ap_{0}+q,\ldots,Ap_{k}+q)$. The study of affine shapes has many applications in computer vision [5]. An example is the identification of planar surfaces (such as facades of buildings) in images taken from different projection angles.

The space of affine shapes with full rank (meaning that the $k+1$ points are not contained in an affine hyperplane of ${\mathbb{R}}^{n}$) can be identified with the real Grassmannian manifold $\mathrm{Gr}(n,k)$ of $n$-dimensional subspaces in ${\mathbb{R}}^{k}$. Indeed, choose a coordinate system and place the coordinates of the vectors $p_{1}-p_{0},\dotsc,p_{k}-p_{0}$ column by column in an $n\times k$ matrix. The kernel of this matrix is then a subspace of ${\mathbb{R}}^{k}$. The full rank assuption means that the kernel has dimension $k-n$. This gives us a point in $\mathrm{Gr}(k-n,k)$, which, by the canonical isomorphism between $\mathrm{Gr}(k-n,k)$ and $\mathrm{Gr}(n,k)$, gives a point in $\mathrm{Gr}(n,k)$.

Since the Grassmannian $\mathrm{Gr}(n,k)$ is a symmetric space, we can use our algorithm to construct $C^{2}$-splines between affine shapes.

Let us study a particular case of affine shapes corresponding to $3+1$ points in ${\mathbb{R}}^{2}$. Our aim is to create a spline between such shapes, defined by (A)–(D) in Figure 7. For every shape, we obtain a $2\times 3$ matrix of full rank, which defines our element in $\mathrm{Gr}(2,3)$. We then use Algorithm 1 to construct a spline interpolating the points in $\mathrm{Gr}(2,3)$.

To visualize the result, notice that $\mathrm{Gr}(2,3)$ is isomorphic to the real projective plane ${\mathbb{R}}P^{2}$. Indeed, this follows since the kernels of the $2\times 3$ matrices are one dimensional, that is, they are directions in ${\mathbb{R}}^{3}$. We may therefore visualise the affine shapes as follows. First, we compute the intersection of the kernel direction with the half-sphere

[TABLE]

We then compute the stereographic projection from the point of coordinate $(0,0,1)$. This function is explicitly given by

[TABLE]

The resulting projected spline curve is plotted in Figure 8.

We stress that interpolation of affine shapes cannot, in general, be accomplished by standard spline interpolation in vector spaces. Indeed, if we introduce local coordinates by expressing the point $p_{2}=ae_{0}+be_{1}$ in the basis spanned by the basis vectors $e_{0}=p_{1}-p_{0}$ and $e_{1}=p_{3}-p_{0}$ with origin at $p_{0}$, then the $b$-coordinate is infinite as the spline passes through the affine shape (C). This is illustrated in Figure 9. More generally, any local coordinate chart would have similar artifacts for some shapes.

5 Outlook

There are a number of possible extensions of our methods and proofs to be studied in future work. Here we list some of them.

• •

We surmise that the conditions of Theorem 2.2 are in fact valid for more general symmetric spaces than Riemannian ones. In particular, we would like to be able to prove that result directly, without resorting to the previous result of Popiel and Noakes, which is the only reason for us to restrict to the Riemannian symmetric space case at this point.

• •

The only structure that we need to implement Algorithm 1 and to state Theorem 2.2 is an invariant connection (see (26a)), which only requires the homogeneous space at hand to be reductive [10, § 4]. We would expect similar results in that case, which would for instance cover arbitrary Lie groups, as well as Stiefel manifolds.

• •

A modification of our algorithm is possible in the more general case of arbitrary Riemannian manifolds (not necessarily symmetric), using the original $C^{2}$ condition of Popiel and Noakes [14]. We are working on an implementation and a proof of convergence.

• •

Our algorithm is currently based on a fixed point iteration. One could instead use a Newton iteration, which would ensure convergence in one step in the Euclidean case.

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 661482, from the Swedish Foundation for Strategic Research under grant agreement ICA12-0052, and from the Knut and Alice Wallenberg Foundation under grant agreement KAW-2014.0354.

Appendix A Identities and bounds on symmetric spaces

In this section we prove several identities and bounds on (Riemannian) symmetric spaces which are used in Appendix B and Appendix C below.

Let $M=G/H_{p}$ be a symmetric space, where $p\in M$ and $\mathfrak{g}=\mathfrak{m}_{p}\oplus\mathfrak{h}_{p}$ is the decomposition described in §​ 2.1, and let $\pi_{p}\vcentcolon\mathfrak{g}\to\mathfrak{m}_{p}$ be the canonical projection. Denote by $\Psi^{p}\vcentcolon G\to M$ and $\tilde{L}_{g}\vcentcolon M\to M$ the maps

[TABLE]

as usual, $L_{g}\vcentcolon G\to G$ and $R_{h}\vcentcolon G\to G$ are defined by

[TABLE]

and we identify $\mathfrak{g}=T_{e}G$ to define

[TABLE]

By an abuse of notation, we will use $g\cdot$ to denote any of the following action maps

[TABLE]

The notation $\langle\cdot,\cdot\rangle_{p}$ is used to denote the natural pairing of one-forms and tangent vectors over a point $p$.

The infinitesimal action (5) and the decomposition $\mathfrak{g}=\mathfrak{m}_{p}\oplus\mathfrak{h}_{p}$ define a principal connection [10], i.e., a $\mathfrak{g}$-valued one-form $\varpi$ on $M$ which fulfills

[TABLE]

for all $v\in T_{p}M$.

By construction, $\mathfrak{m}_{p}=\varpi(T_{p}M)$.

Throughout this section and the following, $\operatorname{dexp}$ denotes the trivialized derivative of the Lie group exponential, that is

[TABLE]

where the last equality follows from differentiating $\exp(\xi)\exp(-\xi)=\mathrm{id}$.

Lemma A.1**.**

Fix $p\in M$, and define $q\vcentcolon\mathfrak{g}\to M$ by $q(\xi)=\exp(\xi)\cdot p$. Then

[TABLE]

Proof.

By the properties of a Lie group action, for $\eta\in\mathfrak{g}$, we have

[TABLE]

We differentiate $q=\Psi^{p}\circ\exp$ by the product rule and get

[TABLE]

where the final equality is because $\pi_{p}(\eta)\cdot p=\eta\cdot p$ for all $\eta\in\mathfrak{g}$. ∎

Lemma A.2**.**

Fix $p\in M$. Let $\xi,\eta\in\mathfrak{m}_{p}$, then

[TABLE]

Proof.

As a function $\mathfrak{g}\to\mathfrak{g}$, $\operatorname{dexp}_{\xi}$ is an analytic function of $\operatorname{ad}_{\xi}$ with Taylor series $\operatorname{dexp}_{\xi}=\sum_{k=0}^{\infty}\frac{1}{(k+1)!}\operatorname{ad}_{\xi}^{k}$. Using the definition (4) of a symmetric space, we have, for any element $\eta\in\mathfrak{m}_{p}$, that

[TABLE]

The effect of $\pi_{p}$ is thus to eliminate terms in $\mathfrak{h}_{p}$, leaving only the terms with an even exponent $k$, so

[TABLE]

Note that, as ${\sinh(\operatorname{ad}_{\xi})}/{\operatorname{ad}_{\xi}}$ is an even function of $\xi$, we have

[TABLE]

∎

A corollary of Appendix A is that for $\xi$ near zero, $\pi_{p}\operatorname{dexp}_{\xi}$ is invertible as a function $\mathfrak{m}_{p}\to\mathfrak{m}_{p}$ with inverse given by $\operatorname{ad}_{\xi}/\sinh{\operatorname{ad}_{\xi}}$.

Lemma A.3**.**

Fix $p\in M$. Let $I_{p}$ be the symmetry function at $p\in M$, i.e., $I_{p}(q)=E_{p}(-\log_{p}(q)),$ and let $\xi\in\mathfrak{m}_{p}$, $\delta p\in T_{p}M$. Then

[TABLE]

Proof.

Let $q(\xi)=\exp(\xi)\cdot p$. Then $q\vcentcolon\mathfrak{m}_{p}\to M$ is a local diffeomorphism and $I_{p}(q(\xi))=q(-\xi)$. Differentiating both $q(\xi)$ and $I_{p}(q(\xi))$ with respect to $\xi$, and using Appendix A, we get

[TABLE]

Since $q$ is a diffeomorphism, we can choose $\eta\in\mathfrak{g}$ such that $T_{\xi}q(\eta)=\exp(\xi)\cdot w$, i.e., $w=\pi_{p}\operatorname{dexp}_{-\xi}(\eta)\cdot p$. Using Appendix A, we also have $w=\pi_{p}\operatorname{dexp}_{\xi}(\eta)\cdot p$, and the result follows. ∎

Lemma A.4**.**

Let $g\in G$, and $p,q\in M$. Then

[TABLE]

Proof.

Let $\sigma(t)=[g\cdot p,q]_{t}$ be the interpolating curve between $g\cdot p$ and $q$. Then $g^{-1}\cdot\sigma(t)=[p,g^{-1}\cdot q]_{t}$, and differentiating at $t=0$, we get

[TABLE]

∎

For the final lemmata, we assume $M=G/H_{p}$ to be a Riemannian symmetric space. This is equivalent to the existence of an $H_{p}$-invariant inner product on $\mathfrak{m}_{p}$, from which we derive a norm $\lVert\cdot\rVert$ on $\mathfrak{m}_{p}$.

Lemma A.5**.**

Let $\xi,\eta\in\mathfrak{m}_{p}$, then

[TABLE]

(Notice $\sin$, not $\sinh$, in the last denominator.)

Proof.

We use Appendix A to set $\pi_{p}\operatorname{dexp}_{\xi}=\frac{\sinh(\operatorname{ad}_{\xi})}{\operatorname{ad}_{\xi}}$. The lemma at hand then follows from considering the infinite series for $\frac{\sinh(\operatorname{ad}_{\xi})}{\operatorname{ad}_{\xi}}$ and its inverse, as well as the definition of the curvature $K$ in equation (22). ∎

Lemma A.6**.**

Fix $p\in M$. Let $\xi_{0},\xi_{1}\in\mathfrak{m}_{p}$, $q_{0}=\exp(\xi_{0})\cdot p$ and $q_{1}=\exp(\xi_{1})\cdot p$. Recall from §​ 3.1 that $d$ is the Riemannian distance between two points on $M$. Let us define

[TABLE]

We then have the inequality

[TABLE]

Proof.

Let $\gamma(t)$ be the geodesic from $\gamma(0)=q_{0}$ to $\gamma(1)=q_{1}$, and let $\alpha(t)\in\mathfrak{m}_{p}$ be such that

[TABLE]

Then $\xi_{0}=\alpha(0)$, $\xi_{1}=\alpha(1)$. We note that, by the triangle inequality, we have

[TABLE]

Differentiating (27) we get

[TABLE]

which we rewrite as

[TABLE]

where we have used that $\pi_{p}(\eta)\cdot p=\eta\cdot p$ for all $\eta\in\mathfrak{g}$. Taking norms on both sides and using that $G$ acts by isometries, we have

[TABLE]

Moreover, by Appendix A,

[TABLE]

We now use the monotonicity of $x/\sin x$, and the bound (28) to obtain

[TABLE]

∎

Appendix B $C^{2}$ continuity

We now prove Theorem 2.2. The proof relies on a result by Popiel and Noakes [14].

Proof of Theorem 2.2.

At non-integer $t$, the spline is $C^{\infty}$. We show that the spline is $C^{2}$ at $t=i\in\{\,1,2,\dotsc,N-1\,\}$, i.e., at the interpolating points $p_{i}$. By [14, Thm. 3], a sufficient and necessary condition for the spline to be $C^{2}$ at $t=i$ is

[TABLE]

where $I_{p_{i}}$ is the symmetry function at $p_{i}$. We proceed by removing the dependence of $TI_{p_{i}}$ in (29) by using the lemmata from Appendix A.

Recall that $q_{i}^{+}=\exp(\xi_{i})\cdot p_{i}$ and $q_{i}^{-}=\exp(-\xi_{i})\cdot p_{i}$. By Appendix A, we have

[TABLE]

Using Appendix A, with $w=\log_{p_{i}}(\exp(-\xi_{i})\cdot q_{i+1}^{-})$, we can write

[TABLE]

Next, we note that

[TABLE]

Using (30) and (31), we can rewrite (29) as

[TABLE]

Acting on the whole equation from the left with $\exp(\xi_{i})$ and rearranging, we get

[TABLE]

We now obtain the $C^{2}$ condition (13) by using the $C^{1}$ condition (12) to replace $v_{i}=\log_{p_{i}}(q_{i}^{+})=-\log_{p_{i}}(q_{i}^{-})$ in the equation above. ∎

Appendix C Convergence of the fixed-point method

We now proceed to prove Theorem 3.1, i.e., that if consecutive interpolation points are sufficiently close then Algorithm 1 converges to a solution. In the proof, we will work with variables in the Lie algebra $\mathfrak{g}$. Let

[TABLE]

and let

[TABLE]

be the Lie algebra elements defined by $\xi_{i}=\langle\varpi,v_{i}\rangle_{p_{i}}.$

Switching to variables in the Lie algebra, the fixed-point iteration (17) becomes a map

[TABLE]

For $1\leq i\leq N-1$, $\phi_{i}(\bm{\xi})=\frac{1}{4}\langle\varpi,w_{i}(\bm{\xi})\rangle_{p_{i}}+\frac{1}{2}\xi_{i}$. To simplify notation, we define $L_{p}\vcentcolon M\to\mathfrak{m}_{p}$ by

[TABLE]

Now, for $0<i<N$ the fixed-point iteration map $\phi_{i}$ can be expressed as

[TABLE]

and both $\phi_{0}$ and $\phi_{N}$ are given by the boundary conditions

[TABLE]

Specifically, for clamped boundary conditions, we get

[TABLE]

where $\xi_{\text{start}}$ and $\xi_{\text{end}}$ are prescribed. For free boundary conditions

[TABLE]

It is convenient to introduce the auxiliary variables

[TABLE]

and to write $\phi_{i}(\bm{\xi})=\phi_{i}^{+}(\bm{\xi})+\phi_{i}^{-}(\bm{\xi}),$ where

[TABLE]

Furthermore, we define the norms

[TABLE]

and

[TABLE]

Recall that $D$ denotes the maximum distance between neighbouring interpolation points as declared in equation (21). Notice that by the triangle inequality, we have

[TABLE]

and similarly for $\lVert\omega_{i}^{-}(\bm{\xi})\rVert$, so

[TABLE]

We are now in position to prove Theorem 3.1. The proof consists in showing the existence of an invariant region, and thereafter bounding the partial derivatives of $\bm{\phi}$ to show that it is a contraction.

C.1 Invariant region

We begin by establishing an invariant region.

Proposition C.1**.**

Let $V>0$ satisfy the inequality

[TABLE]

where $U=D+2V$. If either one of the following boundary conditions is fulfilled:

1. (i)

clamped spline boundary conditions are used, $\lVert\xi_{\text{start}}\rVert\leq V$, $\lVert\xi_{\text{end}}\rVert\leq V$ and $\lVert\bm{\xi}\rVert\leq V$, 2. (ii)

natural spline boundary conditions are used and $\lVert\bm{\xi}\rVert\leq V$,

then we also have $\lVert\bm{\phi}(\bm{\xi})\rVert\leq V$.

To reach this result, we first prove a few bounds for $\phi_{i}^{+}(\bm{\xi})$.

Lemma C.2**.**

Let $1\leq i\leq N-1$, then the following bounds hold,

[TABLE]

where $U_{i}=d(p_{i},p_{i+1})+\lVert\xi_{i+1}\rVert+\lVert\xi_{i}\rVert$.

Proof.

We prove the first inequality. The proof of the second inequality is entirely symmetrical and therefore omitted.

By definition (35),

[TABLE]

The proof uses Appendix A to bound $\omega_{i}^{+}(\bm{\xi})+\xi_{i}=\omega_{i}^{+}(\bm{\xi})-(-\xi_{i})$.

Recall that $G$ acts isometrically on $\mathfrak{m}_{p}$ by the adjoint action

[TABLE]

where $\cdot\omega_{i}^{+}(\bm{\xi})$ and $\xi_{i}$ are vectors in $\mathfrak{m}_{q_{i}^{+}}$. We also have

[TABLE]

and

[TABLE]

By Appendix A, we have

[TABLE]

where

[TABLE]

From the triangle inequality, we have

[TABLE]

Thus $T\leq d(p_{i},p_{i+1})+\lVert\xi_{i}\rVert+\lVert\xi_{i+1}\rVert=U_{i}$, and the claim follows by monotonicity of the function $\frac{x}{\sin x}$. ∎

Proof of §​ C.1.

Using §​ C.1, and the inequalities $\lVert\xi_{i}\rVert\leq\lVert\bm{\xi}\rVert\leq V,$ $d(p_{i},p_{i+1})\leq D$ and $U_{i}\leq U$, we get that

[TABLE]

for all $1\leq i\leq N-1$.

For the boundary velocities:

1. (i)

Clamped spline: $\phi_{0}(\bm{\xi})=\xi_{\text{start}}$, so $\lVert\phi_{0}(\xi)\rVert\leq V$ by our assumption. 2. (ii)

Natural spline: $\phi_{0}(\bm{\xi})=\frac{1}{2}\langle\varpi,\log_{p_{0}}(\exp(-\xi_{1})\cdot p_{1})\rangle_{p_{0}}$, so

[TABLE]

under a weaker condition on $V$. Similarly for $\phi_{N}(\bm{\xi})$.

In conclusion, under the assumption (38), we have

[TABLE]

and \big{\{}\,\bm{\xi}\;\big{|}\;\lVert\bm{\xi}\rVert\leq V\,\big{\}} forms an invariant region. ∎

C.2 Contraction

We now establish that, for $D$ and $\lVert\bm{\xi}\rVert$ sufficiently small, there exists an $\alpha\in(0,1)$ such that, in the operator norm derived from (36) $\lVert T_{\bm{\xi}}\bm{\phi}\rVert\leq\alpha$. We first consider the effect of varying $\xi_{i+1}$ and $\xi_{i-1}$ in $\phi_{i}(\bm{\xi})$.

Proposition C.3**.**

Let $1\leq i\leq N-1$ and denote by $\frac{\partial\phi_{i}(\bm{\xi})}{\partial\xi_{j}}$ the partial derivative of $\phi_{i}$ with respect to $\xi_{j}$, i.e. $\frac{\partial\phi_{i}(\bm{\xi})}{\partial\xi_{j}}$ is a linear operator $\mathfrak{m}_{p_{j}}\to\mathfrak{m}_{p_{i}}$. Then, in the operator norm,

[TABLE]

and

[TABLE]

Proof.

We only prove the first claim. The proof of the second is entirely symmetrical. The term $\xi_{i+1}$ only enters $\phi_{i}(\bm{\xi})$ through the term $\frac{1}{4}\omega_{i}^{+}(\bm{\xi})$, so

[TABLE]

By the definition of $\omega_{i}^{+}(\bm{\xi})$, we have

[TABLE]

Differentiating implicitly on both sides, then acting from the left by $\exp(-\omega_{i}^{+}(\bm{\xi}))$, we get

[TABLE]

where $\delta\omega_{i}^{+}=\frac{\partial\omega_{i}^{+}(\bm{\xi})}{\partial\xi_{i+1}}\delta\xi_{i+1}$. Taking norms on both sides, and using that the group action is isometric, we get

[TABLE]

and

[TABLE]

It follows that

[TABLE]

thus

[TABLE]

The claims in the proposition follow from (39) and Appendix A. ∎

We note that we can also use a Taylor expansion of the right hand side to get the asymptotic bound

[TABLE]

To bound $\frac{\partial\phi_{i}(\bm{\xi})}{\partial\xi_{i}}$, it is again advantageous to do the splitting $\phi_{i}(\bm{\xi})=\phi_{i}^{+}(\bm{\xi})+\phi_{i}^{-}(\bm{\xi})$, and consider each term separately.

While it is possible to bound $\left\lVert\frac{\partial\phi_{i}^{+}(\bm{\xi})}{\partial\xi_{i}}\right\rVert$ by trigonometric and hyperbolic functions of $K\lVert\xi_{i}\rVert$ and $K\lVert\omega_{i}^{+}(\bm{\xi})\rVert$, the required calculations are lengthy, and we restrict ourself to an asymptotic bound.

Proposition C.4**.**

Let $1\leq i\leq N-1$ and denote by $\frac{\partial\phi_{i}^{+}(\bm{\xi})}{\partial\xi_{i}}\vcentcolon\mathfrak{m}_{p_{i}}\to\mathfrak{m}_{p_{i}}$ the partial derivative of $\phi_{i}$ with respect to $\xi_{i}$. In the operator norm

[TABLE]

and

[TABLE]

Again, we will only prove the first claim, the proof of the second is entirely symmetrical. We first prove the following Lemma on the derivative of $\omega_{i}^{+}$.

Lemma C.5**.**

For $1\leq i\leq N-1$, the following equality holds

[TABLE]

Proof.

We use that

[TABLE]

By differentiating implicitly, we get

[TABLE]

where the second equality uses (40). Acting from the left with $\exp(-\omega_{i}^{+}(\bm{\xi}))$, we get

[TABLE]

Applying the connection $\varpi$ to both sides, we get

[TABLE]

where the bracketed expressions are linear operators $\mathfrak{m}_{p_{i}}\to\mathfrak{m}_{p_{i}}$. Now, the map ${\pi_{p_{i}}\operatorname{dexp}_{-\omega_{i}^{+}(\bm{\xi})}}\vcentcolon\mathfrak{m}_{p_{i}}\to\mathfrak{m}_{p_{i}}$ is invertible by Appendix A, and the claim follows. ∎

We now proceed with the proof of §​ C.2.

Proof of §​ C.2.

From the definition of $\phi_{i}^{+}$,

[TABLE]

Using §​ C.2, we have

[TABLE]

Using Appendix A, the series expansions for $\operatorname{Ad}_{\exp(-\omega_{i}^{+}(\bm{\xi}))}=\exp(-\operatorname{ad}_{\omega_{i}^{+}})$, and $\operatorname{dexp}_{-\xi_{i}}$, the expression on the right hand side can be expanded as an infinite series in $\operatorname{ad}_{\omega_{i}^{+}(\bm{\xi})}$ and $\operatorname{ad}_{\xi_{i}}$ applied to $\delta\xi_{i}$. (Recall that, as we saw in the proof of Appendix A, the effect of $\pi_{p_{i}}$ is to cancel out Lie monomials of even degree.) The crucial point is that there is a cancellation in the leading term so that, neglecting terms of fourth order and higher,

[TABLE]

We thus get

[TABLE]

The proposition follows by taking norms. ∎

It is also possible to bound the derivatives of the boundary condition functions.

Proposition C.6**.**

**

1. (i)

(Clamped splines) Let $\phi_{0}$ and $\phi_{N}$ be as in (32). Then $\frac{\partial\phi_{0}(\bm{\xi)})}{\partial\xi_{1}}=0$, $\frac{\partial\phi_{N}(\bm{\xi})}{\partial\xi_{N-1}}=0$. 2. (ii)

(Natural splines) Let $\phi_{0}$ and $\phi_{N}$ be as on (33). Then

[TABLE]

and similar for $\left\lVert\frac{\partial\phi_{N}(\bm{\xi})}{\partial\xi_{N-1}}\right\rVert$.

Proof.

1. (i)

For clamped splines, $\phi_{0}$ and $\phi_{N}$ are constant. 2. (ii)

For natural splines, we have

[TABLE]

Differentiating and taking norms, we get

[TABLE]

where $\delta\phi_{0}=\frac{\partial\phi_{0}(\bm{\xi})}{\partial\xi_{1}}\delta\xi_{1}$. Appendix A now gives

[TABLE]

and the claim follows by a Taylor expansion of $\frac{x}{\sin(x)}$ and $\frac{\sinh(x)}{x}$. The corresponding bound and proof for $\left\lVert\frac{\partial\phi_{N}(\bm{\xi})}{\partial\xi_{N-1}}\right\rVert$ is entirely symmetrical.

∎

C.3 Convergence

We are now ready to prove the convergence theorem.

Proof of Theorem 3.1.

A consequence of §​ C.1 is that $\bm{\phi}$ has an invariant region for small enough $D$. We need to also establish that it is a contraction. It is sufficient to show that there exists an $\alpha\in(0,1)$ such that $\lVert T_{\bm{\xi}}\bm{\phi}\rVert\leq\alpha$ in an invariant region.

We have

[TABLE]

as well as

[TABLE]

By §​ C.2, §​ C.2, and §​ C.2, we can bound the terms appearing with a Taylor expansion to obtain

[TABLE]

We can now use (37) and §​ C.2 to obtain a bound of the form

[TABLE]

It is therefore clear that when $K\max(D,\lVert\bm{\xi}\rVert)$ is sufficiently small, $\bm{\phi}$ is a contraction. When $D$ approaches zero, the smallest $V$ satisfying (38) also approaches zero. Therefore, when $D$ is small enough, there is a $V$ such that

• •

\Big{\{}\,\bm{\xi}\;\Big{|}\;\lVert\bm{\xi}\rVert\leq V\,\Big{\}} is a compact invariant region and

• •

$\bm{\phi}$ is a contraction on that region.

Therefore, the fixed point iteration converges to a solution. ∎

See Figure 10 for an illustration of the final argument in the proof.

Appendix D Python implementation

Here we briefly indicate how to use our Python implementation to compute splines in various geometries. First, one has to install the following package by following the installation instructions:

https://github.com/olivierverdier/bsplinelab

The minimal necessary import is

Here is an example of how to interpolate between random points in the Euclidean case:

⬇

interpolation_points = np.random.rand(10,4)

b = cubic_spline(Symmetric, interpolation_points)

The result is then a $C^{2}$ function b which is equal to the interpolation points at integer points on the interval $[0,9]$.

In order to interpolate on, for instance, a sphere, one would do as follows:

⬇

3 points on the sphere:

interpolation_points = np.array([[1.,0,0], [0,1,0], [0,0,1]])

b = cubic_spline(Symmetric, interpolation_points, geometry=Sphere())

The function b is now a function defined on the interval $[0,2]$ such that it is $C^{2}$, is interpolating the prescribed points at the integer points [math], $1$ and $2$, and is on the sphere for all points in between.

For further information on bsplinelab, we refer to the package documentation and the available code examples.