Geodesic boundary value problems with symmetry

This paper shows how left and right actions of Lie groups on a manifold may be used to complement one another in a variational reformulation of optimal control problems equivalently as geodesic boundary value problems with symmetry. We prove an equivalence theorem to this effect and illustrate it with several examples. In finite-dimensions, we discuss geodesic flows on the Lie groups SO(3) and SE(3) under the left and right actions of their respective Lie algebras. In an infinite-dimensional example, we discuss optimal large-deformation matching of one closed curve to another embedded in the same plane. In the curve-matching example, the manifold $\Emb(S^1, \mathbb{R}^2)$ comprises the space of closed curves $S^1$ embedded in the plane $\mathbb{R}^2$. The diffeomorphic left action $\Diff(\mathbb{R}^2)$ deforms the curve by a smooth invertible time-dependent transformation of the coordinate system in which it is embedded, while leaving the parameterisation of the curve invariant. The diffeomorphic right action $\Diff(S^1)$ corresponds to a smooth invertible reparameterisation of the $S^1$ domain coordinates of the curve. As we show, this right action unlocks an important degree of freedom for geodesically matching the curve shapes using an equivalent fixed boundary value problem, without being constrained to match corresponding points along the template and target curves at the endpoint in time.

1. Introduction. In this paper we are concerned with finding geodesics between points on manifolds. The construction of geodesics is useful for studying problems on manifolds since they can describe the relationship between two points. Within a coordinate patch on a manifold, any point can described relative to a reference point by specifying a direction and a length along the geodesic in that direction. This becomes useful for performing statistics on the coordinate patch, for example. In this paper we consider problems in which the endpoint of the trajectory is only fixed up to the orbit of a Lie group. In low dimensional cases (and we shall describe some examples of these) it is often easy to solve these problems by constructing reduced coordinates which do not change under the action of the Lie group. However, in many cases it is difficult to construct such coordinates, especially if the problem is to be discretised and solved numerically. In this paper we provide a framework that allows one to work with full unreduced coordinates on the manifold, by transforming to an equivalent problem which has the endpoint of the trajectory fixed exactly.
There are many examples of problems where this framework can be applied, but we are motivated by the problem of obtaining diffeomorphisms on R 2 which map one embedded oriented curve Γ A into another embedded oriented curve Γ B (with the same orientation), and which minimise a given metric so that they are geodesics in the diffeomorphism group. The aim is to find a characterisation of curve Γ B with respect to curve Γ A that is independent of parameterisations of the curves. This means that we do not specify a priori the point on Γ A which gets matched to each specific point on Γ B , and so the minimisation is performed over all parameterisations of the curves. In practise the computation is performed using a particular parameterised curve q ∈ Emb(S, R 2 ) (where S is the embedded space, for example, the circle for simple closed curves). In computing the equations of motion, a conjugate momentum p q ∈ T * q Emb(S, R 2 ) is constructed, and the flow taking the initial curve Γ A to the final curve Γ B can be characterised entirely by the initial conditions p q | t=0 for the conjugate momentum. In fact, it turns out that p q | t=0 is normal to the curve, so the flow can be characterised by a one-dimensional signal. Since T * q Emb(S, R 2 ) is a linear space, linear statistics can be computed on p q | t=0 . For example, this may allow one to test the hypothesis that there is a statistical correlation between between the shape of the surface of a biological organ, obtained from a medical scan, and future development of disease.
To discuss the issues further, we formulate the curve matching problem described above, which may be regarded as an control problem in the sense of the problems discussed in [?, ?, ?]: Definition 1.1 (Curve matching problem). Let q(s; t) be a one-parameter family of parameterised simple closed curves in R 2 , with s ∈ [0, 1] being the curve parameter and t ∈ [0, 1] being the parameter for the family. Let u(x; t) be a one-parameter family of vector fields on R 2 . Let η be a diffeomorphism of S 1 . We seek q and u which satisfy [Final state (Target)] q(s; 1) = q B (η(s)), where · V is the chosen norm which defines the space of vector fields V .
The solution of this problem describes a geodesic in the diffeomorphism group which takes the simple closed curve Γ A parameterised by q A to the simple closed curve Γ B parameterised by q B . We represent the shapes of simple closed curves as elements of Emb(S 1 , R 2 )/ Diff + (S 1 ), where Diff(S 1 ) is the group of diffeomorphisms of S 1 and the subscript + indicates the proper subgroup of Diff(S 1 ) that is connected to the identity (i.e. the orientation preserving diffeomorphisms). However, we do not want to calculate on this space; instead, we want to calculate on the full space Emb(S 1 , R 2 ) by minimising over all reparameterisations η(s) ∈ Diff + (S 1 ). We require that q A and q B have the same topology, i.e. that there exists a diffeomorphism g ∈ Diff + (R 2 ) (with smoothness prescribed by the metric · V ) such that q A = g • q B . In particular, if we restrict to simple closed curves, this requires that q A and q B both traverse the curves in the same direction (i.e. both clockwise or both anti-clockwise).
There are two general strategies for solving such problems. The first strategy, used for example in [?], is to use a gradient method (i.e. a modification of the steepest descent method such as the nonlinear conjugate gradient method [?, and references therein]) to minimise the action integral over paths q(s, t) which satisfy the dynamical constraint (this constraint was enforced "softly" via a penalty term in [?]). An alternative method, referred to in [?] as the "Hamiltonian method", is to introduce Lagrange multipliers p(s, t) which enforce the dynamical constraint, and to derive Hamilton's canonical equations for q and p q , following the general derivation described in [?, for example]. Minimisation over the reparameterisation η, together with a conservation law obtained from Noether's theorem, results in the condition that the tangential component of p q vanishes. The aim of the Hamiltonian method is to turn an optimisation problem into an algebraic equation given by the time-1 flow map of Hamilton's canonical equations. One then solves a shooting problem to find initial conditions for the normal component p q which generate solutions to Hamilton's equations that satisfy the boundary condition (3). The difficulty in solving this problem numerically lies in finding a good numerical discretisation of the target constraint condition (3). Various functionals have been proposed which vanish when the constraint condition is satisfied. In [?] a functional was proposed based on singular densities (measures), and in [?] a functional was proposed based on singular vector fields (currents). An alternative spatial discretisation for the current functional based on particle-mesh methods was proposed in [?]. There are several difficulties with these functionals: one is that after numerical discretisation the functionals do not vanish at the minima, and the boundary condition must be replaced by a functional minimising condition. It is also difficult to express the probability distribution of the functional given the distribution of measurement errors; this is important for statistical modelling.
In this paper we consider a transformation of problems of the above type, which results in an alternative formulation that removes the reparameterisation variable η from the target constraint, thereby resulting in a standard two point boundary value problem on T * Emb(S, R 2 ) (with a constraint on the initial conditions plus an additional parameter). This transformation can be applied to a very general class of problems; so we present it in the general case of Lie group actions on a manifold.
The rest of this paper is organised as follows. In Section 2, we formulate the optimal control problem, then transform to the geodesic problem with symmetry and prove that the two problems are equivalent. In Section 3 we give some examples and discuss the application to matching curves and surfaces. Section 4 is the summary and outlook.

Reparameterised geodesic boundary value problems with symmetry.
In this section we describe a general framework for geodesic boundary value problems with symmetry. We define the following Optimal Control Problem.
Definition 2.1 (Geodesic boundary value problem with symmetry). Let Q be a manifold, let G be a Lie group acting on Q from the left, and let H be a (possibly different) Lie group acting on Q from the right that commutes with the left action of G on Q, with corresponding Lie algebras g and h, and corresponding Lie algebra actions X G and X H respectively. Furthermore, let A : g → g * be a positive-definite self-adjoint operator and let ·, · g : g × g * → R be a nondegenerate pairing which defines an inner product on g. We seek a one parameter family q of points on Q parameterised by t ∈ [0, 1], a one parameter family ξ of elements of g for t ∈ [0, 1], and η ∈ H, which minimise where q A , q B are chosen points on Q related by some element g of G so that q A = gq b , and where R η is the right-action of η on Q.
This problem is an example of a Clebsch optimal control problem, as studied from a geometric point of view in [?]. In this case we seek the shortest path in Q from q A to any point q B η, η ∈ H. This means we are seeking the shortest path in Q/H, but are performing the computation on Q. In many cases it is much easier to compute on Q, for example when Q is a vector space. We refer to this process of solving a problem on Q/H by calculating on Q as "un-reduction". One approach to solving this problem is to derive equations of motion for q, ξ and an optimal condition for η and then solving a shooting problem to find η and the initial conditions for ξ which allow equation (6) to be satisfied. We can derive the equations of motion by enforcing the reconstruction relation (4) as a constraint using Lagrange multipliers p q ∈ T * q Q. This approach leads to the following variational principle.
Definition 2.2 (Variational principle for geodesic boundary value problem with symmetry). We seek (p, q) ∈ T * Q and ξ ∈ g for t ∈ [0, 1], and η ∈ H, which satisfy subject to where we allow p q , q, ξ and η to vary.
From this variational principle we can derive the equations of motion, which can be used in solving the shooting problem. Before we do this, we recall the definition of the cotangent-lifted momentum map: Definition 2.3. Given an action of a Lie algebra g on Q, the cotangent-lifted momentum map J : T * Q → g is defined from the formula for all ζ ∈ g. Since we have two Lie algebra actions, we shall write J G for the cotangent-lifted momentum map corresponding to the left action X G of g on Q, and J H for the cotangent-lifted momentum map corresponding to the right action X H of h on Q.
Lemma 2.4 (Equations of motion for geodesic problem). At the optimum, the following equations are satisfied (weakly, for appropriate pairings): Proof. The proof is a direct calculation.
Lemma 2.4 reformulates the geodesic calculation as a shooting problem in which one seeks initial conditions for p q such that q| t=1 = R η q B where η is fixed by the condition (13) which arises from minimising over η and ensures we have the shortest path over Q/H.
Next we show conservation of the momentum map J H ; this will enable us to transfer condition (13) from t = 1 to t = 0.
Proof. The problem in Definition 2.2 is invariant under transformations q → R α q, α ∈ H, which are generated by γ ∈ h, and hence the corresponding momentum map J H (p q ) is conserved, by Noether's theorem. Since it vanishes for t = 1, by equation (13), it must vanish for all times.
Lemma 2.5 states that solutions of the optimal control problem all have vanishing right action momentum map J H (p q ) = 0. This is what facilitates the "unreduction". Namely, we can compute on Q instead of Q/H by keeping J H (p q ) = 0. To obtain the shortest path between two points in Q/H by solving in Q, select a point q ∈ Q which is a member of the equivalence class which is the initial point in Q/H, and find initial conditions for p q such that J H (p q ) = 0; so that the solution to equations (10-12) satisfies q| t=1 = R η q for some η ∈ H. Computationally, there are reasons why solving the problem in this form may be difficult. In Section 3.3, we shall describe how the difficulty arises for the curve matching problem specified in the Introduction. In this paper, we shall introduce a reformulation of the problem for which there is a single fixed value of q| t=1 .
Before introducing the reformulation, we briefly discuss the reduced equation for the Lie algebra variable ξ ∈ g. The latter is the Euler-Poincaré equation for Hamilton's principle with Lagrangian given by the energy ξ, Aξ g /2, where A : g → g * is the positive-definite self-adjoint operator in Definition 2.1 of the geodesic matching problem.
Lemma 2.6 (Reduced equation for geodesic problem). The Lie algebra variable ξ for the geodesic matching problem stated in Definition 2.1 satisfies where we define the operation ad : g × g → g and its dual ad The proof follows from direct computation of the time-derivative of γ, Aξ g and substitution of equations (10-12). See for example, [?, ?].
We will next define a modification of the problem stated in Definition 2.1, which has the advantage that the endpoint conditions do not contain a free reparameterisation variable. This reformulation is more amenable when solving the curve matching problem numerically, for example. We shall proceed to show that solutions of the modified problem can be transformed into solutions of the problem stated in Definition 2.1.
Definition 2.7 (Reparameterised geodesic problem with symmetry). Let Q be a manifold, let G be a Lie group acting on Q from the left, and let H be a (possibly different) Lie group acting on Q from the right that commutes with the left action of G, with corresponding Lie algebras g and h, and corresponding Lie algebra actions X G and X H respectively. Furthermore, let A : g → g * be a positive-definite selfadjoint operator. We seek a one parameter family q of points on Q for t ∈ [0, 1], a one parameter family ξ of elements of g for t ∈ [0, 1], and ν ∈ h, which minimise where ·, · g is the usual inner product on g , subject to the constraints and q A , q B are chosen points on Q.
Note that in this modified definition, we do specify the final boundary condition for q without allowing arbitrary symmetry transformations using H. However, we also introduce an additional variable ν which moves q in the direction of symmetries generated by h. We shall derive the equations of motion associated with this modified problem, and the associated conservation laws. These will lead us to conclude that it possible to construct solutions of the problem in Definition 2.1 out of solutions of the problem in Definition 2.7, and the latter can be solved as a shooting problem in which the boundary conditions are explicitly specified, rather than as an algebraic condition. As before, we can derive the equations of motion forq, ξ and the condition for ν by enforcing the reconstruction relation (15) as a constraint using Lagrange multipliersp q ∈ T * q Q. We then use variational calculus to obtain the equations of motion, as described in the following lemma.
Lemma 2.8 (Equations of motion for reparameterised geodesic problem). At the optimum, the following equations are satisfied in the sense of appropriate pairings: Furthermore, Proof. Similar to the proof of Lemma 2.4.
Proceeding as before, we can transform (21) into an initial condition by making use of the conservation of the right-action momentum map, J H . Lemma 2.9 (Vanishing momentum for reparameterised geodesic problem). The system of equations (19-20) has conserved momentum J H = 0.
Proof. The problem in Definition 2.7 is invariant under transformations q → qα, α ∈ H, which are generated by γ =∈ h, and hence J H is conserved by Noether's theorem. Lemma 2.8 gives 0 = 1 0 J H dt = J H . Next we note that ξ obtained from Definition 2.7 satisfies the same reduced Euler-Poincaré equation as ξ obtained from Definition 2.1.
Proof. Similar to the proof of Lemma 14, noting that the actions X G and X H commute.
This means that we can show that the two problems produce equivalent solutions provided that the initial conditions for ξ are the same in both cases. The following theorem establishes this result.
Proof. First we take ξ and ξ from equations (12) and (20) respectively, and show that ξ = ξ. Since the left and right actions commute, we have that and hence ξ = ξ. Next we verify the equations for q and p q . Taking the time derivative of q, we have and so as required. To check the time evolution equation for p q , we take the inner product with an arbitrary tangent vector v ∈ T q Q, to find: as required. It remains to check the boundary conditions. Trivially, q| t=0 = q A , q| t=1 = q| t=1 ψ| −1 t=1 = q B η, as required. Finally, we need to check the end condition (13). From Corollary 2.9, we have J H (p q )| t=0 = J H (p q )| t=0 = 0, and Lemma 2.5 implies that J H (p q )| t=1 = 0. Hence the boundary conditions are satisfied.
3. Examples. In this section we describe examples of the reparameterised geodesic problem with symmetry, and discuss its applications to the characterisation of the shape of curves and surfaces.
3.1. Example: SO(3). We illustrate our results with the case of the action of SO(3) on itself which gives rise to the equations of a rotating rigid body. We consider the problem in which the end point boundary condition is only determined up to a rotation of the rigid body about its z-axis. Of course, this problem can also be solved by picking reduced coordinates, but we use it as here as a simple example.
Definition 3.1 (Optimal control of a symmetric rigid body). For a given symmetric matrix I, we seek Q(t) which is a one-parameter family of matrices SO (3), ω(t) which is a one-parameter family of matrices in so (3), and R θ which is a rotation in the z-axis by an angle θ, which satisfy This problem is an example of the optimal control problem in Definition 2.1, with the manifold Q being SO(3), the group G being SO(3) acting from the left, and the group H being SO(2) acting from the right. We identify q ≡ Q, ω ≡ ξ, R θ ≡ η, A ≡ I and p q ≡ P, where P ∈ T * SO(3) is the conjugate momentum to Q. Application of Lemmas 2.4 and 2.5 gives the following equations: which are the equations for a rotating rigid body with vanishing z-component of the angular momentum. For this problem, obtaining a solution is simple, since one can define coordinates on T (SO(3)), and remove the coordinates associated with the R θ direction and the corresponding vanishing conserved momentum, and solve a two-part boundary problem for the remaining coordinates. However, we wish to develop a methodology for numerical discretisations of infinite-dimensional problems where it is less clear how to do this. Hence, we define the following formulation which makes use of a time-varying "relabelling" transformation in the R θ direction. Theorem 2.11 states that to obtain solutions to equations (26), we can alternatively solve the following modified problem: Definition 3.2 (Reparameterised optimal control of a symmetric rigid body). We seek Q(t) which is a one-parameter family of matrices in SO(3), ω(t) which is a one-parameter family of matrices in so(3) and ν which is a generator of rotations about the z-axis, which satisfy min ω,ν 1 0 1 2 ω, Iω so(3)×so(3) * dt subject to the constraints where I is a chosen symmetric matrix.
Lemmas 2.8 and 2.9 states that the solution to this problem satisfies the following equations:Q − ωQ + Qν = 0,Ṗ + ω T P − P ν T = 0, Hence, to obtain a solution to equations (26) with boundary conditions (24) and (25) we solve the two-point boundary value problem given by equations (30) with boundary conditions (28)(29). This can be formulated as a shooting problem, in which we seek ν and initial conditions for P with vanishing z-component of angular momentum, such that the end point boundary condition (29) is satisfied. We then construct the reparameterisation matrix R(t) from R(t) = exp(νt), and use equation (22) to reconstruct the solution in the form: 3.2. Example: SE(3). We next describe the example of the action of SE(3) on itself from the left, with SO(2) acting from the right. This example could describe a docking problem of a spacecraft onto a space station. The spacecraft can apply torque to rotate about a central point, or can apply thrust to move itself in the direction in which it is pointing, and we wish to dock the spacecraft using minimal energy. In the language of image registration, this is known as rigid registration. We consider the problem in which the end point boundary condition is only determined up to a rotation of the rigid body about its z-axis. In the spacecraft analogy, this corresponds to a docking procedure which does not require the spacecraft to have any particular orientation about the z-axis when docking. As in the previous example, this problem can also be solved by picking reduced coordinates, but it serves as a prototype for infinite dimensional problems.
Following the notation of [?], we represent the elements q ∈ SE(3), p q ∈ T * q SE(3), ξ ∈ se(3) as 4 × 4 matrices where ω is an antisymmetric matrix, and v ∈ R 3 . The energy cost (metric functional) and reconstruction relation are then given by The start and end conditions for the orientation and poistioin are specified as where R θ is a rotation through any angle θ about the z-axis. If we solve the problem in Definition 2.1, then Lemmas 2.4 and 2.5 give the equationṡ where w is as defined in equation (26). From Lemma 2.5, this quantity vanishes for all t. Theorem 2.11 then states that a solution to these equations can be obtained by solving the following reparameterised equations: with This gives a two-point boundary value problem with a constraint on the initial conditions plus an extra parameter, which can be solved as a shooting problem by finding α and P (subject to the angular momentum constraint) such that Q and r reach their target values Q B and r B . A solution to the problem in Definition 2.1 can then be reconstructed by defining R(t) = exp(αwt), and using the following formulae: 3.3. Curve matching. In this section we return to the problem described in Definition 1.1, and discuss a number of practical issues which are addressed by the formulation discussed in this paper. The aim of solving the problem is to find a characterisation of the simple closed curve Γ B in terms of the reference simple closed curve Γ A , together with a scalar periodic function p(s) which specifies the initial conditions for the normal component of the conjugate momentum p(s; t). In this case, Q is the space Emb(S 1 , R 2 ) of functions q : S 1 → R 2 , G is the group Diff + (R 2 ) of diffeomorphisms of R 2 which acts on Q from the left: Φ L (g, q)(s) = g(q(s)), ∀s ∈ S 1 , and H is the group Diff + (S 1 ) of diffeomorphisms of S 1 which acts on Q from the right: Φ R (g, q)(s) = q(g(s)), ∀s ∈ S 1 . The left and right actions can be expressed succinctly as It is clear from this that the actions of G and H commute with each other. Lemma 2.4 gives the dynamical equations ∂ ∂t q(s; t) = u(q(s; t), t), ∂ ∂t p(s; t) = −(∇u(q(s; t), t)) T · p(s; t) , where ·, · X is the inner product on the vector fields X(Ω) associated with the norm · X , for any test function w ∈ X * . This equation for u has the weak solution which is the singular-solution momentum map, J Sing discussed in [?]. The end condition is p(s; 1) · ∂ ∂s q(s; 1) = 0, ∀s ∈ S 1 , and Lemma 2.5 states that this conserved momentum vanishes for all t. This corresponds to p(s; t) being normal to the curve parameterised by q(s; t). Hence, to find geodesics between Γ A and Γ B , we solve a shooting problem and seek initial conditions p(s; 0) with vanishing tangential component, which fix q(s; 1) = q B (η(s)) for some η ∈ Diff + (S 1 ). Having solved this problem, one can describe Γ B entirely in terms of p(s) = p(s; 0) · n(s) where n is the normal to Γ A . The solution to the problem also provides the distance between the two curves.
There are a number of difficulties with solving such a shooting problem numerically. The parameterisation of the curve must necessarily be discretised numerically, typically by a list of points, as in [?, ?], which can be obtained from a piecewiseconstant representation of q as a function of s [?], or as piecewise linear geometric currents [?]. Having taken the discretisation, the reparameterisation symmetry is broken (although a remnant of it is left behind, as described in [?]) which means that it is difficult to obtain a discrete form of the end condition q(s; 1) = q B (η(s)). As described in the Introduction, this problem has been approached by proposing various functionals which are minimised when q(s; 1) and q(s) overlap the most. However, these functionals produce quite a complicated landscape with local minima, and the case of studying large deformations we have found that they can result in odd artifacts in the shooting process. Also, the changes in these functionals with respect to measurement error are quite technical to quantify which makes statistical inference more complicated. Another difficulty is that of adaptivity. As illustrated in Figure 1, constraining p to be normal to the curve means that any local large deformations give rise to large amounts of stretching which then results in loss of accuracy in the approximation of the functional used to enforce the end condition for q. One possible way to avoid this is to adaptively refine the grid point density in the initial curve during the shooting process.
These difficulties are all removed if, instead, one solves the following problem: Definition 3.3 (Reparameterised curve matching problem). Let q(s; t) be a oneparameter family of parameterised curves in R 2 , with s ∈ [0, 1] being the curve parameter and t ∈ [0, 1] being the parameter for the family. Let u(x; t) be a oneparameter family of vector fields on R 2 . Let ν be a vector field on S 1 . We seek q, u and ν which satisfy where · V is the chosen norm which defines the space of vector fields V .
Lemmas 2.8 and 2.9 state that extrema of this problem can be obtained by solving the shooting problem ∂ ∂t q(s; t) = u(q(s; t), t) + ν(s) ∂ ∂s q(s; t), ∂ ∂t p(s; t) = −p(s; t) · ∇u(q(s; t), t) − ∂ ∂s (ν(s)p(s; t)) , w, u(·, t) V = S 1 p(s; t) · w(q(s; t), t) ds, ∀w ∈ V * , Figure 1. Figure illustrating the way in which deformation takes place in the parameterisation-independent geodesic equations for embedded curves. The initial curve is shown on the left, and the final curve is shown on the right. Since the momentum is constrained to be normal to the curve, and since the velocity is obtained by applying a smoothing kernel to the momentum, the change in the shape emerges as local stretching of the curve, and the discrete points defining the shape become separated.
The aim is to find ν(s) and normal components of p(s; 0) such that these boundary conditions are satisfied. Note that in this modified problem, the boundary condition for q(s; 1) is specified pointwise (i.e. there is no reparameterisation variable η in the boundary condition). This means that the error can be described directly in terms of q(·; 1) − q B 2 for some chosen norm, which can be discussed in terms of measurement errors directly. Theorem 2.11 then states that a solution to the problem described in Definition 1.1 can be reconstructed via q(s; t) = q(η(s; t), t), p(s; t) = p(η(s; t), t) ∂ ∂s η(s; t), η(s; t) = exp(ν(s)t).
This transformation produces a equivalent shooting problem in which the end condition for q is now fixed. See figures 2 and 3 for an example of a solution to this problem. The numerical discretisation and solver methods will be benchmarked and analysed in a future publication.

Figure 2.
Results from a matching algorithm based on the reparameterised equations. The aim is to find the initial conditions for p (with vanishing tangential component) and reparameterisation vector field ν which maps the initial curve (in this case a circle) onto a set of points observed from another curve. The variables p · n and ν are obtained by minimising a functional consisting of the Euclidean distance between certain points on the final curve plus a regularising function which enforces smoothness of p and ν namely the L 2 norm of p · n and the H 3 norm of ν. The reparameterised equations are solved using a splitting method in which the left-action part is discretised using the particle-mesh method described in [?], and the reparameterisation is discretised using a semi-Lagrangian method. The displayed plots are: (Top-Left) the template curve and the final curve together with the points from the target curve, (Top-Right) the momentum co-vectors plotted on the template curve, (Bottom-Left) the vector field on S 1 used to generate the reparameterisation, and (Bottom-Right) the final reparameterisation η plotted against the identity map for reference. The numerical solution was obtained in about three minutes after about 300 iterations of the iterative solver, based on the BFGS algorithm.
4. Summary and outlook. In this paper, we studied an optimal control problem on a Lie group in which the end boundary condition is specified only up to a symmetry. We showed how this problem can be transformed into a modified problem in which the end boundary condition is fixed, but an extra parameter is introduced in the reconstruction relation, and proved that the two problems are equivalent. This  Figure 2. On the left, the template curve is shown together with a reference grid (the grid is purely for visualisation and is not part of the numerical method). On the right, the final curve is shown together with the action of the diffeomorphism on the reference grid.
approach is motivated by the problem of computing a geodesic on the diffeomorphism group in the plane which takes one curve to another. The transformation gives rise to a system of equations for a parameterised curve in which the end boundary condition for each value of the parameter is fixed. This means that when a discrete approximation of the curve is used to solve this problem numerically, the end boundary condition can again be fixed exactly. In particular, when solving a shooting problem, this means that the error between the computed curve and the target curve can be computed simply by measuring the Euclidean distance between points for each value. This method extends straightforwardly to the problem of finding geodesics in the three-dimensional diffeomorphism group which take one parameterised surface to another, with the end boundary condition specified only up to reparameterisations of the target surface. This problem has many applications in, for example, biomedicine, since it allows topologically equivalent surfaces to be characterised by a momentum field distributed on the template surface. Such momentum fields exist in a linear space and so can be manipulated using linear techniques and still a topologically equivalent surface will always be obtained. In the standard approach to planar image registration, the problem of registering a specified closed curve (called the template) at time t = 0 to another (the target) at time t = 1 amounts to deforming the space R 2 in which the template curve is embedded until it matches the fixed target image to within a certain tolerance. Here, we considered a manifold Q comprising the space of closed curves S 1 embedded into the plane R 2 , written as Q = Emb(S 1 , R 2 ). There are two Lie group actions available for manipulating the closed curves in this description. The action G×R 2 → R 2 of the Lie group G = Diff + (R 2 ) by composition from the left deforms the range space R 2 , and thereby drags along a curve embedded in it as GQ = Emb(S 1 , G·R 2 ). This left action produced the singular-solution momentum map, J Sing in equation of an advection term in the q equation, and a continuity term in the p equation, which are very well developed in the finite volume approach. One can also utilise the fact that the left-and right-actions commute by applying a splitting method (so that steps of left-and right-action are alternated) and them use a semi-Lagrangian method for the right-action steps, as was used in the numerical examples in this paper. The commuting actions allow great flexibility and one can even perform all the right-action steps first before applying the left-action steps. We will also investigate discontinuous higher-order polynomial representations of p and q which lead to discontinous Galerkin methods. Since the error in the transformed problem can be quantified in terms of the Euclidean distance between points on the curve for each parameter value, the reparameterised formulation also makes it much easier to perform Bayesian studies in which one observes points on a curve with observation error from some probability distribution, and then one attempts to estimate the probability distribution for the initial conditions of p (and ν) for which specified points on the curve match the actual position of the observed points, after applying the time-1 flow map of the transformed geodesic equations for p and q. This approach provides a considerably simplified observation operator to which algorithms such as the Monte Carlo Markov Chain algorithm can be applied.