A Numerical Algorithm for C2-splines on Symmetric Spaces
AA Numerical Algorithm for C -splineson Symmetric Spaces Klas Modin ∗ , Geir Bogfjellmo † , and Olivier Verdier ‡ Mathematical Sciences, Chalmers and University of Gothenburg, Sweden Instituto de Ciencias Matem´atica, Madrid, Spain Department of Computing, Mathematics and Physics, Western Norway University of AppliedSciences, Bergen, Norway
Cubic spline interpolation on Euclidean space is a standard topic in numer-ical analysis, with countless applications in science and technology. In severalemerging fields, for example computer vision and quantum control, there isa growing need for spline interpolation on curved, non-Euclidean space. Thegeneralization of cubic splines to manifolds is not self-evident, with severaldistinct approaches. One possibility is to mimic the acceleration minimizingproperty, which leads to
Riemannian cubics . This, however, requires the so-lution of a coupled set of non-linear boundary value problems that cannotbe integrated explicitly, even if formulae for geodesics are available. Anotherpossibility is to mimic De Casteljau’s algorithm, which leads to generalizedB´ezier curves . To construct C -splines from such curves is a complicatednon-linear problem, until now lacking numerical methods. Here we providean iterative algorithm for C -splines on Riemannian symmetric spaces, andwe prove convergence of linear order. In terms of numerical tractability andcomputational efficiency, the new method surpasses those based on Rieman-nian cubics. Each iteration is parallel, thus suitable for multi-core implemen-tation. We demonstrate the algorithm for three geometries of interest: the n -sphere, complex projective space, and the real Grassmannian. MSC-2010:
Keywords:
De Casteljau, cubic spline, Riemannian symmetric space,B´ezier curve ∗ [email protected] † [email protected] ‡ [email protected] a r X i v : . [ m a t h . NA ] M a r ontents
1. Introduction 22. Generalized B ´ezier curves on symmetric spaces 43. Numerical algorithm 84. Examples 115. Outlook 14A. Identities and bounds on symmetric spaces 16B. C continuity 20C. Convergence of the fixed-point method 20D. Python implementation 28
1. Introduction
We address the following.
Problem 1.1.
Given a set of points { p i } Ni =0 on a connected manifold M , construct a C path γ : [0 , N ] → M such that γ ( i ) = p i .For M = R n every textbook in numerical analysis teaches cubic splines , i.e., piecewisepolynomials of order 3 with matching first and second derivatives. However, the gener-alization of cubic splines to curved space is non-trivial, essentially because polynomialsare 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 complexprojective space C P n (see § k, n ) (see § M is a Riemannian manifold, a natural generalization of cubic splines are Rie-mannian cubics [12, 2]. They are based on minimizing acceleration under interpolationconstraints, by solving the problemmin γ ( i )= p i (cid:90) N |∇ ˙ γ ( t ) ˙ γ ( t ) | d t. (1)The corresponding Euler–Lagrange equation is a fourth order ODE on M with com-plicated, coupled boundary conditions [12]. In addition to the Euclidean case, explicitsolutions to the initial-value problem are known for so called null Lie quadratics on2O(3) and SO(2 ,
1) equipped with the bi-invariant Riemannian metric [11]. For othercases, even when geodesics are explicitly known, one is left with numerical ODE meth-ods, combined, for example, with a shooting method to match the boundary conditions.This approach is computationally costly and the associated convergence analysis is in-volved. Another possibility is to only approximatelly fulfill the interpolation constraintsby incorporating them in the minimization functional and then use a steepest-descentmethod 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´ezier curves (GBC). They are based on adirect generalization of De Casteljau’s algorithm [13]. Splines can then be constructedby gluing piecewise GBC with matching boundary conditions. In this way, a fully explicitalgorithm for C -splines on Grassmannian and Stiefel manifolds has been developed [8].The C condition, however, first derived by Popiel and Noakes [14], is significantlymore complicated than the C condition. For this reason, there is a lack of numericalalgorithms for C -splines. Nevertheless, C continuity is often required in applications,to avoid discontinuities in accelerations.In this paper we give the first numerical algorithm for C -splines based on GBC.Instead of treating all Riemannian manifolds, we focus on the subclass of symmetricspaces (see § C condition of Popiel and Noakes is, however, still complicated. We nowlist the specific contributions of our paper.1. Using the special homogeneous space structure of symmetric spaces, we give asignificant simplification of the C condition (Theorem 2.2).2. Using the new C condition, we construct an algorithm based on fixed-point iter-ations (Algorithm 1).3. We prove convergence of the algorithm under conditions on the maximal distancebetween consecutive interpolation points (Theorem 3.1).The paper is organized as follows. In § C condition. In § §
4. An outlooktowards future work is given in §
5. The proofs of the C condition and the convergenceresult are long and technical, and therefore given in Appendix A-C. In Appendix D wegive a brief demonstration of an easy-to-use open-source Python package implementingour interpolation algorithm for several geometries.3
G H
ImplementationEuclidean space R n SE( n ) SO( n ) Flat n -sphere S n O( n + 1) O( n ) Sphere
Hyperbolic space H n O( n,
1) O( n ) Hyperbolic
Real Grassmannian Gr( k, n ) O( n ) O( k ) × O( n − k ) Grassmannian
Orthogonal group O( n ) O( n ) × O( n ) O( n )Complex projective space C P n U( n + 1) U( n ) × U(1)
Projective
Complex Grassmannian Gr C ( k, n ) U( n ) U( k ) × U( n − k )Unitary group U( n ) U( n ) × U( n ) U( n )Co-variance matrices P ( n ) GL( n ) O( n ) Table 1: Some examples of Riemannian symmetric spaces. For the details of the cor-responding symmetric space decomposition, see [10, § Python package; we refer to Appendix D for example usage.
2. Generalized B ´ezier curves on symmetric spaces
Let us describe the construction of generalized B´ezier curves on a symmetric space. Theconstruction relies on a notion of interpolation between two points on the manifold,generalizing the notion of a geodesic on a Riemannian manifold [9].
In this section we briefly recall symmetric spaces. A list of the most common examplesof symmetric spaces is given in Table 1. For details, we refer to [16] or [4]. We assumethat a Lie group G acts transitively on M , making M a homogeneous space [16, § action by G × M (cid:51) ( g, p ) (cid:55)→ g · p ∈ M. (2)For every point p ∈ M , we define the isotropy subgroup H p at pH p := { g ∈ G | g · p = p } . Now, the homogeneous space M is called a symmetric space if the Lie algebra g of G has a direct sum decomposition g = h p ⊕ m p , (3)where h p is the Lie algebra of H p , and the summands fulfill the following algebraicconditions [ h p , h p ] ⊂ h p , [ h p , m p ] ⊂ m p , [ m p , m p ] ⊂ h p . (4)If these conditions are satisfied at one point p ∈ M , then they are automatically validat every point with h g · p = g · h p , m g · p = g · m p , where we follow the convention formulated in Appendix A that g · denotes the canonicalaction on the space at hand (here, the adjoint action of a Lie group on its Lie algebra).4 i p i +1 q + i q − i +1 t = 0 . t = 0 .
25. The algorithm proceeds iteratively by choosing points atposition 0 .
25 on each geodesic as indicated in the figure. The precise definitionof the resulting curve is given in (8).By differentiating the action map ( g, p ) (cid:55)→ g · p with respect to the arguments, we alsoobtain, on the one hand, an infinitesimal action of g on M , and, on the other hand, a prolonged action of G on T M . For ξ ∈ g and p ∈ M , we denote the infinitesimal actionby ( ξ, p ) (cid:55)→ ξ · p ∈ T p M. (5)For a fixed p , this map induces an isomorphism between m p and T p M . We now define the interpolating curve η : [0 , → M between two (close enough) points p and p as the unique curve satisfying η ( t ) = exp( tξ ) · p η (1) = p ξ ∈ m p . (6)Let us introduce the notation [ p , p ] t := η ( t ) . In particular, notice that [ p , p ] = p , [ p , p ] = p . (7)In many (but not all) cases, a symmetric space M can be equipped with the structureof a Riemannian manifold, given through an H p -invariant inner product (cid:104)· , ·(cid:105) p on m p . Inthis case, the interpolating curves are geodesics. From here on, we shall assume that M is such a Riemannian symmetric space . We briefly describe the De Casteljau construction, illustrated on Figure 1. This algorithmconstructs B´ezier curves from interpolating curves. We focus on cubic
B´ezier curves.5 q = E p ( v ) · pv = log p ( q )Figure 2: Illustration of the movement function E p and of the logarithm described in § p and q is a geodesic.From four points p , q , q and p , one defines the cubic B´ezier curve γ : [0 , → M by γ ( t ) := (cid:104)(cid:2) [ p , q ] t , [ q , q ] t (cid:3) t , (cid:2) [ q , q ] t , [ q , p ] t (cid:3) t (cid:105) t . (8)Note that, by construction, γ (0) = p and γ (1) = p . For the C condition in § g ∈ G on a point x ∈ M : g · x ∈ M. (9)2. A movement function E , which, to a point p and a velocity v at p , assigns anelement E p ( v ) ∈ G . This function should generate interpolating curves (6), i.e., itshould be of the form E p ( ξ · p ) = exp( ξ ) ξ ∈ m p . (10)3. The logarithm ; given two points p and q in M , it is defined by v ∈ T p M, E p ( v ) · p = q, ⇐⇒ v = log p ( q ) . (11)Throughout the paper we make the blanket assumption that p and q are neverpast conjugate points, so that the logarithm is always well-defined. C condition We shall now consider curves consisting of piecewise cubic B´ezier curves. As we haveseen in § C continuity at the6 i − q + i − p i q − i +1 p i +1 q − i q + i exp( ξ i ) · q + i − exp( − ξ i ) · q − i +1 ξ i − ξ i ξ i − ξ i Figure 3: An illustration of the C condition in Theorem 2.2. The sum of the two vectorsat p i generating the dashed green geodesics should equal the sum of the twovectors at p i generating the dotted blue geodesics. Notice that the action ofexp( ξ i ) on q + i − and q − i +1 does not generate geodesics from these points, since ξ i does not belong to m q + i − or m q − i +1 . For this reason, the condition (13) inTheorem 2.2 does not work for arbitrary Riemannian manifolds.interpolating points, which leads to a linear set of equations. On a general Riemannianmanifold, the corresponding C condition, given by Popiel and Noakes [14], is highlynonlinear, involving the inverse of the derivative of the Riemannian exponential. In thespecial case of symmetric spaces, however, we now show that the C condition simplifiessignificantly, involving only the three operations (9), (10), and (11). Definition 2.1.
Let p = ( p , . . . , p N ) be points on a symmetric space M . A compositecubic B´ezier curve on M is a curve γ : [0 , N ] → M such that for each integer 0 ≤ i ≤ N , γ ( i ) = p i , and such that for each integer 0 ≤ i < N , the restriction γ i := γ | [ i,i +1] is acubic B´ezier curve (as defined in § γ i , we associate two controlpoints , denoted q + i and q − i +1 . A C -continuous composite cubic B´ezier curve is called a cubic spline .The interpolating conditions ensure that γ is continuous (or C ). However, for higherdegree of continuity ( C or C ), 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 γ : [0 , N ] → M be acomposite cubic B´ezier curve according to Definition 2.1. Consider the conditions onthe control points given by log p i ( q + i ) = − log p i ( q − i ) =: v i (12)log p i ( E p i ( − v i ) · q − i +1 ) − log p i ( q + i ) = log p i ( E p i ( v i ) · q + i − ) − log p i ( q − i ) . (13)(i) If condition (12) is fulfilled, then γ is a C -curve. If conditions (12) and (13) are fulfilled, then γ is a C -curve, i.e., a cubic spline.Furthermore, if γ fulfills the C conditions (12) and (13) , then it is uniquely determinedeither by clamped boundary conditions, where γ (cid:48) (0) and γ (cid:48) ( N ) are prescribed, or bynatural boundary conditions, where ∇ γ (cid:48) (0) γ (cid:48) (0) = ∇ γ (cid:48) (0) γ (cid:48) (0) = 0 . Remark 2.3.
Notice how (12) and (13) are generalizations of the first and second orderfinite differences in Euclidean geometry. Indeed, in this case condition (12) reads q + i − p i = − ( q − i − p i ) (14)and condition (13) reads q + i − i − q − i + p i = q − i +1 − q + i + p i . (15)See Figure 3 for a geometric illustration of the C condition (13). Outline of proof.
The details of the proof are rather technical and therefore given inAppendix B. However, for convenience we also provide an outline here:(i) The C condition in the case of Riemannian symmetric spaces is given by Popieland Noakes [14, Sec. 4], expressed in terms of the derivative of the symmetryfunction I p ( q ) = E p ( − log p ( q )).(ii) Using the special symmetric space structure, we show that the derivative of I p canbe expressed solely in terms of the movement function E p (Lemma A.3).(iii) Together with a result on equivariance of the log function (Lemma A.4) we canthen simplify the original Riemannian symmetric space C condition by Popiel andNoakes [14] to the one given by (13).
3. Numerical algorithm
We shall now construct a fixed point iteration algorithm for cubic splines on symmetricspaces, based on the C condition in Theorem 2.2. To this end, we parameterize the con-trol points q − i and q + i by the corresponding tangent vectors v i as defined in Theorem 2.2.Next, let w i ( v i − , v i , v i +1 ) :=log p i (cid:0) E p i ( − v i ) E p i +1 ( − v i +1 ) · p i +1 (cid:1) − log p i (cid:0) E p i ( v i ) E p i − ( v i − ) · p i − (cid:1) . (16)The iteration map over v := ( v , . . . , v N − ) is then given by v (cid:55)→ w ( v ) + 12 v , (17)8 i p i +1 q + i q − i +1 v i v i +1 Figure 4: An illustration of the control points and the corresponding velocities. Wedefine a velocity v i at every point p i . From that velocity, we construct thecontrol point q + i := E p i ( v i ) · p i . Similarly, we use v i +1 to construct the point q − i +1 := E p i +1 ( − v i +1 ) · p i +1 . We can then construct a cubic B´ezier curve fromthe four points p i , q + i , q − i +1 , p i +1 as indicated in Figure 1. If the velocities v i are chosen following Theorem 2.2, the resulting piecewise cubic curve will haveoptimal regularity.where w ( v ) := (cid:0) w ( v , v , v ) , . . . , w N − ( v N − , v N − , v N ) (cid:1) . (18)Pseudo code for the complete algorithm is given in Algorithm 1. We consider two optionsfor boundary conditions.(i) Clamped spline: v ( v ) = v start , v N ( v ) = v end (19)where v start and v end are prescribed constants.(ii) Natural spline: v ( v ) = 12 log p ( q − ) , v N ( v ) = −
12 log p N ( q + N − ) . (20) In this section we give a result on convergence of the fixed point method given byAlgorithm 1. To do so, we first give some preliminary definition.The Riemannian distance between p, q ∈ M is denoted by d ( p, q ). Our proof ofconvergence uses that two consecutive points p i , p i +1 are close. Thus, we define theconstant D = max ≤ i ≤ N − d ( p i , p i +1 ) . (21)The other constant that comes into play is a bound on (essentially) the curvature of M .Indeed, if (cid:107)·(cid:107) p denotes the H p -invariant norm on m p associated with the Riemannian9 lgorithm 1 Computing the control points q + i and q − i The auxiliary functions E and log are defined in (10) and (11). v ← (0 , . . . , repeat ¯ v ← ( v ( v ) , v , v N ( v )) (cid:46) v and v N are determined by (19) or (20) for i ← , . . . , N − do g + i ← E p i (¯ v i ) q + i ← g + i · p i end forfor i ← , . . . , N do g − i ← E p i ( − ¯ v i ) q − i ← g − i · p i end forfor i ← , . . . , N − do δ i ← log p i ( g − i · q − i +1 ) − log p i ( g + i · q + i − ) − v i v i ← v i + δ i end foruntil (cid:107) δ (cid:107) ≤ T OL
The outcome is the control points q + i and q − i , from which one can compute the splinesegments.structure of M , then we define the constant K ≥ K = sup ξ,η,ζ ∈ m p \ (cid:107) [ ξ, [ η, ζ ]] (cid:107)(cid:107) ξ (cid:107)(cid:107) η (cid:107)(cid:107) ζ (cid:107) . (22)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 boundaryconditions converges linearly. If KD and max { K (cid:107) v (cid:107) , K (cid:107) v N (cid:107) } are both sufficientlysmall, then Algorithm 1 with clamped boundary conditions converges linearly.Outline of proof. Again, the proof is long and technical, and therefore given in Ap-pendix C. Here we give a brief outline.The proof is based on showing that the iteration mapping (17) is a contraction. Forconvenience, we use the variables ξ i ∈ m p i instead of v i ∈ T p i M . The iteration mappingis denoted φ .(i) The first step is to give an invariant region of the iteration mapping. This is givenin Proposition C.1, which proves existence of V > D and K ), such that (cid:107) ξ (cid:107) ≤ V implies (cid:107) φ ( ξ ) (cid:107) ≤ V .(ii) The next step is to prove that φ is a contraction, which is established by a seriesof estimates on the partial derivatives on φ (Proposition C.4, Proposition C.3, andProposition C.6). These estimates depend on K .10iii) The final step, in § C.3, is to use the contraction mapping theorem.
4. Examples
Unit quaternions (versors) are used extensively in computer graphics to represent 3-Drotations. They are elements of the form q = q + q i + q j + q k with q + q + q + q = 1 . (23)Thus, we can identify unit quaternions with S . In turn, S is a double cover of SO(3).A point on S therefore gives a rotation matrix, and likewise any rotation matrix cor-responds to two antipodal points on S . The rotation of a vector p = ( p x , p y , p z ) by aunit quaternion q can be compactly written using quaternion multiplication asrot q ( p ) = qpq − where p = p x i + p y j + p z k is thought of as a pure imaginary quaternion.Since S (cid:39) O(4) / O(3) is a Riemannian symmetric space (with respect to the standardRiemannian metric), we can use Algorithm 1 to obtain C continuous spline interpolationbetween rotations. Since geodesics on S are given by great circles, it is straightforwardto derive the mappings E and log.The resulting C -curve interpolating 5 orientations is shown in Figure 5. In the figure,an element γ ( t ) ∈ S is represented by the rotated basis vectors { rot γ ( t ) ( i ) , rot γ ( t ) ( j ) , rot γ ( t ) ( k ) } . The actual interpolation points are marked with bolder lines.Note that S can be canonically identified with SU(2). Note also that SU(2) has asymmetric space structure as the quotient (SU(2) × SU(2)) / SU(2) [10]. Here, however,we considered the symmetric space structure S = SO(4) / SO(3). These two symmetricspace structure are locally equivalent. This is because, SU(2) × SU(2) is a double cover ofSO(4) and S is a double cover of SO(3). Our algorithm gives the same result regardlessof the choice of either of those two symmetric space structures on S . From this pointof view, our algorithm could be used for spline interpolation on any compact Lie group G , using the Cartan–Schouten symmetric space structure ( G × G ) /G . The control of quantum states is an important subproblem in quantum information andthe realization of quantum computers [3]. The objective is to find a time-dependentHamiltonian H ( t ) designed to ‘steer’ a given quantum state | ψ (cid:105) through a sequence ofstates | ψ (cid:105) , | ψ (cid:105) , . . . , | ψ N (cid:105) at given times t , . . . , t N . As instantaneous switching of the11igure 5: Interpolation between 5 orientations, illustrated by the orientation of the frameconsisting of the axes i (red), j (green), and k (blue). The bold frames markthe interpolants.Hamiltonian is not experimentally feasible [1], the interpolating curve | ψ ( t ) (cid:105) should beat least C continuous.If the quantum state space is finite dimensional, corresponding to an ensemble ofqubits, then, in the geometric description of quantum mechanics [6, 7], the phase spaceis given by complex projective space C P n . The quantum control problem can then beseen as a two-step process:1. Find an interpolating curve t (cid:55)→ | ψ ( t ) (cid:105) ∈ C P n such that | ψ (0) (cid:105) = | ψ (cid:105) and | ψ ( t k ) (cid:105) = | ψ k (cid:105) .2. Using the homogeneous structure C P n (cid:39) U( n + 1) / U( n ) × U(1), lift | ψ ( t ) (cid:105) toa curve t (cid:55)→ g ( t ) ∈ U( n + 1) such that g (0) = e and π ( g ( t )) = | ψ ( t ) (cid:105) . Thetime-dependent Hamiltonian is then given by H ( t ) = − i ˙ g ( t ) g ( t ) − .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 C P , which isisomorphic to a sphere (called the Bloch sphere ). Using 6 interpolation points, theresulting interpolating curve is visualized on the Bloch sphere in Figure 6.
The space of affine shapes is the quotient set whose elements consist of k + 1 labelled( i.e. ordered) points in R n defined up to affine transformations. Thus, ( p , . . . , p k ) ∈ ( R n ) k +1 and (˜ p , . . . , ˜ p k ) ∈ ( R n ) k +1 define the same affine shape if there exists an affinetransformation p (cid:55)→ Ap + q such that (˜ p , . . . , ˜ p k ) = ( Ap + q, . . . , Ap k + q ). The studyof affine shapes has many applications in computer vision [5]. An example is the identi-fication of planar surfaces (such as facades of buildings) in images taken from differentprojection angles.The space of affine shapes with full rank (meaning that the k + 1 points are notcontained in an affine hyperplane of R n ) can be identified with the real Grassmannianmanifold Gr( n, k ) of n -dimensional subspaces in R k . Indeed, choose a coordinate systemand place the coordinates of the vectors p − p , . . . , p k − p column by column in an n × k matrix. The kernel of this matrix is then a subspace of R k . The full rank assuption The lifting is, of course, not unique. It is natural to minimize the change in Hamiltonian ˙ H , assuggested by Brody, Holm, and Meier [1]. Figure 6: Interpolation between 6 qubit states visualized on the Bloch sphere. Theresulting curve is C and therefore generated by a time-dependent Hamiltonian iH ( t ) ∈ u ( n + 1) where t (cid:55)→ ˙ H ( t ) is continuous.means that the kernel has dimension k − n . This gives us a point in Gr( k − n, k ), which, bythe canonical isomorphism between Gr( k − n, k ) and Gr( n, k ), gives a point in Gr( n, k ).Since the Grassmannian Gr( n, k ) is a symmetric space, we can use our algorithm toconstruct C -splines between affine shapes.Let us study a particular case of affine shapes corresponding to 3 + 1 points in R .Our aim is to create a spline between such shapes, defined by (A)–(D) in Figure 7. Forevery shape, we obtain a 2 × , , ,
3) is isomorphic to the real projective plane R P . Indeed, this follows since the kernels of the 2 × R . We may therefore visualise the affine shapes as follows.First, we compute the intersection of the kernel direction with the half-sphere (cid:8) ( x, y, z ) ∈ R (cid:12)(cid:12) z ≤ (cid:9) . (24)We then compute the stereographic projection from the point of coordinate (0 , , x, y, z ) (cid:55)→ (cid:16) x − z , y − z (cid:17) . (25)The resulting projected spline curve is plotted in Figure 8.We stress that interpolation of affine shapes cannot, in general, be accomplished bystandard spline interpolation in vector spaces. Indeed, if we introduce local coordinates131 2 3 (a)
10 2 3 (b)
20 1 3 (c)
30 1 2 (d)
Figure 7: (A)–(D) shows different sets of labelled points p , p , p , p in R . Each suchset defines an affine shape, corresponding to an element in Gr(3 , p = ae + be in the basis spanned by the basis vectors e = p − p and e = p − p with origin at p , then the b -coordinate is infinite as the splinepasses through the affine shape (C). This is illustrated in Figure 9. More generally, anylocal 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 infuture work. Here we list some of them. • We surmise that the conditions of Theorem 2.2 are in fact valid for more generalsymmetric spaces than Riemannian ones. In particular, we would like to be ableto prove that result directly, without resorting to the previous result of Popiel andNoakes, which is the only reason for us to restrict to the
Riemannian symmetricspace case at this point. • The only structure that we need to implement Algorithm 1 and to state Theo-rem 2.2 is an invariant connection (see (26a)), which only requires the homoge-neous space at hand to be reductive [10, § • A modification of our algorithm is possible in the more general case of arbitrary
Riemannian manifolds (not necessarily symmetric), using the original C conditionof Popiel and Noakes [14]. We are working on an implementation and a proof ofconvergence. 14 .Figure 8: Using that Gr(3 , (cid:39) R P , we visualize the spline interpolating the affineshapes (A)–(B)–(C)–(D)–(A) in Figure 7 using a stereographic projection. Aswe identify any two opposite points of the circle, the curve is closed. Figure 9: This figure illustrates why it is not possible, in general, to use a local coordinatechart to construct interpolation between affine shapes. In this example, theaffine shape spline is described in a coordinate system obtained by placingthe origin at p , with basis vectors e = p − p and e = p − p , and thenexpressing p = ae + be in this basis. Then ( a, b ) constitutes local coordinatesfor Gr(3 , b –coordinate becomes infinite there. This effect is only apparent inthese local coordinates: nothing is infinite in the actual computed spline curve,as shown in Figure 8. Any other choice of local coordinates suffers analogousproblems. 15 Our algorithm is currently based on a fixed point iteration. One could instead usea Newton iteration, which would ensure convergence in one step in the Euclideancase.
Acknowledgments
This project has received funding from the European Union’s Horizon 2020 research andinnovation programme under grant agreement No 661482, from the Swedish Foundationfor Strategic Research under grant agreement ICA12-0052, and from the Knut and AliceWallenberg Foundation under grant agreement KAW-2014.0354.
A. Identities and bounds on symmetric spaces
In this section we prove several identities and bounds on (Riemannian) symmetric spaceswhich are used in Appendix B and Appendix C below.Let M = G/H p be a symmetric space, where p ∈ M and g = m p ⊕ h p is the decom-position described in § π p : g → m p be the canonical projection. Denote byΨ p : G → M and ˜ L g : M → M the mapsΨ p ( g ) = ˜ L g ( p ) = g · p, as usual, L g : G → G and R h : G → G are defined by L g ( h ) = R h ( g ) = gh and we identify g = T e G to define T e L g : g → T g G. By an abuse of notation, we will use g · to denote any of the following action maps˜ L g : M → M,T ˜ L g : T M → T M Ad g := T g R g − T e L g : g → g . The notation (cid:104)· , ·(cid:105) p is used to denote the natural pairing of one-forms and tangentvectors over a point p .The infinitesimal action (5) and the decomposition g = m p ⊕ h p define a principalconnection [10], i.e., a g -valued one-form (cid:36) on M which fulfills (cid:104) (cid:36), v (cid:105) p · p = v (consistency) (26a) (cid:104) (cid:36), g · v (cid:105) g · p = g · (cid:104) (cid:36), v (cid:105) p (equivariance) (26b)for all v ∈ T p M .By construction, m p = (cid:36) ( T p M ). 16hroughout this section and the following, dexp denotes the trivialized derivative ofthe Lie group exponential, that is T ξ exp = T e R exp( ξ ) dexp ξ = T e L exp( ξ ) dexp − ξ , where the last equality follows from differentiating exp( ξ ) exp( − ξ ) = id. Lemma A.1.
Fix p ∈ M , and define q : g → M by q ( ξ ) = exp( ξ ) · p . Then T ξ q ( δξ ) = exp( ξ ) · ( π p dexp − ξ ( δξ ) · p ) , ∀ δξ ∈ g . Proof.
By the properties of a Lie group action, for η ∈ g , we have T g Ψ p ( g · η ) = g · ( η · p )We differentiate q = Ψ p ◦ exp by the product rule and get T ξ q ( δξ ) = T exp( ξ ) Ψ p ◦ T ξ exp( δξ )= exp( ξ ) · (cid:0) dexp − ξ ( δξ ) · p (cid:1) = exp( ξ ) · ( π p dexp − ξ ( δξ ) · p )where the final equality is because π p ( η ) · p = η · p for all η ∈ g . Lemma A.2.
Fix p ∈ M . Let ξ, η ∈ m p , then π p dexp − ξ ( η ) = sinh(ad ξ )ad ξ ( η ) = π p dexp ξ ( η ) . Proof.
As a function g → g , dexp ξ is an analytic function of ad ξ with Taylor seriesdexp ξ = (cid:80) ∞ k =0 1( k +1)! ad kξ . Using the definition (4) of a symmetric space, we have, forany element η ∈ m p , that ad kξ η ∈ (cid:40) m p , k even h p , k oddThe effect of π p is thus to eliminate terms in h p , leaving only the terms with an evenexponent k , so π p dexp ξ η = ∞ (cid:88) l =0 l + 1)! ad lξ η = sinh(ad ξ )ad ξ η. Note that, as sinh(ad ξ ) / ad ξ is an even function of ξ , we have π p dexp ξ η = π p dexp − ξ η. A corollary of Lemma A.2 is that for ξ near zero, π p dexp ξ is invertible as a function m p → m p with inverse given by ad ξ / sinh ad ξ .17 emma A.3. Fix p ∈ M . Let I p be the symmetry function at p ∈ M , i.e., I p ( q ) = E p ( − log p ( q )) , and let ξ ∈ m p , δp ∈ T p M . Then T exp( ξ ) · p I p (exp( ξ ) · δp ) = − exp( − ξ ) · δp. Proof.
Let q ( ξ ) = exp( ξ ) · p . Then q : m p → M is a local diffeomorphism and I p ( q ( ξ )) = q ( − ξ ). Differentiating both q ( ξ ) and I p ( q ( ξ )) with respect to ξ , and using Lemma A.1,we get T ξ q ( η ) = exp( ξ ) · (cid:0) π p dexp − ξ ( η ) · p (cid:1) T q ( ξ ) I p i ◦ T ξ q ( η ) = − exp( − ξ ) · (cid:0) π p dexp ξ ( η ) · p (cid:1) . Since q is a diffeomorphism, we can choose η ∈ g such that T ξ q ( η ) = exp( ξ ) · w , i.e. , w = π p dexp − ξ ( η ) · p . Using Lemma A.2, we also have w = π p dexp ξ ( η ) · p , and the resultfollows. Lemma A.4.
Let g ∈ G , and p, q ∈ M . Then log g · p ( q ) = g · log p ( g − · q ) . Proof.
Let σ ( t ) = [ g · p, q ] t be the interpolating curve between g · p and q . Then g − · σ ( t ) =[ p, g − · q ] t , and differentiating at t = 0, we get g − · log g · p ( q ) = log p ( g − · q ) . 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 m p , from whichwe derive a norm (cid:107)·(cid:107) on m p . Lemma A.5.
Let ξ, η ∈ m p , then (cid:107) π p dexp ξ ( η ) (cid:107) = (cid:13)(cid:13)(cid:13)(cid:13) sinh(ad ξ )ad ξ ( η ) (cid:13)(cid:13)(cid:13)(cid:13) ≤ sinh( K (cid:107) ξ (cid:107) ) K (cid:107) ξ (cid:107) (cid:107) η (cid:107) , (cid:107) ( π p dexp ξ ) − ( η ) (cid:107) = (cid:13)(cid:13)(cid:13)(cid:13) ad ξ sinh(ad ξ ) ( η ) (cid:13)(cid:13)(cid:13)(cid:13) ≤ K (cid:107) ξ (cid:107) sin( K (cid:107) ξ (cid:107) ) (cid:107) η (cid:107) . (Notice sin , not sinh , in the last denominator.)Proof. We use Lemma A.2 to set π p dexp ξ = sinh(ad ξ )ad ξ . The lemma at hand then followsfrom considering the infinite series for sinh(ad ξ )ad ξ and its inverse, as well as the definitionof the curvature K in equation (22). Lemma A.6.
Fix p ∈ M . Let ξ , ξ ∈ m p , q = exp( ξ ) · p and q = exp( ξ ) · p . Recallfrom § d is the Riemannian distance between two points on M . Let us define U := 12 ( (cid:107) ξ (cid:107) + (cid:107) ξ (cid:107) + d ( q , q )) . e then have the inequality (cid:107) ξ − ξ (cid:107) ≤ KU sin( KU ) d ( q , q ) . Proof.
Let γ ( t ) be the geodesic from γ (0) = q to γ (1) = q , and let α ( t ) ∈ m p be suchthat exp( α ( t )) · p = γ ( t ) . (27)Then ξ = α (0), ξ = α (1). We note that, by the triangle inequality, we have (cid:107) α ( t ) (cid:107) = d ( p, γ ( t )) ≤ min { d ( p, q ) + d ( q , γ ( t )) , d ( p, q ) + d ( q , γ ( t )) }≤ min (cid:110) (cid:107) ξ (cid:107) + td ( q , q ) , (cid:107) ξ (cid:107) + (1 − t ) d ( q , q ) (cid:111) ≤ (cid:16) (cid:107) ξ (cid:107) + (cid:107) ξ (cid:107) + d ( q , q ) (cid:17) = U. (28)Differentiating (27) we getexp( α ( t )) · (dexp − α ( t ) ˙ α ( t ) · p ) = ˙ γ ( t ) , which we rewrite as ( π p dexp − α ( t ) ˙ α ( t )) · p = exp( − α ( t )) · ˙ γ ( t ) , where we have used that π p ( η ) · p = η · p for all η ∈ g . Taking norms on both sides andusing that G acts by isometries, we have (cid:107) π p dexp − α ( t ) ˙ α ( t ) (cid:107) = (cid:107) exp( − α ( t )) · ˙ γ ( t ) (cid:107) = (cid:107) ˙ γ ( t ) (cid:107) = d ( q , q ) . Moreover, by Lemma A.5, (cid:107) ˙ α ( t ) (cid:107) ≤ K (cid:107) α ( t ) (cid:107) sin( K (cid:107) α ( t ) (cid:107) ) d ( q , q ) . We now use the monotonicity of x/ sin x , and the bound (28) to obtain (cid:107) α (1) − α (0) (cid:107) ≤ (cid:90) (cid:107) ˙ α ( t ) (cid:107) d t ≤ (cid:90) K (cid:107) α ( t ) (cid:107) sin( K (cid:107) α ( t ) (cid:107) ) d ( q , q )d t ≤ KU sin( KU ) d ( q , q ) . . C 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 ∞ . We show that the spline is C at t = i ∈ { , , . . . , N − } , i.e., at the interpolating points p i . By [14, Thm. 3], asufficient and necessary condition for the spline to be C at t = i is T q + i I p i (log q + i q − i +1 ) = − log q − i ( q + i − ) − q − i ( p i ) , (29)where I p i is the symmetry function at p i . We proceed by removing the dependence of T I p i in (29) by using the lemmata from Appendix A.Recall that q + i = exp( ξ i ) · p i and q − i = exp( − ξ i ) · p i . By Lemma A.4, we havelog q + i ( q − i +1 ) = log exp( ξ i ) · p i ( q − i +1 )= exp( ξ i ) · log p i (exp( − ξ i ) · q − i +1 ) , log q − i ( q + i − ) = log exp( − ξ i ) · p i ( q + i − )= exp( − ξ i ) · log p i (exp( ξ i ) · q + i − ) . Using Lemma A.3, with w = log p i (exp( − ξ i ) · q − i +1 ), we can write T q + i I p i (log q + i q − i +1 ) = T exp( ξ i ) · p i I p i (exp( ξ i ) · log p i (exp( − ξ i ) · q − i +1 )= − exp( − ξ i ) · log p i (exp( − ξ i ) · q − i +1 ) . (30)Next, we note that log q − i ( p i ) = exp( − ξ i ) · v i . (31)Using (30) and (31), we can rewrite (29) as − exp( − ξ i ) · log p i (exp( − ξ i ) · q − i +1 ) = − exp( − ξ i ) · log p i (exp( ξ i ) · exp( ξ i − ) · p i − ) − − ξ i ) · v i . Acting on the whole equation from the left with exp( ξ i ) and rearranging, we getlog p i (exp( − ξ i ) · q − i +1 ) − v i = log p i (exp( ξ i ) · q + i − ) + v i . We now obtain the C condition (13) by using the C condition (12) to replace v i =log p i ( q + i ) = − log p i ( q − i ) in the equation above. C. Convergence of the fixed-point method
We now proceed to prove Theorem 3.1, i.e., that if consecutive interpolation points aresufficiently close then Algorithm 1 converges to a solution. In the proof, we will workwith variables in the Lie algebra g . Let m p := m p × · · · × m p N , ξ := ( ξ , . . . , ξ N ) ∈ m p , be the Lie algebra elements defined by ξ i = (cid:104) (cid:36), v i (cid:105) p i . Switching to variables in the Lie algebra, the fixed-point iteration (17) becomes a map φ : m p → m p , φ ( ξ ) = ( ξ , . . . , ξ N ) → ( φ ( ξ ) , . . . , φ N ( ξ )) . For 1 ≤ i ≤ N − φ i ( ξ ) = (cid:104) (cid:36), w i ( ξ ) (cid:105) p i + ξ i . To simplify notation, we define L p : M → m p by L p ( q ) := (cid:104) (cid:36), log p ( q ) (cid:105) p . Now, for 0 < i < N the fixed-point iteration map φ i can be expressed as φ i ( ξ ) = 12 ξ i ++ 14 (cid:16) L p i (exp( − ξ i ) exp( − ξ i +1 ) · p i +1 ) − L p i (exp(+ ξ i ) exp(+ ξ i − ) · p i − ) (cid:17) , and both φ and φ N are given by the boundary conditions φ ( ξ ) = (cid:104) (cid:36), v ( ξ ) (cid:105) , φ N ( ξ ) = (cid:104) (cid:36), v N ( ξ ) (cid:105) . Specifically, for clamped boundary conditions, we get φ ( ξ ) = ξ start , φ N ( ξ ) = ξ end (32)where ξ start and ξ end are prescribed. For free boundary conditions φ ( ξ ) = 12 (cid:104) (cid:36), log p (exp( − ξ ) · p ) (cid:105) p ,φ N ( ξ ) = − (cid:104) (cid:36), log p N (exp( ξ N − ) · p N − ) (cid:105) p N . (33)It is convenient to introduce the auxiliary variables ω + i ( ξ ) = (cid:104) (cid:36), w + i (cid:105) p i = L p i (exp( − ξ i ) exp( − ξ i +1 ) · p i +1 ) ,ω − i ( ξ ) = (cid:104) (cid:36), w − i (cid:105) p i = L p i (exp(+ ξ i ) exp(+ ξ i − ) · p i − ) , (34)and to write φ i ( ξ ) = φ + i ( ξ ) + φ − i ( ξ ) , where φ + i ( ξ ) = 14 ω + i ( ξ ) + 14 v i ,φ − i ( ξ ) = − ω − i ( ξ ) + 14 v i . (35)Furthermore, we define the norms (cid:107) ξ (cid:107) := max ≤ i ≤ N (cid:107) ξ i (cid:107) = max ≤ i ≤ N (cid:107) v i (cid:107) (36)21nd (cid:107) ω ( ξ ) (cid:107) := max ≤ i ≤ N − max( (cid:107) ω + i ( ξ ) (cid:107) , (cid:107) ω − i ( ξ ) (cid:107) ) . Recall that D denotes the maximum distance between neighbouring interpolationpoints as declared in equation (21). Notice that by the triangle inequality, we have (cid:107) ω + i ( ξ ) (cid:107) = d ( q + i , q − i +1 ) ≤ d ( p i , p i +1 ) + (cid:107) ξ i (cid:107) + (cid:107) ξ i +1 (cid:107) ≤ D + 2 (cid:107) ξ (cid:107) and similarly for (cid:107) ω − i ( ξ ) (cid:107) , so (cid:107) ω ( ξ ) (cid:107) ≤ D + 2 (cid:107) ξ (cid:107) . (37)We are now in position to prove Theorem 3.1. The proof consists in showing theexistence of an invariant region, and thereafter bounding the partial derivatives of φ toshow that it is a contraction. C.1. Invariant region
We begin by establishing an invariant region.
Proposition C.1.
Let
V > satisfy the inequality KU sin( KU ) ( U − V ) ≤ V, (38) where U = D + 2 V . If either one of the following boundary conditions is fulfilled: (i) clamped spline boundary conditions are used, (cid:107) ξ start (cid:107) ≤ V , (cid:107) ξ end (cid:107) ≤ V and (cid:107) ξ (cid:107) ≤ V , (ii) natural spline boundary conditions are used and (cid:107) ξ (cid:107) ≤ V ,then we also have (cid:107) φ ( ξ ) (cid:107) ≤ V . To reach this result, we first prove a few bounds for φ + i ( ξ ). Lemma C.2.
Let ≤ i ≤ N − , then the following bounds hold, (cid:107) φ + i ( ξ ) (cid:107) ≤ KU i sin( KU i ) ( d ( p i , p i +1 ) + (cid:107) ξ i +1 (cid:107) ) , (cid:107) φ − i ( ξ ) (cid:107) ≤ KU i − sin( KU i − ) ( d ( p i , p i − ) + (cid:107) ξ i − (cid:107) ) , where U i = d ( p i , p i +1 ) + (cid:107) ξ i +1 (cid:107) + (cid:107) ξ i (cid:107) .Proof. We prove the first inequality. The proof of the second inequality is entirelysymmetrical and therefore omitted.By definition (35), φ + i ( ξ ) = 14 ( ω + i ( ξ ) + ξ i ) . ω + i ( ξ ) + ξ i = ω + i ( ξ ) − ( − ξ i ).Recall that G acts isometrically on m p by the adjoint action (cid:107) ω + i ( ξ ) + ξ i (cid:107) = (cid:107) exp( ξ ) · ( ω + i ( ξ ) + ξ i ) (cid:107) = (cid:107) exp( ξ i ) · ω + i ( ξ ) + ξ i (cid:107) , where · ω + i ( ξ ) and ξ i are vectors in m q + i . We also haveexp( − ξ i ) · q + i = p i and exp(exp( ξ i ) · ω + i ( ξ )) · q + i = exp( ξ i ) exp( ω + i ( ξ )) exp( − ξ i ) · q + i = exp( ξ i ) exp( ω + i ( ξ )) · p i = exp( − ξ i +1 ) · p i +1 = q + i +1 . By Lemma A.6, we have (cid:13)(cid:13)(cid:13) exp( ξ i ) · ω + i ( ξ ) − ( − ξ i ) (cid:13)(cid:13)(cid:13) ≤ KT sin( KT ) d ( p i , q − i +1 )where 2 T := (cid:107) exp( ξ i ) ω i (cid:107) + (cid:107) ξ i (cid:107) + d ( p i , q − i +1 )= d ( q + i , q − i +1 ) + (cid:107) ξ i (cid:107) + d ( p i , q − i +1 ) . From the triangle inequality, we have d ( q + i , q − i +1 ) ≤ d ( p i , p i +1 ) + (cid:107) ξ i (cid:107) + (cid:107) ξ i +1 (cid:107) ,d ( p i , q − i +1 ) ≤ d ( p i , p i +1 ) + (cid:107) ξ i +1 (cid:107) . Thus T ≤ d ( p i , p i +1 ) + (cid:107) ξ i (cid:107) + (cid:107) ξ i +1 (cid:107) = U i , and the claim follows by monotonicity of thefunction x sin x . Proof of Proposition C.1.
Using Lemma C.2, and the inequalities (cid:107) ξ i (cid:107) ≤ (cid:107) ξ (cid:107) ≤ V,d ( p i , p i +1 ) ≤ D and U i ≤ U , we get that (cid:107) φ i ( ξ ) (cid:107) ≤ (cid:107) φ − i ( ξ ) (cid:107) + (cid:107) φ − i ( ξ ) (cid:107) ≤ KU sin( KU ) ( D + V )for all 1 ≤ i ≤ N − φ ( ξ ) = ξ start , so (cid:107) φ ( ξ ) (cid:107) ≤ V by our assumption.(ii) Natural spline: φ ( ξ ) = (cid:104) (cid:36), log p (exp( − ξ ) · p ) (cid:105) p , so (cid:107) φ ( ξ ) (cid:107) ≤
12 ( d ( p , p ) + (cid:107) ξ (cid:107) ) ≤
12 ( D + V ) ≤ V under a weaker condition on V . Similarly for φ N ( ξ ).23n conclusion, under the assumption (38), we have (cid:107) φ ( ξ ) (cid:107) = max ≤ i ≤ N (cid:107) φ i ( ξ ) (cid:107) ≤ V and (cid:8) ξ (cid:12)(cid:12) (cid:107) ξ (cid:107) ≤ V (cid:9) forms an invariant region. C.2. Contraction
We now establish that, for D and (cid:107) ξ (cid:107) sufficiently small, there exists an α ∈ (0 ,
1) suchthat, in the operator norm derived from (36) (cid:107) T ξ φ (cid:107) ≤ α . We first consider the effect ofvarying ξ i +1 and ξ i − in φ i ( ξ ). Proposition C.3.
Let ≤ i ≤ N − and denote by ∂φ i ( ξ ) ∂ξ j the partial derivative of φ i with respect to ξ j , i.e. ∂φ i ( ξ ) ∂ξ j is a linear operator m p j → m p i . Then, in the operatornorm, (cid:13)(cid:13)(cid:13)(cid:13) ∂φ i ( ξ ) ∂ξ i +1 (cid:13)(cid:13)(cid:13)(cid:13) ≤ (cid:13)(cid:13) ( π p i dexp − ω + i ( ξ ) ) − (cid:13)(cid:13)(cid:13)(cid:13) π p i +1 dexp ξ i +1 (cid:13)(cid:13) ≤ K (cid:107) ω + i ( ξ ) (cid:107) sin( K (cid:107) ω + i ( ξ ) (cid:107) ) sinh( K (cid:107) ξ i +1 (cid:107) ) K (cid:107) ξ i +1 (cid:107) , and (cid:13)(cid:13)(cid:13)(cid:13) ∂φ i ( ξ ) ∂ξ i − (cid:13)(cid:13)(cid:13)(cid:13) ≤ (cid:107) ( π p i dexp − ω − i ( ξ ) ) − (cid:107)(cid:107) π p i − dexp ξ i − (cid:107)≤ K (cid:107) ω − i ( ξ ) (cid:107) sin( K (cid:107) ω − i ( ξ ) (cid:107) ) sinh( K (cid:107) ξ i − (cid:107) ) K (cid:107) ξ i − (cid:107) . Proof.
We only prove the first claim. The proof of the second is entirely symmetrical.The term ξ i +1 only enters φ i ( ξ ) through the term ω + i ( ξ ), so ∂φ i ( ξ ) ∂ξ i +1 = 14 ∂ω + i ( ξ ) ∂ξ i +1 . (39)By the definition of ω + i ( ξ ), we haveexp( ω + i ( ξ )) · p i = exp( − ξ i ) exp( − ξ i +1 ) · p i +1 . Differentiating implicitly on both sides, then acting from the left by exp( − ω + i ( ξ )), weget (dexp − ω + i ( ξ ) δω + i ) · p i = − exp( − ω + i ( ξ )) exp( − ξ i ) exp( − ξ i +1 ) · (dexp ξ i +1 δξ i +1 ) · p i +1 , where δω + i = ∂ω + i ( ξ ) ∂ξ i +1 δξ i +1 . Taking norms on both sides, and using that the group actionis isometric, we get (cid:107) (dexp − ω + i ( ξ ) δω + i ) · p i (cid:107) = (cid:107) (dexp ξ i +1 δξ i +1 ) · p i +1 (cid:107) , (cid:107) π p i dexp − ω + i ( ξ ) δω + i (cid:107) = (cid:107) π p i +1 dexp ξ i +1 δξ i +1 (cid:107) . It follows that (cid:107) δω + i (cid:107) ≤ (cid:107) ( π p i dexp − ω + i ( ξ ) ) − (cid:107)(cid:107) π p i +1 dexp ξ i +1 (cid:107)(cid:107) δξ i +1 (cid:107) , thus (cid:13)(cid:13)(cid:13)(cid:13) ∂ω + i ( ξ ) ∂ξ i +1 (cid:13)(cid:13)(cid:13)(cid:13) ≤ (cid:107) ( π p i dexp − ω + i ( ξ ) ) − (cid:107)(cid:107) π p i +1 dexp ξ i +1 (cid:107) . The claims in the proposition follow from (39) and Lemma A.5.We note that we can also use a Taylor expansion of the right hand side to get theasymptotic bound (cid:13)(cid:13)(cid:13)(cid:13) ∂φ i ( ξ ) ∂ξ i +1 (cid:13)(cid:13)(cid:13)(cid:13) ≤
14 + 124 K (cid:107) ω + i ( ξ ) (cid:107) + 124 K (cid:107) ξ i (cid:107) + O ( K (cid:107) ω + i ( ξ ) (cid:107) , K (cid:107) ξ (cid:107) ) . To bound ∂φ i ( ξ ) ∂ξ i , it is again advantageous to do the splitting φ i ( ξ ) = φ + i ( ξ ) + φ − i ( ξ ),and consider each term separately.While it is possible to bound (cid:13)(cid:13)(cid:13) ∂φ + i ( ξ ) ∂ξ i (cid:13)(cid:13)(cid:13) by trigonometric and hyperbolic functions of K (cid:107) ξ i (cid:107) and K (cid:107) ω + i ( ξ ) (cid:107) , the required calculations are lengthy, and we restrict ourself toan asymptotic bound. Proposition C.4.
Let ≤ i ≤ N − and denote by ∂φ + i ( ξ ) ∂ξ i : m p i → m p i the partialderivative of φ i with respect to ξ i . In the operator norm (cid:13)(cid:13)(cid:13)(cid:13) ∂φ + i ( ξ ) ∂ξ i (cid:13)(cid:13)(cid:13)(cid:13) ≤ K (cid:18) (cid:107) ω + i ( ξ ) (cid:107) + 16 (cid:107) ξ i (cid:107) + 12 (cid:107) ξ i (cid:107)(cid:107) ω + i ( ξ ) (cid:107) (cid:19) + O ( K (cid:107) ω + i ( ξ ) (cid:107) , K (cid:107) ξ i (cid:107) ) and (cid:13)(cid:13)(cid:13)(cid:13) ∂φ − i ( ξ ) ∂ξ i (cid:13)(cid:13)(cid:13)(cid:13) ≤ K (cid:18) (cid:107) ω − i ( ξ ) (cid:107) + 16 (cid:107) ξ i (cid:107) + 12 (cid:107) ξ i (cid:107)(cid:107) ω − i ( ξ ) (cid:107) (cid:19) + O ( K (cid:107) ω − i ( ξ ) (cid:107) , K (cid:107) ξ i (cid:107) )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 ω + i . Lemma C.5.
For ≤ i ≤ N − , the following equality holds ∂ω + i ( ξ ) ∂ξ i δξ i = − ( π p i dexp − ω + i ( ξ ) ) − π p i Ad exp( − ω + i ( ξ )) dexp − ξ i δξ i . roof. We use that exp( ω + i ( ξ )) · p i = exp( − ξ i ) · q − i +1 . (40)By differentiating implicitly, we getexp( ω + i ( ξ )) · (dexp − ω + i ( ξ ) δω + i ) · p i = − (dexp − ξ i δξ i ) · exp( − ξ i ) · q − i +1 = − (dexp − ξ i δξ i ) · exp( ω + i ( ξ )) · p i , where the second equality uses (40). Acting from the left with exp( − ω + i ( ξ )), we get(dexp − ω + i ( ξ ) δω + i ( ξ )) · p i = − (Ad exp( − ω + i ( ξ )) dexp − ξ i δξ ) · p i . Applying the connection (cid:36) to both sides, we get (cid:16) π p i dexp − ω + i ( ξ ) (cid:17) δω + i = − (cid:16) π p i Ad exp( − ω + i ( ξ )) dexp − ξ i (cid:17) δξ i , where the bracketed expressions are linear operators m p i → m p i . Now, the map π p i dexp − ω + i ( ξ ) : m p i → m p i is invertible by Lemma A.2, and the claim follows.We now proceed with the proof of Proposition C.4. Proof of Proposition C.4.
From the definition of φ + i , ∂φ + i ( ξ ) ∂ξ i ( δξ i ) = 14 (cid:18) δξ i + ∂ω + i ( ξ ) ∂ξ i ( δξ i ) (cid:19) Using Lemma C.5, we have ∂φ + i ( ξ ) ∂ξ i δξ i = 14 (cid:16) δξ i − ( π p i dexp − ω + i ( ξ ) ) − π p i Ad exp( − ω + i ( ξ )) dexp − ξ i δξ i (cid:17) . (41)Using Lemma A.2, the series expansions for Ad exp( − ω + i ( ξ )) = exp( − ad ω + i ), and dexp − ξ i ,the expression on the right hand side can be expanded as an infinite series in ad ω + i ( ξ ) and ad ξ i applied to δξ i . (Recall that, as we saw in the proof of Lemma A.2, the effectof π p i is to cancel out Lie monomials of even degree.) The crucial point is that there isa cancellation in the leading term so that, neglecting terms of fourth order and higher, δξ i − ( π p i dexp − ω + i ( ξ ) ) − π p i Ad exp( − ω + i ( ξ )) dexp − ξ i δξ i = −
13 ad ω + i ( ξ ) ( δξ i ) −
16 ad ξ i ( δξ i ) −
12 ad ω + i ( ξ ) ad ξ i ( δξ i )+ O (cid:16) K (cid:107) ω + i ( ξ ) (cid:107) (cid:107) δξ i (cid:107) , K (cid:107) ξ i (cid:107) (cid:107) δξ i (cid:107) (cid:17) . (42)We thus get ∂φ + i ( ξ ) ∂ξ i = − (cid:18)
13 ad ω + i ( ξ ) + 16 ad ξ i + 12 ad ω + i ( ξ ) ad ξ i + O ( K (cid:107) ω + i ( ξ ) (cid:107) , K (cid:107) ξ i (cid:107) ) (cid:19) . The proposition follows by taking norms. 26t is also possible to bound the derivatives of the boundary condition functions.
Proposition C.6. (i) (Clamped splines) Let φ and φ N be as in (32) . Then ∂φ ( ξ ) ) ∂ξ = 0 , ∂φ N ( ξ ) ∂ξ N − = 0 . (ii) (Natural splines) Let φ and φ N be as on (33) . Then (cid:13)(cid:13)(cid:13)(cid:13) ∂φ ( ξ ) ∂ξ (cid:13)(cid:13)(cid:13)(cid:13) ≤
12 + 13 K (cid:107) φ (cid:107) + 112 K (cid:107) ξ (cid:107) + O ( K (cid:107) ξ (cid:107) ) and similar for (cid:13)(cid:13)(cid:13) ∂φ N ( ξ ) ∂ξ N − (cid:13)(cid:13)(cid:13) .Proof. (i) For clamped splines, φ and φ N are constant.(ii) For natural splines, we haveexp(2 φ ( ξ )) · p = exp( − ξ ) · p . Differentiating and taking norms, we get2 (cid:107) π p dexp − φ δφ (cid:107) = (cid:107) π p dexp ξ δξ (cid:107) , where δφ = ∂φ ( ξ ) ∂ξ δξ . Lemma A.5 now gives (cid:107) δφ (cid:107) ≤
12 2 K (cid:107) φ (cid:107) sin(2 K (cid:107) φ (cid:107) ) sinh( K (cid:107) ξ (cid:107) ) K (cid:107) ξ (cid:107) (cid:107) δξ (cid:107) , and the claim follows by a Taylor expansion of x sin( x ) and sinh( x ) x . The correspondingbound and proof for (cid:13)(cid:13)(cid:13) ∂φ N ( ξ ) ∂ξ N − (cid:13)(cid:13)(cid:13) is entirely symmetrical. C.3. Convergence
We are now ready to prove the convergence theorem.
Proof of Theorem 3.1.
A consequence of Proposition C.1 is that φ has an invariant re-gion for small enough D . We need to also establish that it is a contraction. It is sufficientto show that there exists an α ∈ (0 ,
1) such that (cid:107) T ξ φ (cid:107) ≤ α in an invariant region.We have (cid:107) T ξ φ (cid:107) = max ≤ i ≤ N (cid:107) T ξ φ i (cid:107)
27s well as (cid:107) T ξ φ (cid:107) = (cid:13)(cid:13)(cid:13)(cid:13) ∂φ ( ξ ) ∂ξ (cid:13)(cid:13)(cid:13)(cid:13) , (cid:107) T ξ φ i (cid:107) = (cid:13)(cid:13)(cid:13)(cid:13) ∂φ i ( ξ ) ∂ξ i − (cid:13)(cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13)(cid:13) ∂φ i ( ξ ) ∂ξ i (cid:13)(cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13)(cid:13) ∂φ i ( ξ ) ∂ξ i +1 (cid:13)(cid:13)(cid:13)(cid:13) , (cid:107) T ξ φ N (cid:107) = (cid:13)(cid:13)(cid:13)(cid:13) ∂φ N ( ξ ) ∂ξ N − (cid:13)(cid:13)(cid:13)(cid:13) . By Proposition C.3, Proposition C.4, and Proposition C.6, we can bound the termsappearing with a Taylor expansion to obtain (cid:107) T ξ φ i (cid:107) ≤ (cid:18) K (cid:107) ω − i ( ξ ) (cid:107) + 16 K (cid:107) ξ i (cid:107) (cid:19) + 14 K (cid:18) (cid:107) ω − i ( ξ ) (cid:107) + 16 (cid:107) ξ i (cid:107) + 12 (cid:107) ξ i (cid:107)(cid:107) ω − i ( ξ ) (cid:107) (cid:19) + 14 K (cid:18) (cid:107) ω + i ( ξ ) (cid:107) + 16 (cid:107) ξ i (cid:107) + 12 (cid:107) ξ i (cid:107)(cid:107) ω + i ( ξ ) (cid:107) (cid:19) + 14 (cid:18) K (cid:107) ω + i ( ξ ) (cid:107) + 16 K (cid:107) ξ i (cid:107) (cid:19) + O ( K (cid:107) ω (cid:107) , K (cid:107) ξ (cid:107) ) ≤
12 + 14 K (cid:107) ω (cid:107) + 16 K (cid:107) ξ (cid:107) + 14 K (cid:107) ξ (cid:107)(cid:107) ω (cid:107) + O ( K (cid:107) ω (cid:107) , K (cid:107) ξ (cid:107) ) . We can now use (37) and Proposition C.6 to obtain a bound of the form (cid:107) T ξ φ (cid:107) ≤
12 + O ( K D , K (cid:107) ξ (cid:107) ) . It is therefore clear that when K max( D, (cid:107) ξ (cid:107) ) is sufficiently small, φ 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 • (cid:110) ξ (cid:12)(cid:12)(cid:12) (cid:107) ξ (cid:107) ≤ V (cid:111) is a compact invariant region and • φ 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. D. Python implementation
Here we briefly indicate how to use our
Python implementation to compute splines invarious geometries. First, one has to install the following package by following theinstallation instructions: https://github.com/olivierverdier/bsplinelab .0 0.1 0.2 0.30.00.20.40.6 condition (38) (cid:107) T ξ φ (cid:107) < (cid:107) ξ (cid:107) < KVKV KD (cid:107) ξ (cid:107) Figure 10: An illustration of the convergence argument. The upper curved region (lightblue) denotes the values (
KD, KV ) satisfying (38). The lower rectangular re-gion (light orange) illustrates the condition K max( D, (cid:107) ξ (cid:107) ) sufficiently small.“Sufficiently small” is not quantified, and the region is just an illustration.To ensure a solution, KV has to be in the upper curved region, and the linebelow KV has to be contained in the lower rectangular region.The minimal necessary import is from bspline . interpolation import Symmetric , cubic_spline
Here is an example of how to interpolate between random points in the Euclideancase:
The result is then a C function b which is equal to the interpolation points at integerpoints on the interval [0 , from bspline . geometry import Sphere
The function b is now a function defined on the interval [0 ,
2] such that it is C , isinterpolating the prescribed points at the integer points 0, 1 and 2, and is on the spherefor all points in between.For further information on bsplinelab , we refer to the package documentation andthe available code examples. 29 eferences [1] D. C. Brody, D. D. Holm, and D. M. Meier, Quantum splines, Phys. Rev. Lett. (2012), 100501.[2] P. Crouch and F. S. Leite, Geometry and the dynamic interpolation problem,
Amer-ican Control Conference, 1991 , pp. 1131–1136, IEEE, 1991.[3] D. Dong and I. R. Petersen, Quantum control theory and applications: a survey,
IET Control Theory & Applications (2010), 2651–2671.[4] S. Helgason, Differential Geometry and Symmetric Spaces , American MathematicalSociety, 1962.[5] D. Kendall, D. Barden, T. Carne, and H. Le,
Shape and Shape Theory , Wiley, 2009.[6] T. W. B. Kibble, Relativistic models of nonlinear quantum mechanics,
Comm. Math.Phys. (1978), 73–82.[7] T. W. B. Kibble, Geometrization of quantum mechanics, Comm. Math. Phys. (1979), 189–201.[8] K. A. Krakowski, L. Machado, F. Silva Leite, and J. Batista, A modified Casteljaualgorithm to solve interpolation problems on Stiefel manifolds, J. Comput. Appl.Math. (2017), 84–99.[9] J. M. Lee,
Riemannian manifolds , Springer-Verlag, New York, 1997.[10] H. Z. Munthe-Kaas and O. Verdier, Integrators on homogeneous spaces: Isotropychoice and connections,
Found. Comput. Math. (2016), 899–939.[11] L. Noakes, Duality and Riemannian cubics, Adv. Comput. Math. (2006), 195–209.[12] L. Noakes, G. Heinzinger, and B. Paden, Cubic splines on curved spaces, IMA J.Math. Control Inform. (1989), 465–473.[13] F. C. Park and B. Ravani, B´ezier Curves on Riemannian Manifolds and Lie Groupswith Kinematics Applications, J. Mech. Design (1995), 36.[14] T. Popiel and L. Noakes, B´ezier curves and C interpolation in Riemannian mani-folds, J. Approx. Theory (2007), 111–127.[15] C. Samir, P.-A. Absil, A. Srivastava, and E. Klassen, A Gradient-Descent Methodfor Curve Fitting on Riemannian Manifolds,