About simple variational splines from the Hamiltonian viewpoint

In this paper, we study simple splines on a Riemannian manifold $Q$ from the point of view of the Pontryagin maximum principle (PMP) in optimal control theory. The control problem consists in finding smooth curves matching two given tangent vectors with the control being the curve's acceleration, while minimizing a given cost functional. We focus on cubic splines (quadratic cost function) and on time-minimal splines (constant cost function) under bounded acceleration. We present a general strategy to solve for the optimal hamiltonian within the PMP framework based on splitting the variables by means of a linear connection. We write down the corresponding hamiltonian equations in intrinsic form and study the corresponding hamiltonian dynamics in the case $Q$ is the $2$-sphere. We also elaborate on possible applications, including landmark cometrics in computational anatomy.


Introduction
Cubic Riemannian splines and their higher order extensions have been instrumental for interpolation and statistics on manifolds. The bibliography is vast, see e.g. [85], [50], [49], [26], [46], [47], [77], [95], [43]. In this note we consider only the simplest case, namely, of splines having the tangent bundle as state space, the acceleration vector being the control. We remark that a tangent vector gives a minimal model for a short process. For instance, the rendezvous problem in robotics and in space science ( [45], [91], [70]) consists of planning a path with prescribed initial and end tangent vectors. For instance, to achieve a smooth docking of a service spacecraft to the International Space Station.
There is a potential use of simple splines also in sport science, computer animation, recognition, and video/movies repair, see [15] and references therein.
Let us now present our general framework, which belongs to the class of so-called mechanical control problems. These were first studied via Geometric Mechanics in Andrew Lewis thesis [66], consolidating earlier work by several authors, specially R. Brocket, J. Baillieul, A. Van der Schaft, P. Crouch and A. Bloch. The standard reference is the book by Bullo and Lewis [24], where the Lagrangian viewpoint is mostly used.
One considers a configuration space Q endowed with a Riemannian (kinetic energy) metric g. The organism or device under study is controlled by forces that produce an acceleration u, taken here as the control variable. If ∇ denotes the Levi-Civita connection, the state equation on the tangent bundle T Q is The aim is to connect two tangent vectors (or 'short processes') vq 0 , wq 1 ∈ T Q minimizing a cost functional. In this paper, we focus on the following two special cases which have received special attention.
Cubic splines on a Riemannian manifold, which were introduced around 1990 ( [81,34]), correspond to minimizing the cost functional: with prescribed time T . Cubic splines have been extensively used in computational anatomy.
The other special case is the time minimal problem under bounded acceleration. It consists in connecting two vectors in minimum time, under the restriction |u| ≤ A , where A is a prescribed bound. ( There is no accessibility issue: for any arbitrarily small bound on |u|, under mild hypothesis on Q it is possible to concatenate any two tangent vectors by a smooth curve with non vanishing velocity and acceleration norm ≤ [99]. The time minimal-bounded acceleration problem is in general (though not always) equivalent to the L ∞ control problem considered recently by Noakes and Kaya [58], [80], where they ask for a trajectory that minimizes the sup of the norms of the accelerations, with fixed transition time.
From now on we shall refer to cubic splines also as L 2 -splines, due to the form of the underlying cost functional (2), and to time minimal-bounded acceleration splines as L ∞ -splines, due to the previous discussion.
Via Pontryagin's maximum principle, every cost functional associated to (1) yields a Hamiltonian system in T * (T Q). In this paper, we study the resulting Hamiltonian systems in detail.

Motivations
Splines in S m .
According to Lyle Noakes, "the problem of interpolating and approximating spherical data in the m-dimensional manifold S m is much more widespread than might at first be thought" [78]. Indeed, see [98], [41], [18], [37] for information on spherical statistics.
In order to match two infinitesimal processes on T S m with simple variational splines, it seems to us that it is sufficient to understand the cases m = 1, 2, 3. This is due to homogeneity under SO(m + 1) action: any two tangent vectors vq 1 , wq 2 on 'big' T S N are actually tangent to a isometrically imbedded S m with m ≤ 3.
(For the analogous problem on R N , q1 can be the origin, v1, q1 − q2 and v2 define at most a three dimensional subspace. ) We have collected a number of references on cubic splines on S 2 . Somewhat surprisingly, as far as we know, the reduction of the SO(3)-symmetry of T * (T S 2 ) by Hamiltonian methods was still awaiting. We present here a reduction procedure, but we must confess that it works well only outside the zero section.

Splines in Diff
The systems studied in this paper can be taken as finite-dimensional toy models for the following infinite-dimensional one. Following the notations of [51], let D a domain in R d and Diff(D), X respectively the group of diffeomorphisms φ and vector fields u(x) with appropriate boundary conditions. The idea is to consider control systems as in (1) but now on the infinite dimensional Q = Diff. Upon right translations, one has T Diff = Diff × X , T * Diff = Diff × X * . With due care to functional analysis, X * is the space of momentum densities m, with m = m · dx ⊗ dV. Consider a Sobolev metric on Diff ( [22], [48], [72], [12], [13], [16], [17], [11]). One can write where G(x, y) is the Green function for the inverse of the linear partial differential operator that yields m when applied to u. EPDiff (geodesics on Diff), the celebrated Euler-Poincaré partial differential equation in terms of the momentum density is given by and comes from a noble tradition going back to Arnold's interpretation of Euler's incompressible fluid equations as geodesics in the infinite dimensional Lie group of volume preserving diffeomorphisms [75]. One of the striking facts is that often the solutions of EPDiff tend to concentrate on pulson submanifolds In turn, from singular momentum solutions one recovers the vector fields via Mario Michelli [71,72] implemented geodesic equations for landmarks. His formulae for for Christoffel symbols and curvatures coefficients in terms of cometrics, can be used to implement a code for landmark splines using our methodology. In the infinite dimensional case cubic splines have only recently been considered [89], [90]. We will discuss the open problem of relating landmark splines to the infinite dimensional problem in the final section.

Contents of the paper
In section 2, we study the control problem associated to the state equation (1) from the point of view of Pontryagin maximum principle (PMP). The main theoretical results of the paper are intrinsic formulas for the optimal hamiltonian and the resulting hamiltonian equations (Propositions 1 and 2). The key idea is to use split variables for T * (T Q) coming from taking horizontal and vertical components w.r.t. the underlying linear connection on T Q. A dual splitting of T (T Q) was already explored by Lewis and Murray and is presented in full details in Bullo and Lewis [25]. In [7] we generalize the construction to T * A, for A → Q a vector bundle with a connection. In split variables, the symplectic form is no longer canonical (it contains curvature terms, Proposition 1) but the parametric hamiltonian appearing from the PMP has a simple form. We can then find the optimal hamiltonian for cubic and time-minimal splines easily and derive the corresponding hamiltonian equations in intrinsic form (Proposition 2). We show that, in the particular case Q = S n a n-sphere, we recover the higher dimensional case analogue to the hamiltonian system for cubic splines derived by Crouch and Leite [34,33] for S 2 .
In section 3, we study dynamical aspects of the resulting hamiltonian equations on T * (T Q) in the particular case of cubic splines on Q = S 2 . Our main technical tool in this study is reduction by the natural rotation symmetry. The reduced system has two degrees of freedom. Two known families of solutions for cubic splines on S 2 are reinterpreted.
They correspond to equilibria and partial equilibria of the reduced system. i) "Figure eights" formed by two kissing circles with geodesic curvature κg = 1, run uniformly in time. They correspond to unstable, loxodromic, fixed points of the reduced Hamiltonian. ii) Equators, run cubically on time, correspond to "partially fixed" points of the reduced system. Moreover, in section 3.3 we discuss some numerical simulations. Poincaré sections indicate that the reduced system is non-integrable, but has regions rich of invariant tori. In section 3.4, we discuss the limitations of the reduction procedure (coming from the rotation action not being free at zero tangent vectors) and provide the explicit relation to the general split variables approach.
Section 4 contains some further comments and presents suggestions for further research. In particular, since it is well known that landmark geodesics lift to geodesic solutions on the full Diff [71], [23], we pose the question if landmark splines can be lifted to splines on Diff.
In appendix A, we present a simple Fortran program used for the reconstruction of trajectories in S 2 . It uses the Runge-Kutta-Fehlberg routine RK78 ( [39]) which is of standard use in Celestial Mechanics [5]. In appendix B, we discuss how to interpret the controls entering the state equation (1) in terms of curvatures of the underlying curve, in the particular case Q = Σ ⊂ R 3 is a convex surface. These results are applied in section 3.

Optimal control on T Q as state space
As mentioned in the introduction, our main object of study will be the following optimal control problem. The state space is T Q, where (Q, g) is a Riemannian manifold, and the control is represented by a (force or acceleration) vector field u ∈ X (Q). The state equations (1) can be written in first-order form as v =ẋ, ∇ẋv = u, (x(t), v(t)) ∈ T Q.
Above, ∇ẋ represents covariant derivative with respect to the Levi-Civita connection ∇ on T Q so that, for u = 0, we get the equations for the geodesic flow on T Q. We impose the boundary conditions that (x(t), v(t)) is fixed to be given (x0, v0), (x f , v f ) ∈ T Q at initial and final times.
To recover the second-order equation (1), we notice that the first equation in (7) above says that we are dealing with a second order problem for x(t) ∈ Q. More precisely, the above equations correspond to the second order vector field U ∈ X (2nd) (T Q) given by where hor∇ denotes the horizontal lift w.r.t. ∇ and u vert the natural vertical lift TxQ → In standard local coordinates for T Q, the state equations reduce toẋ ij are the Cristoffel symbols of ∇ and sum over repeated indices is understood. Following the introduction further, we shall also consider an optimization component in the problem: γ(t) = (x(t), v(t)) must also minimize a cost functional of the form with C : T (T Q) → R a given cost function. The control u is also (possibly) subjected to a constraint of the form Applying Pontryagin's maximum principle.
Our general strategy to attack the above optimization problem will be to apply PMP and transform it into a hamiltonian system on T * (T Q). We briefly describe here the principle in our present situation, however we assume that the reader is acquainted with the general recipe for PMP (otherwise we suggest [87] for a tutorial book, [84] for the fundamental reference in the area).
The general key idea is, for each state variable, to introduce a new co-state variable. This leads one to consider T * (T Q) endowed with its canonical symplectic form ω0 and consider the u-family of hamiltonians Hu ∈ C ∞ (T * (T Q)) given by where ·, · denotes the natural pairing between covectors and vectors. The PMP then states that the solution to our optimal control problem is a trajectory (for some suitable initial conditions to be found) of the hamiltonian system (T * (T Q), ω0, H * ) where Hu is the optimal hamiltonian function.
In local coordinates, denoting y k and z k the conjugated coordinates to x k and v k respectively, we have It is easy to optimize this local hamiltonian in the u k 's. In general, though, it is not so easy to find the optimal hamiltonian in an intrinsic way (i.e. globaly w.r.t. the manifold Q).
In the following subsections we propose a method to achieve this and to write down the corresponding hamilton equations in intrinsic form. This method is general and makes use of 'global splittings' of variables in which the Hu becomes simple (and hence easier to optimize) but the symplectic structure is no longer in canonical form (it incorporates curvature terms).
Other methods simplifying the optimization of Hu can be available in particular cases, in which case we provide the dictionary between the two.

Hamiltonian equations from PMP in split variables
In this subsection, we introduce global split variables to solve the optimal hamiltonian described above. These results can be seen as the Hamiltonian analogue of some results given in Lagrangian form in the supplementary materials of Bullo and Lewis' book [25].
In section 4 we outline further developments: a splitting for T * A, with A being an affine bundle with connection [7].
The main observation is that the linear connection ∇ on q : T Q → Q allows one to decompose tangent vectors into horizontal and vertical components: In more differential geometric terms, ∇ induces an Ehresmann connection for the submersion q : T Q → Q given by the bundle projection. The vertical component of X ∈ T (x,v) (T Q) can be written as Xv = Θ∇(X) vert with Θ∇ ∈ Ω 1 (T Q, T Q) a vector valued 1-form encoding the vertical projection.
Since both horizontal and vertical spaces can be identified with TxQ by taking horizontal and vertical lifts, the above decomposition defines a global diffeomorphism 1 The dual decomposition similarly induces a splitting diffeomorphism The covectors p and α above define split variables for every co-vector P ∈ T * (x,v) (T Q). For clarity, let us examine these diffeomorphisms in local coordinates. Let (pi,αj,ṽ k ,x k ) be canonical coordinates in T * T Q relative to standard ones (ṽ k ,x k ) on T Q and let (pi, αj, v k , x k ) be natural coordinates on q * T * Q ⊕T Q q * T * Q. Then, The next proposition shows the effect of using global split variables in the symplectic form and the underlying general form of Hamilton's equations.
, the Hamiltonian vector field XH is given in local coordinates bẏ Proof. Both the l.h.s. and the r.h.s. of equations (11) and (12) define global differential forms on q * T * Q ⊕T Q q * T * Q. To prove (i) and (ii) it is then enough to show that, when restricted to any coordinate chart, the corresponding local expressions of the l.h.s. and of the r.h.s. coincide. Now, the formulas after the ≡ symbols evidently correspond to the local coordinate expressions for the r.h.s'. Hence, we only need to show that the local expressions of ψ * ∇ θT Q and ψ * ∇ ΩT Q coincide with the given ones. Let us then choose coordinates for T * T Q and q * T * Q ⊕T Q q * T * Q as in eq. (10). We have that θT Q ≡pidx i +αidṽ i and thus (i) follows directly by computing the pullback ψ * ∇ θT Q following the change of coordinates (10). For (ii), we observe that , so that we need to apply −d to the known local expression for θ∇. By direct computation using the following local expression for the curvature tensor We recognize that the l.h.s'. of the last two equations in (iii) correspond to covariant derivatives. The equations for theṗi are more intricate, but we will show below that they simplify for spline control problems with cost functions of a special form.

Optimal spline Hamiltonians
Let us now show how the split variables simplify the hamiltonian of our problem. Indeed, recalling Hu defined in (8), then is a (typically convex) function c : R+ → R of the norm square of the control u. Cubic splines have cost functional (2) and are thus a particular case of the above with c being a linear function. Time-minimal splines are also a particular case with c ≡ −1 a constant (recall that, in this case, u is constraint by (3)).
It is now easy to find the optimal value H * ,∇ of Hu,∇: where g −1 is the cometric, the optimal value of the control is and Leg(c) denotes the Legendre-Fenchel dual of c [86]: We now show that a dramatic simplification in theṗi equation results from the connection preserving the metric. Proposition 2. Hamilton's equations for H * ,∇ in the case of cost functions of the form C = c(g(u, u)) can be intrinsically written aṡ Locally, they reaḋ Proof. From item (iii) in Proposition 1 and the definition of H * ≡ H * ,∇ given in (14), we immediately getẋ The equation forv reduces to Deriving with respect to αi both sides of the eq. (15) defining Leg(c) one gets We are thus only left with the equation (13) forṗ. Transporting the term The proof will be finished when we show that the term between brackets in the r.h.s. vanishes. This, in turn, follows by virtue of the fact that ∇ preserves the metric g, and hence∇ preserves g −1 ≡ (g ab (x)), so that Finally, this identity directly implies the desired vanishing:

Cubic and time minimal splines
Let us examine the particular cases mentioned in the introduction. For cubic splines, the cost functional is (2) so The advantage of the split variables is now apparent: the optimal hamiltonian is immediately given by where, for any covector α, one defines α by g(α , v) = α(v) for all vectors v. Likewise, the optimal Hamiltonian for the time minimal problem with constraint (3) is also easily shown to be The equations of motion in both cases are (16) with the corresponding u * from (17) and (18). For instance, taking β = 1 in the cubic splines problem, then ∇ẋẋ = u * = α .
Deriving this equation covariantly two more times and using the equations of motion for α and p we get 4 ∇ recovering the equations found by L. Noakes, G. Heinzinger and B. Paden [81], and P. Crouch and F. S. Leite [34].
Before moving on to splines on spheres, we present some comments about landmark cometrics. For N = 1 one has the euclidian metric on R d , for which L 2 splines are cubic polynomials on each coordinate. The L ∞ problem has been addressed in [27] for any value of d (actually d = 2, 3 is enough). The next simplest nontrivial case is d = 1, N = 2. One observes that the underlying geodesic problem (i.e. when the control u = 0) is integrable, with Hamiltonian In the spline problem (u = 0) one also has invariance under translations on the line, so there will be a conserved momentum and it will be reducible to 3 degrees of freedom. Numerical experiments for landmark splines are in order. In "Mario's formulas", partial derivatives are computed on cometric entries [71]. At every computation step, done at the current landmark locations, there is only one matrix inversion, of the cometric matrix, which has a block structure. For simulations we suggest using the Cauchy kernel G(x1, x2) = 1/(1+|x1 −x2| 2 ) that has a weaker decay at infinity than kernels involving the exponential.

Cubic splines on Q = S n and extrinsic vs intrinsic description
In this case, there is an alternative approach to finding the optimal hamiltonian which uses extrinsic variables coming from the embedding S n ⊂ R n+1 . We shall show below how to explicitly relate the two descriptions.
We first follow Dong-Eui Chang [28], fix a sphere S n (r) of radius r and consider the following state equations in R 2(n+1) To avoid confusions with scalar velocities used later, we use boldface for velocity vectors from now on. The idea is that these equations have T S n (r) ⊂ R 2(n+1) as invariant submanifolds and induce the correct state equations (7) on the sphere. The natural coordinates on R 2(n+1) restrict to T S n (r) yielding (local, but almost global) coordinates that we call extrinsic variables and denote by (x, v). Following the notation of section 2.1, we consider the cotangent bundle T * (T R n+1 ) = R 4(n+1) with coordinates (x, v,p,α) and canonical symplectic form In the case of cubic splines, the parametric Hamiltonian (8) in the ambient R 4(n+1) iŝ where the controls are restricted to the tangent planes: u ⊥ x. Notice that the ambient scalar product · allows us to identify vectors and covectors. Let us consider the projections onto the plane perpendicular to x, so that (x, v,p , α ) define extrinsic variables for T * (T S n (r)) ⊂ R 4(n+1) . It is immediate to deduce that, upon restriction Hu =Ĥu| T * (T S n (r)) , the optimal control is u * = 1 β α and that the optimal Hamiltonian reads The relation between the extrinsic variables (x, v,p , α ) and the split variables (x, v, p, α) of section 2.1 is given by the following: The split variables (9) for T * (T S n ) are given by Of course, expressing the Hamiltonian (22) in terms of the split variables we get which is the general optimal hamiltonian (17). We stress the evident simplification operated on the Hamiltonian (22) when passing to split variables.

Hamiltonian equations for cubic splines in S n and Crouch-Leite equations
To write down the equations of motion (16) in this particular case, let us first notice that for a curve (x(t), w(t)) ∈ T S n (r) described in extrinsic variables, we have ∇ẋw =ẇ + (w ·ẋ)x.
Indeed, the above implies (∇ẋw) · x = 0 by virtue of w · x = 0, so that the covariant derivative remains tangent to the sphere. Secondly, the curvature tensor of the sphere can be expressed in terms of the ambient inner product as Introducing the following change of notation and recalling u * = α (≡ α) in the cubic spline case, it immediately follows that eqs. (16) Notice that, by the previous general results, the above system automatically implies the well known nonlinear equation (19) for cubic splines.
Remark 1. The above system of equations reproduces the system derived by Crouch and Leite [34,33] for cubic splines in the case of the (n = 2)-sphere. Notice that our results imply, in particular, that this system is Hamiltonian for any n. We shall come back to these equations for S 2 in the next section.
In this section, we want to analyze the solutions of the hamiltonian equations (16) coming from the PMP in the particular case of Q = S 2 . In this case, the SO(3) symmetry plays an important role: we can perform symplectic reduction to decrease the dimension of the system. In order to simplify the reduction-reconstruction procedure, we will use a description of non-zero tangent vectors which is based on the Gauss map for a convex hypersurface 5 and detailed in Appendix B. In section 3.4 the dictionary between the above variables and the general split ones of Proposition 2 is described. We also discuss the artificial singularity at v = 0 introduced by the above identification.

Hamiltonian equations on T * (T S 2 − 0)
Let us consider a 2-sphere of radius r and take M := T S 2 − 0 to be the manifold given by all non-zero tangent vectors to the sphere. Following appendix B, the map where {e1, e2, e3} ∈ R 3 denote the standard basis vectors, is a diffeomorphism (c.f. (50)).
(One should not confuse the scalar v ∈ R+ with notation previously used for vectors.) Since the underlying surface Σ = S 2 is a sphere, this diffeomorphism preserves the natural SO(3)-actions (on M is given by the tangent lift of the action by rotations on S 2 ). Following equation (54) in appendix B, the state equations (1) for curves in T S 2 with non-zero velocity (i.e. lying in M ) readv and (u1, u2) ∈ R 2 being the controls. Note that u1 represents the tangential acceleration and u2 = v 2 κg, where κg is the geodesic curvature. The skew-symmetric matrix X ∈ so(3) can be conveniently represented as We recall that, for cubic splines, the cost function is while for the time-minimal problem: C = 1 and u = (u 1 , u 2 ) is constraint by

Applying Pontryagin's principle
The cotangent bundle of the state space M is where we used the standard left trivialization of the cotangent bundle of SO(3). We denote by a ∈ R the conjugate to v ∈ R + and (R, M 1 , M 2 , M 3 ) ∈ SO(3) × R 3 the left trivialized covectors on the rotation group. The parametric hamiltonian (8) yields in this case Now, finding the optimal Hamiltonian is a trivial task. For cubic splines so that For time minimal splines,

Hamiltonian equations for reduction-reconstruction
The symplectic structure on T * M T * R + × (SO(3) × R 3 ) is the product of the canonical one on the first factor and the very well known one on the second factor (e.g. from the rigid body problem). It is then straightforward to derive the hamiltonian equations coming from H * . Moreover, we observe that H * does not depend on R so that it descends to a reduced hamiltonian function on The induced Poisson brackets on M red are also well known The equations on T * M can be thus split into the reconstruction equations for R(t) ∈ SO(3),Ṙ 13 The function is a Casimir and restricts the dynamics of M to a momentum sphere.
For the rest of this section, we shall concentrate on the case of cubic splines. In this case, H * is given by (30) and, then, the reduced equations The study of time minimal case on S 2 will be submitted elsewhere [61].

Dynamics of cubic splines from reduced system's fixed points
Equators: linearly accelerating geodesics The simplest fixed points for eqs. (35) in the virtual momentum sphere correspond to the values M 1 = M 3 = 0, M 2 = µ (we allow µ positive or negative). Then, so to speak, the 'poles' on the virtual momentum sphere are in the second coordinate M 2 . The variables a and v follow, respectively, a linear and a quadratic function of time: Since Ω = (0, v(t)/r, 0) it is easy to reconstruct R via (32). In fact, setting R(t o ) = I, in the x 1 − x 3 plane there is a family of trajectories passsing at t = t o trough the north pole of the physical sphere (radius r): which honors the name "cubic" splines.

More fixed points of the reduced system
Fixed points can be parametrized by v ∈ R + since stationary points of (35) must satisfy a = 0, M 2 = r M 2 3 /(β v 3 ) and M parallel to ( 0, v/r, M 3 /(βv 2 ) ). A simple algebraic manipulation yields Proposition 4. For each v > 0 there are two equilibria with µ = √ 2 βv 3 /r and These fixed points correspond to relative equilibria in the unreduced system. From (29) we have u * 2 = M 3 /(βv) while from the general state equation on T S 2 one deduces u * 2 = κ g v 2 where κ g denotes the geodesic curvature of the underlying curve γ(t) ∈ S 2 (c.f. appendix B). Since M 3 = ±β v 3 r , we get Recall that on a sphere of radius r, the parallel of latitude θ has geodesic curvature κ g = tan θ/r and thus θ = π/4. Moreover, we observe that Proposition 5. (Figure eights.) The reconstructed curves in S 2 with R(0) = I, corresponding to the two equilibria parametrized by v > 0 as above, are two orthogonal (touching) circles making a 45 • angle with the equatorial plane. They are given by Proof. Since u * 2 = M 3 βv , M 3 = ± βv 3 r it follows that u * 2 = ± v 2 r anḋ

So we have steady rotations with angular velocity
Recall that for an unit vector (u x , u y , u z ) the rotation matrix R(θ) with R(0) = I is given by In hindsight, we could allow v < 0 in (37), so we can describe both twin circles in both directions. We have therefore four solutions, each twin pair starting at the north pole (0, 0, r) with velocity vector (v, 0, 0).

Discrete symmetries
They are in correspondence with expected geometric symmetries in the full system (in S 2 ). i) Reflecting a solution curve γ(t) over the equator that it is tangent at a given point. Hence, given a solution, one gets infinitely many others (but two successive reflections correspond to the action of an SO(3) element on the original curve). ii) Velocity reversal 6 iii) Time reversal: it implies that there is a symmetry between stable and unstable manifolds that perhaps could be numerically explored. Proposition 6. Discrete symmetries.
(ii) velocity reversal: These symmetries may be useful in finding periodic orbits via global calculus of variations, and perhaps find their stability using symplectic techniques [29], [68].

The fixed points are focus-focus singularities
In order to linearize about the equilibria it is convenient to take spherical coordinates on the momentum sphere, The reduced system is confined to the symplectic manifold M µ := T * R + ×S 2 µ , where S 2 µ is the momentum sphere of radius |µ| (and recall that T * R + = {(v, a) : v > 0}).
We can also define z = sin φ so that the symplectic form on M µ becomes and the reduced optimal Hamiltonian writes as (recall z = sin φ) The equilibria are We add v to the parameters µ, r, β. It turns out that the matrix that linearizes the Hamiltonian system given by (41) and (42) does not depend on µ and is the same for both equilibria: Furthermore, its characteristic polynomial does not depend on β.
Proposition 7. The eigenvalues at the fixed points (36) (equivalently (43)) are loxodromic (focus-focus type) In T * T S 2 , the union for all v = 0 of these special circle solutions with κ g = 1/r forms a center manifold C of dimension 4.
In the reduced space we have local unstable and stable (spiralling) manifolds of dimension two. They lift to 6-dimensional stable and unstable manifolds W s C , W u C inside T * T S 2 . This dimension count is coherent with dimC = 6 + 6 − 8 = 4.
Several global dynamical question can now be posed: on the unreduced system, take initial conditions near the focus-focus equilibrium. What happens with their solutions and with the corresponding unreduced solutions?
More precisely, understanding the global behavior of W u and W s is in order. Do they intersect transversally?
Are equators in the 'periphery' of phase space?
The equations of motion corresponding to the symplectic form (41) and the Hamiltonian (42) are given bẏ Assuming that v = 0 is a regularizable singularity (taking into account the various symmetries of Proposition 6, translated to these coordinates), we have a, v, θ ∈ R , |z| ≤ 1 .
The horizontal lines z = ±1 are invariant, equivalent to M 1 = M 3 = 0, M 2 = ±µ. We know from the previous discussion that reconstruction yields the equators in the unreduced system. The coordinate a runs uniformly in time from left to right at z = −1 and from right to left at z = +1, namely a(t) = −sign(z)µt/r + a o . As we expect, v is quadratic on time, with leading term −sign(z)µt 2 /(2rβ).
As for θ, for |t| sufficiently large the second term in the equation forθ can be dropped out. Thus for such large |t| we have θ(t) ∼ −sign(z)µt 3 /(6rβ).
This means that except possibly at intermediate times, the horizontal invariant θ lines in the plane (θ, z) run in opposite ways 7 for z = ±1.

Simulations of S 2 cubic splines
Numerical work on sphere splines include (we apologize for omissions), [52], [78] (a survey for the computational geometry community), [88] (a gradient descent method). Here we present some experiments using the reduction to two degrees of freedom. Besides the 'figure eights' of Prop.5, the other family of solutions that seems to have fundamental dynamical importance are the equators, described in 3.2.
Understanding the dynamics near the equators is important both conceptually and numerically (see figs 6.1 and 6.2 in [79]). Linearization does not help. It is easy to see from the general cubic spline equation (19) that for any Riemannian metric geodesics whose accelerations vary linearly are L 2 splines. Our original (uninformed) guess was that, for any Riemannian manifold, splines would tend to accelerating geodesics as t → ±∞.
Our numerical experiments indicate that this is not the case for cubic splines on S 2 . They strongly suggest that the system is non-integrable. However, we found zones in the reduced phase space having invariant tori.

Invariant tori
Figs.6 and 7 depict invariant tori. They were found, among several others with complicated Lagrangian projections, by (obsessive) trial and error experimentation. They live on remote regions of phase space, neither close to the loxodromic equilibria nor to the equators z = ±1. There is a variety of confined trajectories whose 'morphology' merit further study.

Trajectories emanating from the focus-focus
The simulations suggest that (at least some of them) are reaching a neighborhood of an equator. Will they eventually recur back to the focus-focus loxodromic equilibrium (figure eights of the unreduced)? Our simulations suggest that as t → ∞ those trajectories seem to 'orbit' around the equator z = −1, but they stay at a 'safe' distance to it (see Figure 8).
We assert that this is not an artifact of the integrator. Here's an heuristic argument. From (48) it follows that when v → ∞ thenθ ∼ v/r, and θ also diverges (only more so). Thereforeż = O(v −2 ) → 0 is an ever oscillatory way. This suggests that z perhaps could stabilize close to z = −1, but at a distance of it.
A study of (48) is in order. The behavior of the systems for moderate values of v seems quite unpredictable.

Zero velocities and the relation to the general approach via split variables
The price we had to pay in the description given so far of cubic splines in S 2 is that the scalar velocity v appears in denominators of (35), and it may vanish along a solution in finite time. Nonetheless, we expect that this is regularizable. In other words, we posit that troubles at v = 0 are (unfortunate) artifacts of our parameterization of T S 2 which excludes zero velocities and, ultimately, of the reduction procedure we implemented. (Notice that the lifted SO(3) action on T S 2 is not principal; it is so when restricted to M = T S 2 − 0.) We must then go back to the beautiful equations (23) for the unreduced system (taking n = 2 there). As mentioned earlier, these were originally derived by Crouch and Leite [34,33] (we also obtained them from our general split variables recipe in section 2.2). Clearly, unreduced solutions pass through v = x 1 = 0 as smoothly as anywhere else.
Regularization consists in lifting a reduced solution such that v(t o ) = 0 to an unreduced solution. In Crouch-Leite system there is no trouble passing though x 1 (t o ) = 0, and projecting back the continued unreduced solution 8 . In Proposition 8 below we will provide the explicit Poisson map relating the unreduced split variables (x o , x 1 , x 2 , x 3 ) to the reduced ones (a, v, M 1 , M 2 , M 3 ).

Regularizing the reduced systems directly?
Before moving on to relating the unreduced and reduced system, we mention the following heuristic argument. Take the generic situation that x 2 (t o ) = 0 when x 1 (t o ) = 0. Then This has a dramatic consequence in our reduction: e 1 = x 1 /|x 1 | sudenly flips from −x 2 (t o )/|x 2 (t o )| to +x 2 (t o )/|x 2 (t o )|. In this case we argue that the continuation of (35) beyond t o could be done using the velocity reversal symmetry of the problem.
We expect that this "infinite infinitesimal rotation" around e 3 will amount to an instantaneous rotation by π of the tangent plane. In other words, the vectors e 1 , e 2 will instantaneously change sign at t o , that is, Indeed, in view of (52) and we posit that κ g should be π times the delta function at s o = s(t o ).
We need to compute the composition The symplectomorphism (i) is the ψ ∇ of section 2.1, with underlying change of variables given in proposition 3 (in the particular case of the (n = 2)-sphere). The symplectomorphism (ii) is induced by the diffeomorphism (24) between T S 2 − 0 and M = R + × SO(3) (see also appendix B). The Poisson map (iii) is just the projection which forgets the R ∈ SO(3) (once we use left trivialized covectors: A straightforward tracking the above maps yields our final result (for simplicity we took r = β = 1).
Proposition 8. The above Poisson map taking the split variables (x, v, p, α) for T * (T S 2 ), as described in proposition 3, to the reduced variables (a, v, M 1 , M 2 , M 3 ) ∈ M red = T * R + × R 3 is given by

Comments and further questions
Control systems on anchored vector bundles More generally than in (1), which is a control problem with state space A = T Q, one could consider control problems with state space A (x, a) being a vector (or, more generally, affine) bundle q : A → Q with a connection ∇, and state equations of the forṁ with ρ : A → T M a given ('anchor') map. In this paper, we have been considering the particular case A = T Q and ρ = id. Also note that for u = 0 (uncontrolled problem) we recover the geodesic equations relative to (ρ, ∇). Examples of such systems arise from nonholonomic control problems [21] and control on (almost) algebroids [55] (for background on algebroids, see also [69], [100], [35]). We observe that control problems with state space T Q with a Levi-Civita connection can be recast, via the dual connection, to a control problem with state space T * Q. This is a usefull observation for landmark splines, since the problem is best described in terms of a cometric.
In general, the PMP leads to the cotangent bundle T * A and to the problem of finding the optimal hamiltonian. The connection ∇ allows us to obtain split variables generalizing those of section 2.1, Proposition 1 generalizes to this more general setting and provides a formula for the symplectic structure in split variables (also containing curvature terms). Moreover, the optimal hamiltonian is easy to find in these variables just as in 2.1. This will be detailed in [7].

Higher order splines and natural curvatures
For future work one can think of higher order splines as a control problem with state equation corresponding to The state space can be taken to be A = J k Q, the space of k-jets on Q. The cost functional can depend on u and (possibly) lower order covariant derivatives D (i) qq , i = 1, · · · k. An optimal curve γ(t) should connect two prescribed k-jets j (k) | xo , j (k) | x 1 (in computational anatomy it is often required that γ(t) passes through a number of intermediary points at prescribed times). In this paper, we have treated the case k = 1.
We remark that such a system is related to two very nice papers by Gay-Balmaz, Holm, Meier, Ratiu, and Vialard [46,47] on higher order lagrangians.
Our general approach via PMP applies here as well. One obtains T * (J k (Q)) and has to find the optimal hamiltonian. At this point, one notices that J k (Q) fits into a tower of affine bundles see [65]. One can thus find split variables for T * (J k (Q)) recursively using a given connection on J 1 Q = T Q. The higher order curvatures of the underlying curve γ(t) ∈ Q will be related to components of the control, similarly to what we obtained in appendix B for Q a surface (and k = 1) and, more generally, to what happens in elastica (see e.g. [64], [56,57]). We plan to pursue this in a sequel paper.

Is accessibility an issue?
For the general context of accessibility in mechanical control problems, see [9]. In Alan Weinstein's Ph.D. dissertation [99], about cut and conjugate loci on Riemannian manifolds, there is basic lemma stating that, if the manifold is complete, connected, and of finite volume, then any two unit tangent vectors can be joined by a smooth curve, parametrized by arc length, with geodesic ends, having geodesic curvature smaller than any arbitrarily small bound. Using this result it is easy to show accessibility for the time-minimal, bounded acceleration problem. We wonder if the same is true in higher order, namely joining two given 2-jets by a curve with arbitrarily small "jerk".

Controllability on vector bundles
In the seminal paper by Lewis and Murray [67] on configuration controllability of mechanical systems, the concept of symmetric product of vector fields was introduced. Their results were extended to mechanical systems with constraints and symmetries, see [31]. For a geometric interpretation, see [8]. Can the techniques be used in the general context of control problems on vector bundles with connection? Note that the control appears in a fraction rank(A)/(n + rank(A)) of the equations (further, the system can be sub-actuated). For results on controlabillity of affine connection mechanical systems, see [10].

Diffusion PCA
More generally, on any framework where the phase space is a cotangent bundle T * P of a manifold with a bundle structure P → B (and a connection), a splitting of variables will be useful. In the case of principal bundles G → P → P/G, the reduction of T * P goes back to Kummer [62,63] in the 1980's. For instance, a theory for diffusion principal component analysis (PCA) was developed by Sommer [93], based upon stochastic development via Eells-Elworthy-Malliavin construction of Brownian motion [38], [94]. A Hamiltonian system on the cotangent bundle of the frame bundle T * (Fr(Q)), governs the most probable paths 9 .

Interpreting the terms in the (simple splines) Hamiltonian equations
Peter Michor observed at the workshop that the extra terms in the equation forṗ i in Proposition 1 could be related to the concepts of symmetrized force and shape stress in his work with Michelli and Mumford [73].
For certain applications, Noakes has argued that L ∞ could be better than L 2 . Indeed, from the mathematical side, a drawback of cubic splines is that for manifolds of negative curvature the velocity can become infinite in finite time 9 A code is available in https://github.com/stefansommer. [82]. One can anticipate this behavior from equations (13). They contain curvatures -signs matter. In contradistinction, under bounded acceleration constraint, the scalar velocity grows at most linearly, so in all cases trouble is avoided by default. Would that be physiologically reasonable? Some simple experiments with n = d = 1 shows that, for cubic splines, the acceleration can attain high values of during the prescribed time interval, while the time minimal bounded acceleration can do the job in not a much longer time, depending on the concrete problem at hand. For robotics applications, or for an athlete, disastrous consequences could happen if the norm of the control force exceeds a given bound at some instant, see [96], [92].

Singular reduction
SO (3) does not act freely on T S 2 : trouble happens when the scalar velocity v vanishes (ie., the zero section of T S 2 ), and this propagates to non-freeness of the action of SO(3) on the symplectic manifold T * (T S 2 ). More generally, one may consider a vector bundle A → Q with a G-action, that is not free on the base (hence on the zero section). A procedure to do the singular Hamiltonian reduction of T * A is in order, a research direction that we hope to address in the future 10 .

Applying Morales-Ramis theory
Let us go back to cubic splines in S 2 as described in section 3. Since the kissing circles unstable periodic orbits are explicitly known, one may hope to prove nonintegrability using the Morales-Ramis approach [74]. However, linearizing (23) and doing the required Galois theory for the time periodic linear equations would be, no doubt, a tour-de-force. On could also attempt to show nonintegrability linearizing around the equator solutions.

Controlled Lagrangians
In [20,19] the concept of controlled Lagrangian is introduced, for mechanical systems L = T − V of natutal type. The control forces here keep the conservative nature of the controlled system. This is achieved by conveniently shaping the kinetic and/or potential energy. The modified system is still a closed-loop system, and the controlled system is Lagrangian by construction. Energy methods are used to find control gains that yield closed-loop stability. It would be interesting to see if such methods could be used to match tangent vectors.

Splines in infinite dimensional Riemannian geometry
This is a special edition about a meeting on infinite-dimensional Riemannian geometry, so we now try to link the present work to the infinite dimensional setting. As it is customary, the idea is to use the present study of systems with underlying finite dimensional configuration spaces Q as simplified models for the cases in which Q is an infinite dimensional Riemannian manifold. In this direction, the Levi-Civita connection and the curvature of Sobolev metrics on Q = Diff have been studied by several authors, see eg. [14], [72], [48], [59,6,60], [75]. Below we enumerate some related questions.
(i) PDEs for splines in shape space. This means optimal control problems with state space A = T Diff. In order to describe T (T Diff) and T * (T Diff), one can take advantage of the fact that Diff is a group and, thus, T Diff = Diff × X . What are the corresponding PDEs for L 2 and L ∞ splines? They should involve not only the momentum density m but another density p corresponding to a Pontryagin multiplier (alternatively, a PDE for m involving three time derivatives).

B State equations on convex surfaces and the Gauss map
In this appendix, we elaborate on a description of non-zero tangent vectors on a convex surface Σ which uses the Gauss map. We use it to provide an alternative form of the state equations (7) on T Σ. This description is used in section 3 in the particular case of Σ = S 2 to exploit the rotational symmetry. Let Σ ⊂ R 3 be a closed smooth convex surface. The Gauss map induces a diffeomorphism between T Σ − 0 and R+ × SO (3): where R ∈ SO(3) is constructed as follows: points q ∈ Σ correspond uniquely to external unit normal vectors to the surface, which we denote e3. Now, a nonzero tangent vector vq corresponds uniquely to a pair (v, e1) with vq = v e1 , v > 0 , and e1 · e3 = 0 , |e1| = |e3| = 1.
We use a redundant vector e2 = e3 × e1 to construct the matrix R with columns e1, e2, e3. Therefore, a control problem with state space T Σ corresponds to a control problem on SO(3) × R+, provided we exclude the zero section 11 . Let us now move on to rewritting the state equations (7) in our present situation. Recall the Darboux formulas for a curve γ(s) ∈ Σ ( = d/ds) e 1 = κg e2 + κn e3 , e 2 = −κg e1 + τg e3 , e 3 = −κn e1 − τg e2 where κg is the geodesic curvature, κn the normal curvature, and τg the geodesic torsion of γ. These formulas can be rewritten aṡ The normal curvature is not freely controllable since it corresponds to the force that constrains the curve to stay in the surface. Indeed, taking derivatives in the ambient space,γ =v e1 + v 2 e 1 =v e1 + v 2 (κg e2 + κn e3) = ∇γγ + v 2 κn e3 with κn = e 1 · e3 = −e 3 · e1 := B(e1, e1) where B is the second fundamental form of the surface. On the other hand, the intrinsic description of the state equations, using the Levi-Civita connection, reads ∇γγ = u1 e1 + u2 e2 (51) where u1, u2 are the controls. The previous simple calculation thus showed that u1 =v and u2 = v 2 κg.
But the geodesic torsion κg also admits the following interesting formula found by Darboux where φ is the angle between the unit tangent vector e1 to the curve and a principal direction on the surface. We then conclude that the state equations can be written aṡ v = u1 ,Ṙ = R X , X =      Note that v is growing quadratically with respect to a. The reconstructed trajectory is approaching a neighborhood of an equator. It remains to be seen if it stays there or returns to a vicinity of the reduced equilibrium.