Sobolev metrics on shape space of surfaces

Let $M$ and $N$ be connected manifolds without boundary with $\dim(M)<\dim(N)$, and let $M$ compact. Then shape space in this work is either the manifold of submanifolds of $N$ that are diffeomorphic to $M$, or the orbifold of unparametrized immersions of $M$ in $N$. We investigate the Sobolev Riemannian metrics on shape space: These are induced by metrics of the following form on the space of immersions: $$ G^P_f(h,k) = \int_{M} \g(P^f h, k)\, \vol(f^*\g)$$ where $\g$ is some fixed metric on $N$, $f^*\g$ is the induced metric on $M$, $h,k \in \Gamma(f^*TN)$ are tangent vectors at $f$ to the space of embeddings or immersions, and $P^f$ is a positive, selfadjoint, bijective scalar pseudo differential operator of order $2p$ depending smoothly on $f$. We consider later specifically the operator $P^f=1 + A\Delta^p$, where $\Delta$ is the Bochner-Laplacian on $M$ induced by the metric $f^*\bar g$. For these metrics we compute the geodesic equations both on the space of immersions and on shape space, and also the conserved momenta arising from the obvious symmetries. We also show that the geodesic equation is well-posed on spaces of immersions, and also on diffeomorphism groups. We give examples of numerical solutions.


Introduction
Many procedures in science, engineering, and medicine produce data in the form of shapes. If one expects such a cloud to follow roughly a submanifold of a certain type, then it is of utmost importance to describe the space of all possible submanifolds of this type (we call it a shape space hereafter) and equip it with a significant metric which is able to distinguish special features of the shapes. Most of the metrics used today in data analysis and computer vision are of an ad-hoc and naive nature; one embeds shape space in some Hilbert space or Banach space and uses the distance therein. Shortest paths are then line segments, but they leave shape space quickly. Riemannian metrics on shape space itself are a better solution. They lead to geodesics, to curvature and diffusion. Eventually one also needs statistics on shape space like means of clustered subsets of data (called Karcher means on Riemannian manifolds) and standard deviations. Here curvature will play an essential role; statistics on Riemannian manifolds seems hopelessly underdeveloped just now.
1.1. The shape spaces used in this work. Thus, initially, by a shape we mean a smoothly embedded surface in N which is diffeomorphic to M . The space of these shapes will be denoted B e = B e (M, N ) and viewed as the quotient (see [18] for more details) More generally, a shape will be an element of the Cauchy completion (i.e., the metric completion for the geodesic distance) of B i (M, N ) with respect to a suitably chosen Riemannian metric. This will allow for corners. In practice, discretization for numerical algorithms will hide the need to go to the Cauchy completion.

1.2.
Where this work comes from. In [20], Michor and Mumford have investigated a variety of Riemannian metrics on the shape space B i (S 1 , R 2 ) = Imm(S 1 , R 2 )/ Diff(S 1 ) of unparametrized immersion of the circle into the plane. In [19, section 3.10] they found that the simplest such metric has vanishing geodesic distance; this is the metric induced by L 2 (arc length) on Imm(S 1 , R 2 ): In [18] they found that the vanishing geodesic distance phenomenon for the L 2metric occurs also in the more general shape space Imm(M, N )/ Diff(M ) where S 1 is replaced by a compact manifold M and Euclidean R 2 is replaced by Riemannian manifold N ; it also occurs on the full diffeomorphism group Diff(N ), but not on the subgroup Diff(N, vol) of volume preserving diffeomorphisms, where the geodesic equation for the L 2 -metric is the Euler equation of an incompressible fluid. In [20, sections 3, 4 and 5] three classes of metrics were investigated: Almost local metrics on planar curves, Sobolev metrics on planar curves, and metrics induced from Sobolev metrics on the diffeomorphism group of the plane. The results about almost local metrics from [20, section 3] were generalized by the authors to the case of surfaces in [2]. Now we take up the investigations from [20, section 4]. The immersion-Sobolev metric considered there is The interesting special case p = 1 and A → ∞ has been studied in [22,25] and in [24] where an isometry to an infinite dimensional Grassmannian with the Fubini-Study metric was described. In this case, the metric reduces to: The cases p = 1, 2 and A → ∞ have also been treated in [15], where estimates on the geodesic distance are proven and the metric completion of the space of curves is characterized.
In this work we generalize the immersion-Sobolev metrics from [20, section 4] to higher dimensions and to non-flat ambient space, namely to the shape space B i (M, N ) = Imm(M, N )/ Diff(M ) of surfaces of type M in N ; here M is a compact orientable connected manifold of smaller dimension than N , for example a sphere S m , m < dim(N ).
metric. The G 0 -horizontal vectors in T f Imm(M, N ) are just those vector fields along f which are pointwise g-normal to f (M ); we will call them normal fields.
All of the metrics we will look at will be of the form (see section 6): where P f : T f Imm → T f Imm is a positive bijective operator depending smoothly on f , which is selfadjoint unbounded in the Hilbert space T f Imm with inner product G 0 f . We will assume that P is in addition equivariant with respect to reparametrizations, i.e.
The G P -horizontal vectors will be those h ∈ T f Emb(M, N ) = C ∞ (M, N ) such that P f h is normal.
The tangent map of the quotient map Emb(M, N ) → B i (M, N ) is then an isometry when restricted to the horizontal spaces, just as in the finite dimensional situation. Riemannian submersions have a very nice effect on geodesics: the geodesics on the quotient space B i are exactly the images of the horizontal geodesics on the top space Imm; by a horizontal geodesic we mean a geodesic whose tangent lies in the horizontal bundle. The induced metric is invariant under the action of Diff(M ) and therefore induces a unique metric on B i . See for example [2, section 1]. Later in section 8 we shall consider the special case P f = 1 + A∆ p .
1.4. Inner versus outer metrics. The metrics studied in this work are induced from Imm(M, N ) on shape space. One might call them inner metrics since the differential operator governing the metric is defined intrinsic to M . Intuitively, these metrics can be seen as describing some elastic or viscous behaviour of the shape itself.
In contrast to these metrics, there are also metrics induced from Diff(N ) on shape space. (The widely used LDDMM algorithm uses such a metric.) The differential operator governing these metrics is defined on all of N , even outside of the shape. Intuitively, these metrics can be seen as describing some elastic or viscous behaviour of the ambient space N that gets deformed as the shape changes. One might call these metrics outer metrics.

Contributions of this work.
• This work is the first to treat Sobolev inner metrics on spaces of immersed surfaces and on higher dimensional shape spaces. • It contains the first description of how the geodesic equation can be formulated in terms of gradients of the metric with respect to itself when the ambient space is not flat. To achieve this, a covariant derivative on some bundles over immersions is defined. This covariant derivative is induced from the Levi-Civita covariant derivative on ambient space. • The geodesic equation is formulated in terms of this covariant derivative.
Well-posedness of the geodesic equation is shown under some regularity assumptions that are verified for Sobolev metrics. Well-posedness also follows for the geodesic equation on diffeomorphism groups, where this result has not yet been obtained in that full generality. • To derive the geodesic equation, a variational formula for the Laplacian operator is developed. The variation is taken with respect to the metric on the manifold where the Laplacian is defined. This metric in turn depends on the immersion inducing it. • It is shown that Sobolev inner metrics separate points in shape space when the order of the differential operator governing the metric is high enough.
(The metric needs to be as least as strong as the H 1 -metric.) Thus Sobolev inner metrics overcome the degeneracy of the L 2 -metric. • The path-length distance of Sobolev inner metrics is compared to the Fréchet distance. It would be desirable to bound Féchet distance by some Sobolev distance. This however remains an open problem. • Finally it is demonstrated in some examples that the geodesic equation for the H 1 -metric on shape space of surfaces in R 3 can be solved numerically.
Big parts of this work can also be found, partly in more details, in the doctoral theses of Martin Bauer [4] and Philipp Harms [11].

Content of this work
This work progresses from a very general setting to a specific one in three steps. In the beginning, a framework for general inner metrics is developed. Then the general concepts carry over to more and more specific inner metrics.
• First, shape space is endowed with a general inner metric, i.e with a metric that is induced from a metric on the space of immersions, but that is unspecified otherwise. It is shown how various versions of the geodesic equation can be expressed using gradients of the metric with respect to itself and how conserved quantities arise from symmetries. (This is section 4.) • Then it is assumed that the inner metric is defined via an elliptic pseudodifferential operator. Such a metric will be called a Sobolev-type metric.
The geodesic equation is formulated in terms of the operator, and existence of horizontal paths of immersions within each equivalence class of paths is proven. (This is section 6.) Then estimates on the path-length distance are derived. Most importantly it is shown that when the operator involves high enough powers of the Laplacian, then the metric does not have the degeneracy of the L 2 -metric. (This is section 7.) • Motivated by the previous results it is assumed that the elliptic pseudodifferential operator is given by the Laplacian and powers of it. Again, the geodesic equation is derived. The formulas that are obtained are ready to be implemented numerically. (This is section 8.) The remaining sections cover the following material: • Section 3 treats some differential geometry of surfaces that is needed in this work. It is also a good reference for the notation that is used. The biggest emphasis is on a rigorous treatment of the covariant derivative.
Some material like the adjoint covariant derivative is not found in standard text books. • Section 5 contains formulas for the variation of the metric, volume form, covariant derivative and Laplacian with respect to the immersion inducing them. These formulas are used extensively later. • Section 9 covers the special case of flat ambient space. The geodesic equation is simplified and conserved momenta for the Euclidean motion group are calculated. Sobolev-type metrics are compared to the Fréchet metric which is available in flat ambient space. • Section 10 treats diffeomorphism groups of compact manifolds as a special case of the theory that has been developed so far. • In section 11 it is shown in some examples that the geodesic equation on shape space can be solved numerically.

Differential geometry of surfaces and notation
In this section the differential geometric tools that are needed to deal with immersed surfaces are presented and developed. The most important point is a rigorous treatment of the covariant derivative and related concepts.
The notation of [17] is used. Some of the definitions can also be found in [12]. A similar exposition in the same notation is [2].

Basic assumptions and conventions.
Assumption. It is always assumed that M and N are connected manifolds of finite dimensions m and n, respectively. Furthermore it is assumed that M is compact, and that N is endowed with a Riemannian metric g.
In this work, immersions of M into N will be treated, i.e. smooth functions M → N with injective tangent mapping at every point. The set of all such immersions will be denoted by Imm(M, N ). It is clear that only the case dim(M ) ≤ dim(N ) is of interest since otherwise Imm(M, N ) would be empty.
Immersions or paths of immersions are usually denoted by f . Vector fields on Imm(M, N ) or tangent vectors with foot point f , i.e., vector fields along f , will be called h, k, m, for example. Subscripts like f t = ∂ t f = ∂f /∂t denote differentiation with respect to the indicated variable, but subscripts are also used to indicate the foot point of a tensor field.

Tensor bundles and tensor fields. The tensor bundles
will be used. Here T r s M denotes the bundle of ( r s )-tensors on M , i.e.
and f * T N is the pullback of the bundle T N via f , see [17, section 17.5]. A tensor field is a section of a tensor bundle. Generally, when E is a bundle, the space of its sections will be denoted by Γ(E).
To clarify the notation that will be used later, some examples of tensor bundles and tensor fields are given now. S k T * M = L k sym (T M ; R) and Λ k T * M = L k alt (T M ; R) are the bundles of symmetric and alternating ( 0 k )-tensors, respectively. Ω k (M ) = Γ(Λ k T * M ) is the space of differential forms, X(M ) = Γ(T M ) is the space of vector fields, and is the space of vector fields along f .

3.3.
Metric on tensor spaces. Let g ∈ Γ(S 2 >0 T * N ) denote a fixed Riemannian metric on N . The metric induced on M by f ∈ Imm(M, N ) is the pullback metric where X, Y are vector fields on M . The dependence of g on the immersion f should be kept in mind. Let g can be extended to the cotangent bundle T * M = T 0 1 M by setting for α, β ∈ T * M . The product metric extends g to all tensor spaces T r s M , and g r s ⊗ g yields a metric on T r s M ⊗ f * T N .

3.4.
Traces. The trace contracts pairs of vectors and co-vectors in a tensor product: A special case of this is the operator i X inserting a vector X into a co-vector or into a covariant factor of a tensor product. The inverse of the metric g can be used to define a trace contracting pairs of co-vecors. Note that Tr g depends on the metric whereas Tr does not. The following lemma will be useful in many calculations: Lemma.
(In the expression under the trace, B and C are seen as maps T M → T * M .) Proof. Express everything in a local coordinate system u 1 , . . . , u m of M .
Note that only the symmetry of C has been used.
The volume of the immersion is given by The integral is well-defined since M is compact. If M is oriented the volume density may be identified with a differential form.
3.6. Metric on tensor fields. A metric on a space of tensor fields is defined by integrating the appropriate metric on the tensor space with respect to the volume density: for B, C ∈ Γ(T r s M ), and N ). The integrals are well-defined because M is compact. 3.7. Covariant derivative. Covariant derivatives on vector bundles as explained in [17, sections 19.12, 22.9] will be used. Let ∇ g , ∇ g be the Levi-Civita covariant derivatives on (M, g) and (N, g), respectively. For any manifold Q and vector field X on Q, one has Usually the symbol ∇ will be used for all covariant derivatives. It should be kept in mind that ∇ g depends on the metric g = f * g and therefore also on the immersion f . The following properties hold [17, section 22.9]: (1) ∇ X respects base points, i.e. π • ∇ X h = π • h, where π is the projection of the tangent space onto the base manifold.
(5) For any manifold Q and smooth mapping q : Q → Q and Y y ∈ T y Q one has ∇ T q.Yy h = ∇ Yy (h • q). If Y ∈ X(Q 1 ) and X ∈ X(Q) are q-related, then The two covariant derivatives ∇ g X and ∇ g X can be combined to yield a covariant derivative ∇ X acting on C ∞ (Q, T r s M ⊗ T N ) by additionally requiring the following properties [17, section 22.12]: , a derivation with respect to the tensor product. (8) ∇ X commutes with any kind of contraction (see [17, section 8.18]). A special case of this is Property (1) is important because it implies that ∇ X respects spaces of sections of bundles. For example, for Q = M and f ∈ C ∞ (M, N ), one gets 3.8. Swapping covariant derivatives. Some formulas allowing to swap covariant derivatives will be used repeatedly. Let f be an immersion, h a vector field along f and X, Y vector fields on M . Since ∇ is torsion-free, one has [17, section 22.10]: (1) Furthermore one has [17, section 24.5]: where R g ∈ Ω 2 N ; L(T N, T N ) is the Riemann curvature tensor of (N, g).
These formulas also hold when f : R × M → N is a path of immersions, h : R × M → T N is a vector field along f and the vector fields are vector fields on R × M . A case of special importance is when one of the vector fields is (∂ t , 0 M ) and the other (0 R , Y ), where Y is a vector field on M . Since the Lie bracket of these vector fields vanishes, (1) and (2) yield 3.9. Second and higher covariant derivatives. When the covariant derivative is seen as a mapping , then the second covariant derivative is simply ∇∇ = ∇ 2 . Since the covariant derivative commutes with contractions, ∇ 2 can be expressed as Higher covariant derivates are defined accordingly as ∇ k , k ≥ 0. g r s+1 (∇B, C) = g r s (B, ∇ * C). In the same way, ∇ * can be defined when ∇ is acting on Γ(T r s M ⊗ f * T N ). In either case it is given by , where the trace is contracting the first two tensor slots of ∇B. This formula will be proven now: Proof. The result holds for decomposable tensor fields β ⊗ B ∈ Γ(T r s+1 M ) since = g r s − div(β )B + Tr g (∇β)B − Tr g (∇(β ⊗ B)), C = g r s 0 − Tr g (∇(β ⊗ B)), C Here it has been used that ∇ X g = 0, that ∇ X commutes with any kind of contraction and acts as a derivation on tensor products [17, section 22.12] and that div(X) = Tr(∇X) for all vector fields X [17, section 25.12]. To prove the result for β ⊗ B ∈ Γ(T r s+1 M ⊗ f * T N ) one simply has to replace g r s by g r s ⊗ g.

3.11.
Laplacian. The definition of the Laplacian used in this work is the Bochner-Laplacian. It can act on all tensor fields B and is defined as 3.12. Normal bundle. The normal bundle Nor(f ) of an immersion f is a subbundle of f * T N whose fibers consist of all vectors that are orthogonal to the image of f : If dim(M ) = dim(N ) then the fibers of the normal bundle are but the zero vector. Any vector field h along f ∈ Imm can be decomposed uniquely into parts tangential and normal to f as h = T f.h + h ⊥ , where h is a vector field on M and h ⊥ is a section of the normal bundle Nor(f ).
3.13. Second fundamental form and Weingarten mapping. Let X and Y be vector fields on M . Then the covariant derivative ∇ X T f.Y splits into tangential and a normal parts as S is the second fundamental form of f . It is a symmetric bilinear form with values in the normal bundle of f . When T f is seen as a section of T * M ⊗ f * T N one has S = ∇T f since The trace of S is the vector valued mean curvature Tr g (S) ∈ Γ Nor(f ) .

Shape space
Briefly said, in this work the word shape means an unparametrized surface. (The term surface is used regardless of whether it has dimension two or not.) This section is about the infinite dimensional space of all shapes. First some spaces of parametrized and unparametrized surfaces are described, and it is shown how to define Riemannian metrics on them. The geodesic equation and conserved quantities arising from symmetries are derived.
The agenda that is set out in this section will be pursued in section 6 when the arbitrary metric is replaced by a Sobolev-type metric involving a pseudo-differential operator and later in section 8 when the pseudo-differential operator is replaced by an operator involving powers of the Laplacian. where denotes the c ∞ -completed bornological tensor product of locally convex vector spaces [14, section 5.7, section 4.29]. Note that L(T Imm; T Imm) is not isomorphic to T * Imm ⊗ T Imm since the latter bundle corresponds to multilinear mappings with finite rank.
It is worth to write down more explicitly what some of these bundles of multilinear mappings are. The tangent space to Imm is given by . Thus T f Imm is the space of vector fields along the immersion f . Now the cotangent space to Imm will be described. The symbol ⊗ C ∞ (M ) means that the tensor product is taken over the algebra C ∞ (M ).
is of interest for the definition of a Riemannian metric on Imm. (The subscripts sym and alt indicate symmetric and alternating multilinear maps, respectively.) Letting ⊗ S denotes the symmetric tensor product and ⊗ S the c ∞ -completed bornological symmetric tensor product, one has Imm is a section of the bundle L 2 sym (T Imm; R) such that at every f ∈ Imm, G f is a symmetric positive definite bilinear mapping Each metric is weak in the sense that G f , seen as a mapping f Imm is injective. (But it can never be surjective.) 4.2. Covariant derivative ∇ g on immersions. The covariant derivative ∇ g defined in section 3.7 induces a covariant derivative over immersions as follows. Let Q be a smooth manifold. Then one identifies h ∈ C ∞ Q, T Imm(M, N ) and As described in section 3.7 one has the covariant derivative This covariant derivative is torsion-free by section 3.8, formula (1). It respects the metric g but in general does not respect G.
It is helpful to point out some special cases of how this construction can be used. The case Q = R will be important to formulate the geodesic equation. The expression that will be of interest in the formulation of the geodesic equation is ∇ ∂t f t , which is well-defined when f : R → Imm is a path of immersions and f t : R → T Imm is its velocity.
Another case of interest is Q = Imm. Let h, k, m ∈ X(Imm). Then the covariant derivative ∇ m h is well-defined and tensorial in m. Requiring ∇ m to respect the grading of the spaces of multilinear maps, to act as a derivation on products and to commute with compositions of multilinear maps, one obtains as in section 3.7 a covariant derivative ∇ m acting on all mappings into the natural bundles of multilinear mappings over Imm. In particular, ∇ m P and ∇ m G are well-defined for

Metric gradients.
The metric gradients H, K ∈ Γ L 2 (T Imm; T Imm) are uniquely defined by the equation where h, k, m are vector fields on Imm and the covariant derivative of the metric tensor G is defined as in the previous section. (This is a generalization of the definition used in [20] that allows for a curved ambient space N = R n .) Existence of H, K has to proven case by case for each metric G, usually by partial integration. For Sobolev metrics, this will be proven in sections 8.2 and 8.3.
Assumption. Nevertheless it will be assumed for now that the metric gradients H, K exist.

Geodesic equation on immersions.
Theorem. Given H, K as defined in the previous section and ∇ as defined in section 4.2, the geodesic equation reads as This is the same result as in [20, section 2.4], but in a more general setting.
Proof. Let f : (−ε, ε) × [0, 1] × M → N be a one-parameter family of curves of immersions with fixed endpoints. The variational parameter will be denoted by s ∈ (−ε, ε) and the time-parameter by t ∈ [0, 1]. In the following calculation, let Remember that the covariant derivative on Imm that has been introduced in section 4.2 is torsion-free so that one has Thus the first variation of the energy of the curves is If f (0, ·, ·) is energy-minimizing, then one has at s = 0 that

Geodesic equation on immersions in terms of the momentum.
In the previous section the geodesic equation for the velocity f t has been derived. In many applications it is more convenient to formulate the geodesic equation as an equation is an element of the smooth cotangent bundle, also called smooth dual, which is given by It is strictly smaller than T * Imm since at every f ∈ Imm the metric G f : T f Imm → T * f Imm is injective but not surjective. It is called smooth since it does not contain distributional sections of f * T N , whereas T * f Imm does. Theorem. The geodesic equation for the momentum p ∈ T * Imm is given by where H is the metric gradient defined in section 4.3 and ∇ is the covariant derivative action on mappings into T * Imm as defined in section 4.2.
Proof. Let G f denote G composed with the path f : R → Imm, i.e.
This equation is equivalent to Hamilton's equation restricted to the smooth cotangent bundle: Here ω denotes the restriction of the canonical symplectic form on T * Imm to the smooth cotangent bundle and E is the Hamiltonian which is only defined on the smooth cotangent bundle.
4.6. Shape space. Diff(M ) acts smoothly on Imm(M, N ) and Emb(M, N ) by composition from the right. For Imm, the action is given by the mapping The tangent prolongation of this group action is given by the mapping Shape space is defined as the orbit space with respect to this action. That means that in shape space, two mappings differing only in their parametrization will be regarded the same.
Theorem. Let M be compact and of dimension ≤ n. Then Emb(M, N ) is the total space of a smooth principal fiber bundle with structure group Diff(M ), whose base manifold is a Hausdorff smooth Fréchet manifold denoted by The proof for immersions can be found in [6] and the one for embeddings in [14, section 44.1]. As with immersions and embeddings, the notation B i , B e will be used when it is clear that M and N are the domain and target of the mappings. 4.7. Riemannian metrics on shape space. We start with a metric G on Imm. The mapping π : Imm → B i is a submersion of smooth manifolds, that is, T π : Hor = Hor(π, G) := V (π) ⊥ ⊂ T Imm. It need not be a complement to V (recall that the metric is weak; the complement could be in a suitable completion of the tangent space). For all metrics in this paper it will turn out to be a complement, however. Then any vector h ∈ T Imm can be decomposed uniquely in vertical and horizontal components as This definition extends to the cotangent bundle as follows: An element of T * Imm is called horizontal when it annihilates all vertical vectors, and vertical when it annihilates all horizontal vectors.
In the setting described so far, the mapping is an isomorphism of vector spaces for all f ∈ Imm. This isomorphism will be used to describe the tangent space to B i . If both Imm and B i are Riemannian manifolds and if this isomorphism is also an isometry for all f ∈ Imm, then π is called a Riemannian submersion. In that case, the metric G on Imm is Diff(M )-invariant. This means that G = (r ϕ ) * G for all ϕ ∈ Diff(M ), where r ϕ denotes the right action of ϕ on Imm that was described in section 4.6. This condition can be spelled out in more details using the definition of r ϕ as follows: The following theorem establishes the converse statement: Theorem. Given a Diff(M )-invariant Riemannian metric on Imm, there is a unique Riemannian metric on the quotient space B i such that the quotient map π : Imm → B i is a Riemannian submersion.
Proof. If the horizontal bundle Hor f is a complement to V f then T f π : Hor f → T π(f ) B i is an isomorphism (off the orbifold singularities of B i ) and we can induce the metric on T π(f ) B i which is independent of the choice of f in the fiber over π(f ) by the the Diff(M )-invariance of the metric. If it is not a complement one has to consider the metric quotient norm. See for example [19, section 3].
Assumption. It will always be assumed that a Diff(M )-invariant metric G on Imm(M, N ) is given and that shape space B i is endowed with the unique metric such that the quotient map is a Riemannian submersion.

4.8.
Riemannian submersions and geodesics. It follows from the general theory of Riemannian submersions that horizontal geodesics in the top space correspond nicely to geodesics in the quotient space: Theorem. Let c : [0, 1] → Imm be a geodesic. Theorem. Assuming that every curve in B i can be lifted to a horizontal curve in Imm, the geodesic equation on shape space is equivalent to where f is a horizontal curve in Imm, where H, K are the metric gradients defined in section 4.3 and where ∇ is the covariant derivative defined in section 4.2.
This is a consequence of the Diff(M )-invariance of the metric G and the conservation of the reaparametrization momentum. A general proof can be found in [11, section 3.14].
It will be shown in section 6.9 that curves in B i can be lifted to horizontal curves in Imm for the very general class of Sobolev type metrics. Thus all assumptions and conclusions of the theorem hold. 4.10. Geodesic equation on shape space in terms of the momentum. As in the previous section, theorem 4.8 will be applied to the Riemannian submersion π : Imm → B i . But this time, the formulation of the geodesic equation in terms of the momentum will be used, see section 4.5. As will be seen in section 6.11, this is the most convenient formulation of the geodesic equation for Sobolev-type metrics.
Theorem. Assuming that every curve in B i can be lifted to a horizontal curve in Imm, the geodesic equation on shape space is equivalent to the set of equations Here f is a curve in Imm, H is the metric gradient defined in section 4.3, and ∇ is the covariant derivative defined in section 4.2. f is horizontal because p is horizontal.

Variational formulas
Recall that many operators like implicitly depend on the immersion f . In this section their derivative with respect to f which is called their first variation will be calculated . These formulas will be used to calculate the metric gradients that are needed for the geodesic equation.
This section is based on [2], see also [11]. Some but not all of the formulas were known before [5,18]. More variational formulas can be found in [5,23,4]. 5.1. Paths of immersions. All of the differential-geometric concepts introduced in section 3 can be recast for a path of immersions instead of a fixed immersion. This allows to study variations of immersions. So let f : R → Imm(M, N ) be a path of immersions. By convenient calculus [14], f can equivalently be seen as is an immersion for each t. The bundles over M can be replaced by bundles over R × M : The covariant derivative ∇ Z h is now defined for vector fields Z on R × M and sections h of the above bundles. The vector fields (∂ t , 0 M ) and (0 R , X), where X is a vector field on M , are of special importance. In later sections they will be identified with ∂ t and X whenever this does not pose any problems. Let Then by property 5 from section 3.7 one has for vector fields X, Y on M This shows that one can recover the static situation at t by using vector fields on R × M with vanishing R-component and evaluating at t.

5.2.
Directional derivatives of functions. The following ways to denote directional derivatives of functions will be used, in particular in infinite dimensions. Given a function F (x, y) for instance, D (x,h) F will be written as a shorthand for ∂ t | 0 F (x + th, y).
Here (x, h) in the subscript denotes the tangent vector with foot point x and direction h. If F takes values in some linear space, this linear space and its tangent space will be identified.

Setting for first variations.
In all of this chapter, let f be an immersion and f t ∈ T f Imm a tangent vector to f . The reason for calling the tangent vector f t is that in calculations it will often be the derivative of a curve of immersions through f . Using the same symbol f for the fixed immersion and for the path of immersions through it, one has in fact that

Variation of equivariant tensor fields. Let the mapping
F : Imm(M, N ) → Γ(T r s M ) take values in some space of tensor fields over M , or more generally in any natural bundle over M , see [13].
Lemma. If F is equivariant with respect to pullbacks by diffeomorphisms of M , i.e.
for all ϕ ∈ Diff(M ) and f ∈ Imm(M, N ), then the tangential variation of F is its Lie-derivative: This allows us to calculate the tangential variation of the pullback metric and the volume density, for example.

Variation of the metric.
Lemma. The differential of the pullback metric Here Sym denotes the symmetric part of the tensor field C of type ( 0 2 ) given by Proof. Let f : R×M → N be a path of immersions. Swapping covariant derivatives as in section 3.8, formula (3) one gets Splitting f t into its normal and tangential part yields follows either from the equivariance of g with respect to pullbacks by diffeomorphisms (see section 5.4) or directly from

Variation of the inverse of the metric.
Lemma. The differential of the inverse of the pullback metric is given by Proof.
Variation of the volume density.
Lemma. The differential of the volume density is given by Proof. Let g(t) ∈ Γ(S 2 >0 T * M ) be any curve of Riemannian metrics. Then This follows from the formula for vol(g) in a local oriented chart (u 1 , . . . u m ) on M : Now one can set g = f * g and plug in the formula This immediately proves the first formula: Expanding this further yields the second formula: Here it has been used that Note that by 5.4, the formula for the tangential variation would have followed also from the equivariance of the volume form with respect to pullbacks by diffeomorphisms.

5.8.
Variation of the covariant derivative. In this section, let ∇ = ∇ g = ∇ f * g be the Levi-Civita covariant derivative acting on vector fields on M . Since any two covariant derivatives on M differ by a tensor field, the first variation of ∇ f * g is tensorial. It is given by the tensor field D (f,ft) ∇ f * g ∈ Γ(T 1 2 M ).
Lemma. The tensor field D (f,ft) ∇ f * g is determined by the following relation holding for vector fields X, Y, Z on M : Proof. The defining formula for the covariant derivative is Taking the derivative D (f,ft) yields Then the result follows by replacing all Lie brackets in the above formula by covariant derivatives using [X, Y ] = ∇ X Y − ∇ Y X and by expanding all terms of the form X (D (f,ft g)(Y, Z) using

5.9.
Variation of the Laplacian. The Laplacian as defined in section 3.11 can be seen as a smooth section of the bundle L(T Imm; T Imm) over Imm since for every f ∈ Imm it is a mapping The right way to define a first variation is to use the covariant derivative defined in section 4.2.
Lemma. For ∆ ∈ Γ L(T Imm; T Imm) , f ∈ Imm and f t , h ∈ T f Imm one has Proof. Let f be a curve of immersions and h a vector field along f . One has ∆ : Imm → L(T Imm; T Imm), ∆ • f = ∆ f * g : R → Imm → L(T Imm; T Imm).
Using property 3.7.5 one gets The term Tr g (∇ ∂t ∇ 2 h) will be treated further. Let X, Y be vector fields on M that are constant in time. When they are seen as vector fields on R × M then ∇ ∂t X = ∇ ∂t Y = 0. Using the formulas from section 3.8 to swap covariant derivatives one gets The Lie bracket is since (now without the slight abuse of notation) Putting together all terms one obtains It remains to calculate Tr g (D (f,ft) ∇). Using the variational formula for ∇ from section 5.8 one gets for any vector field Z and a g-orthonormal frame s i Therefore Tr g (D (f,ft) ∇) = − ∇ * (D (f,ft) g) + 1 2 d Tr g (D (f,ft) g) .

Sobolev-type metrics
Assumption. Let P be a smooth section of the bundle L(T Imm; T Imm) over Imm such that at every f ∈ Imm the operator is an elliptic pseudo differential operator that is symmetric and positive with respect to the H 0 -metric on Imm, Note that an elliptic symmetric operator is self-adjoint by [21, 26.2]. Then P induces a metric on the set of immersions, namely The metric G P is positive definite since P is assumed to be positive with respect to the H 0 -metric. In this section, the geodesic equation on Imm and B i for the G P -metric will be calculated in terms of the operator P and it will be proven that it is well-posed under some assumptions.
Assumption. It will be assumed that P is invariant under the action of the reparametrization group Diff(M ) acting on Imm (M, N ), i.e. P = (r ϕ ) * P for all ϕ ∈ Diff(M ).
For any f ∈ Imm and ϕ ∈ Diff(M ) this means Applied to h ∈ T f Imm this means The invariance of P implies that the induced metric G P is invariant under the action of Diff(M ), too. Therefore it induces a unique metric on B i as explained in section 4.7 6.2. The adjoint of ∇P . The following construction is needed to express the metric gradient H which is part of the geodesic equation. H f arises from the metric G f by differentiating it with respect to its foot point f ∈ Imm. Since G is defined via the operator P , one also needs to differentiate P f with respect to its foot point. As for the metric, this is accomplished by the covariant derivate. For P ∈ Γ L(T Imm; T Imm) and m ∈ T Imm one has ∇ m P ∈ Γ L(T Imm; T Imm) , ∇P ∈ Γ L(T 2 Imm; T Imm) .

See section 4.2 for more details.
Assumption. It is assumed that there exists a smooth adjoint Adj(∇P ) ∈ Γ L 2 (T Imm; T Imm) of ∇P in the following sense: The existence of the adjoint needs to be checked in each specific example, usually by partial integration. For the operator P = 1 + A∆ p , the existence of the adjoint will be proven and explicit formulas will be calculated in sections 8.2 and 8.3.
Lemma. If the adjoint of ∇P exists, then its tangential part is determined by the invariance of P with respect to reparametrizations: Adj(∇P )(h, k) = g(∇P h, k) − g(∇h, P k) = grad g g(P h, k) − g(P h, ∇k) + g(∇h, P k) for f ∈ Imm, h, k ∈ T f Imm.
Proof. Let X be a vector field on M . Then g m , g(∇P h, k) − g(∇h, P k) vol(g). 6.3. Metric gradients. As explained in section 4.4, the geodesic equation can be expressed in terms of the metric gradients H and K. These gradients will be computed now. We shall use that P f is invertible on the space of smooth sections. This follows because P f is an elliptic, self-adjoint, and positive operator, see the beginning of the proof of theorem 6.6 for a detailed argument.
Lemma. If Adj(∇P ) exists, then also H and K exist and are given by Proof. For vector fields m, h, k on Imm one has (1) One immediately gets the K-gradient by plugging in the variational formula 5.7 for the volume form: To calculate the H-gradient, one rewrites equation (1)  Thus the H-gradient is given by The highest order term grad g g(P h, k) cancels out when taking into account the formula for the tangential part of the adjoint from section 6.2: − g(P h, k). Tr g (S) .
Plugging in the formulas for H, K derived in the last section yields the following theorem.
Theorem. The geodesic equation for a Sobolev-type metric G P on immersions is given by Theorem. The geodesic equation written in terms of the momentum for a Sobolevtype metric G P on Imm is given by: 6.6. Well-posedness of the geodesic equation. It will be proven that the geodesic equation for a Sobolev-type metric G P on Imm is well-posed under some assumptions on P . These assumptions are satisfied for the operator 1 + A∆ p considered in section 8. It will also be shown that (π, exp) is a diffeomorphism from a neighbourhood of the zero section in T Imm to a neighbourhood of the diagonal in Imm × Imm.
Before we can state the theorem, we have to introduce Sobolev completions of the relevant spaces of mappings. More information can be found in [21], [9], and in [8]. We consider Sobolev completions of Γ(E), where E → M is a vector bundle. First we choose a fixed (background) Riemannian metricĝ on M and its covariant derivative ∇ M . We equip E with a (background) fiber Riemannian metricĝ E and a compatible covariant derivative∇ E . Then the Sobolev space H k (E) is the Hilbert space completion of the space of smooth sections Γ(E) in the Sobolev norm This Sobolev space does not depend on the choices ofĝ, ∇ M ,ĝ E and∇ E since M is compact: The resulting norms are equivalent.
We shall need the following results (see [8], e.g.): Sobolev lemma. L 2 (T Imm; T Imm) Imm Imm Imm, respectively. Viewed locally in trivializations of these bundles, ⊥ are pseudo-differential operators of order 2p in h, k separately. As mappings in the footpoint f they are non-linear, and it is assumed that they are a composition of operators of the following type: (a) Local operators of order l ≤ 2p, i.e., nonlinear differential operators (b) Linear pseudo-differential operators of degrees l i , such that the total (top) order of the composition is ≤ 2p.

Assumption 2.
For each f ∈ Imm(M, N ), the operator P f is an elliptic pseudodifferential operator of order 2p for p > 0 which is positive and symmetric with respect to the H 0 -metric on Imm, i.e. Moreover, in each Sobolev completion Imm k+2p , the Riemannian exponential mapping exp P exists and is smooth on a neighborhood of the zero section in the tangent bundle, and (π, exp P ) is a diffeomorphism from a (smaller) neigbourhood of the zero section to a neighborhood of the diagonal in Imm k+2p × Imm k+2p . All these neighborhoods are uniform in k > dim(M )/2 + 1 and can be chosen Proof. By assumption 1 the mapping P f h is of order 2p in f and in h where f is the footpoint of h. Therefore f → P f extends to a smooth section of the smooth Sobolev bundle where T Imm k | Imm k+2p denotes the space of all H k tangent vectors with foot point a H k+2p immersion, i.e., the restriction of the bundle T Imm k → Imm k to Imm k+2p ⊂ Imm k .
This means that P f is a bounded linear operator It is injective since it is positive. As an elliptic operator, it is an unbounded operator on the Hilbert completion of T f Imm with respect to the H 0 -metric, and a Fredholm operator H k+2p → H k for each k. It is selfadjoint elliptic, thus by [21, theorem 26.2] it has vanishing index. Since it is injective, it is thus also surjective.
By the implicit function theorem on Banach spaces, f → P −1 f is then a smooth section of the smooth Sobolev bundle L T Imm k | Imm k+2p ; T Imm k+2p → Imm k+2p As an inverse of an elliptic pseudodifferential operator, P −1 f is also an elliptic pseudo-differential operator of order −2p.
Using the module property of Sobolev spaces and counting the order of all remaining terms in the geodesic equation 6.4, one obtains that the Christoffel symbols extend to a smooth (C ∞ ) section of the smooth Sobolev bundle is a smooth quadratic mapping T Imm → T Imm which extends to smooth quadratic mappings T Imm k+2p → T Imm k+2p for each k ≥ dim(2) 2 + 1. The geodesic equation can be reformulated using the linear connection C g : T N × N T N → T T N (horizontal lift mapping) of ∇ g , see [17, section 24.2]: The right-hand side is a smooth vector field on T Imm k+2p , the geodesic spray. Note that the restriction to T Imm k+1+2p of the geodesic spray on T Imm k+2p equals the geodesic spray there. By the theory of smooth ODE's on Banach spaces, the flow of this vector field exists in T Imm k+2p and is smooth in time and in the initial condition.
Consider a C ∞ initial condition h 0 ∈ T Imm with foot point f 0 ∈ Imm. Suppose the trajectory Fl k t (h 0 ) of geodesic spray through these initial conditions in T Imm k+2p maximally exists for t ∈ (−a k , b k ), and the trajectory Fl k+1 t (h 0 ) in T Imm k+1+2p maximally exists for t ∈ (−a k+1 , b k+1 ) with b k+1 < b k , say. By uniqueness of solutions one has Fl k+1 t (h 0 ) = Fl k t (h 0 ) for t ∈ (−a k+1, b k+1 ). We now write∇ for the covariant derivative induced by g on N and the background metriĉ g on M . Let X be a vector field on M . Applying ∇ g X =:∇ X to the geodesic equation and swapping covariant derivatives yields: Thus we can omit X and rewrite (A) as an equation for T f . We aim to rewrite equation (A) as a linear first order equation for the highest derivative whose coefficients are given by Fl k t (h 0 ) and thus exist beyond (−a k+1 , b k+1 ). For this we have to pass to one (ore more) canonical chart for Imm and the induced trivializations of all bundles as before and in assumption (1). Then f itself has values in a vector space and we may regard (A) as a vector valued 1-form on M . So we rewrite (A) as: We claim that (B) consists of: in f, f t ; order here means that the expression prolongs continuously to the corresponding Sobolev spaces.
To see this we claim that the highest derivatives of order 2p + 1 of f and f t appear only linearly in (A). This claim follows from assumption 1: (a) For a local operator we can apply the chain rule: The highest derivative of f appears only linearly. (b) For a linear pseudo differential operator A of order k the commutator [∇, A] is a pseudo-differential operator of order k again. On the left hand side of (B) we write T f =: u and∇f t =: v. On the right hand side of (B) we write∇ 2p T f =∇ 2p u and∇ 2p+1 f t =∇ 2p v for the highest derivatives only. Then the system (B) becomes: The coefficients f, f t in (C) exist for t ∈ (−a k , b k ) as Fl k t (h 0 ). Then (C) is a bounded and smooth inhomogeneous linear ODE for (u, v) ∈ Ω 1 (M, T Imm k+2p ), i.e., in a Banach space. This equation therefore has a solution (u(t), v(t)) for all t for which the coefficients exists, thus for all t ∈ (a k , b k ), which is unique for the initial values u 0 = T f 0 and v 0 =∇h 0 . The limit lim t b k+1 (u(t), v(t)) exists in Ω 1 (M, T Imm k+2p ) and by continuity it equals (∇T f,∇f t ) for t = b k+1 . Thus the flow line Fl k t (h 0 ) was not maximal and can be continued. So assuming b k+1 < b k leads to a contradiction, and thus (−a k+1 , b k+1 ) = (−a k , b k ). Iterating this procedure one concludes that the flow line Fl ∞ t (h 0 ) exists in k≥ dim(M ) 2 +1 T Imm k+2p = T Imm.
It remains to check the properties of the Riemannian exponential mapping exp P . It is given by exp P f (h) = c(1) where c(t) is the geodesic emanating from value f with initial velocity h. Let k 0 > dim(M )/2 + 1 and k ≥ k 0 . On each space T Imm k+2p , the properties claimed follow from local existence and uniqueness of solutions to the flow equation of the geodesic spray, from the form of the geodesic equation when it is written down in a chart, namely linearity in f tt and bilinearity in f t , and from the inverse function theorem which holds on each of the Sobolev spaces Imm k+2p . See for example [17, 22.6 and 22.7,] for a detailed proof which works without any change in notation.
Imm k0+2p contains Imm k+2p for k > k 0 . Since the spray on Imm k0+2p restricts to the spray on each Imm k+2p , the exponential mapping exp P and the inverse (π, exp P ) −1 on Imm k0+2p restrict to the corresponding mappings on each Imm k+2p . Thus the neighborhoods of existence are uniform in k and can be chosen is constant along f .
Note that the horizontal bundle consists of vector fields that are normal to f when P = Id, i.e. for the H 0 -metric on Imm.
Let us work out the G P -decomposition of h into vertical and horizontal parts. This decomposition is written as Thus one considers the operators The operator P f is unbounded, positive and symmetric on the Hilbert completion of T f Imm with respect to the H 0 -metric since one has Let σ P f and σ P f denote the principal symbols of P f and P f , respectively. Take any x ∈ M and ξ ∈ T * x M \ {0}. Then σ P f (ξ) is symmetric, positive definite on (T f (x) N, g). This means that one has for any h, k ∈ T f (x) N that The principal symbols σ P f and σ P f are related by . Therefore P f is again elliptic, thus it is selfadjoint, so its index (as operator H k+2p → H k ) vanishes. It is injective (since positive) with vanishing index (since self-adjoint elliptic, by [21, theorem 26.2]) hence it is bijective and thus invertible by the open mapping theorem. Thus it has been proven: Lemma. The decomposition of h ∈ T Imm into its vertical and horizontal components is given by 6.9. Horizontal curves. To establish the one-to-one correspondence between horizontal curves in Imm and curves in shape space that has been described in theorem 4.8, one needs the following property: Lemma. For any smooth path f in Imm(M, N ) there exists a smooth path ϕ in Diff(M ) with ϕ(0, . ) = Id M depending smoothly on f such that the pathf given byf (t, x) = f (t, ϕ(t, x)) is horizontal: Thus any path in shape space can be lifted to a horizontal path of immersions.
The basic idea is to write the path ϕ as the integral curve of a time dependent vector field. This method is called the Moser-trick (see [18,Section 2.5]).
Proof. Since P is invariant, one has (r ϕ ) * P = P or P f •ϕ (u • ϕ) = (P f u) • ϕ for ϕ ∈ Diff(M ). In the following f • ϕ will denote the map f (t, ϕ(t, x)), etc. One looks for ϕ as the integral curve of a time dependent vector field ξ(t, x) on M , given by ϕ t = ξ • ϕ. The following expression must vanish for all x ∈ M and X x ∈ T x M : Since T ϕ is surjective, T ϕ.X exhausts the tangent space T ϕ(x) M , and one has This holds for all x ∈ M , and by the surjectivity of ϕ, one also has that at all x ∈ M . This means that the tangential part P f (∂ t f ) + P f (T f.ξ) vanishes.
Using the time dependent vector field and its flow ϕ achieves this.
6.10. Geodesic equation on shape space. By the previous section and theorem 4.8, geodesics in B i correspond exactly to horizontal geodesics in Imm. The equations for horizontal geodesics in the space of immersions have been written down in section 4.9. Here they are specialized to Sobolev-type metrics: Theorem. The geodesic equation on shape space for a Sobolev-type metric G P is equivalent to the set of equations where f is a horizontal path of immersions.
These equations are not handable very well since taking the horizontal part of a vector to Imm involves inverting an elliptic pseudo-differential operator, see section 6.8. However, the formulation in the next section is much better. 6.11. Geodesic equation on shape space in terms of the momentum. The geodesic equation in terms of the momentum has been derived in section 4.10 for a general metric on shape space. Now it is specialized to Sobolev-type metrics using the formula for the H-gradient from section 6.3.
As in section 6.5 the momentum G P (f t , ·) is identified with P f t ⊗ vol(g). By definition, the momentum is horizontal if it annihilates all vertical vectors. This is the case if and only if P f t is normal to f . Thus the splitting of the momentum in horizontal and vertical parts is given by This is much simpler than the splitting of the velocity in horizontal and vertical parts where a pseudo-differential operator has to be inverted, see section 6.8. Thus the following version of the geodesic equation on shape space is the easiest to solve.
Theorem. The geodesic equation on shape space is equivalent to the set of equations for a path of immersions f : The equation for geodesics on Imm without the horizontality condition is Tr g (S) ⊗ vol(g), see section 6.5. It has been proven in section 4.9 that the vertical part of this equation is satisfied automatically when the geodesic is horizontal. Nevertheless this will be checked by hand because the proof is much simpler here than in the general case.
If f t is horizontal then by definition P f t is normal to f . Thus one has for any X ∈ X(M ) that which is exactly the vertical part of the geodesic equation.

Geodesic distance on shape space
It came as a big surprise when it was discovered in [19] that the Sobolev metric of order zero induces vanishing geodesic distance on shape space B i . It will be shown that this problem can be overcome by using higher order Sobolev metrics. The proof of this result is based on bounding the G P -length of a path from below by its area swept out. The main result is in section 7. 7.1. Geodesic distance on shape space. Geodesic distance on B i is given by where the infimum is taken over all F : [0, 1] → B i with F (0) = F 0 and F (1) = F 1 . L Bi G P is the length of paths in B i given by Letting π : Imm → B i denote the projection, one has when f : [0, 1] → Imm is horizontal. In the following sections, conditions on the metric G P ensuring that dist Bi G P separates points in B i will be developed.
Theorem. The distance dist Bi H 0 induced by the Sobolev L 2 metric of order zero vanishes. Indeed it is possible to connect any two distinct shapes by a path of arbitrarily short length.
This result was first established by Michor and Mumford for the case of planar curves in [19]. A more general version can be found in [18], where the same result is proven also on diffeomorphism groups.  Lemma. Let G P be a Sobolev type metric that is at least as strong as the H 0metric, i.e. there is a constant C 1 > 0 such that Then one has the area swept out bound for any path of immersions f : The proof is an adaptation of the one given in [2, section 7.3] for almost local metrics. Proof.
Lemma. Let G P be a Sobolev type metric that is at least as strong as the H 1metric, i.e. there is a constant C 2 > 0 such that Then the mapping √ Vol : (B i , dist Bi G P ) → R ≥0 is Lipschitz continuous, i.e. for all F 0 and F 1 in B i one has: For the case of planar curves, this has been proven in [20, section 4.7]. Proof.
7.6. Non-vanishing geodesic distance. Using the estimates proven above and the fact that the area swept out separates points at least on B e , one gets the following result: Theorem. The Sobolev type metric G P induces non-vanishing geodesic distance on B e if it is stronger or as strong as the H 1 -metric, i.e. if there is a constant C > 0 such that Proof. By lemma 7.4 we have Now we use the Lipschitz continuity 7.5 of √ Vol and that area swept out separates points on B e .

Sobolev metrics induced by the Laplace operator
The results on non-vanishing geodesic distance from the previous section lead us to consider operators P that are induced by the Laplacian operator: for a constant A > 0. (See section 3.11 for the definition of the Laplacian that is used in this work.) At every f ∈ Imm, P f is a positive, selfadjoint and bijective operator of order 2p acting on T f Imm = Γ(f * T N ). Note that ∆ depends smoothly on the immersion f via the pullback-metric f * g, so that the same is true of P . P is invariant under the action of the reparametrization group Diff(M ). It induces the Sobolev metric When A = 1 we write H p := G 1+∆ p .
In this section we will calculate explicitly for P = 1 + A∆ p the geodesic equation and conserved momenta that have been deduced in section 6 for a general operator P . The hardest part will be the partial integration needed for the adjoint of ∇P . As a result we will get explicit formulas that are ready to be implemented numerically.
8.1. Other choices for P . Other choices for P are the operator P = 1+A(∇ * ) p ∇ p corresponding to the metric and other operators that differ only in lower order terms. Since these operators all have the same principal symbol, they induce equivalent metrics on each tangent space T f Imm. It would be interesting to know if the induced geodesic distances on B i are equivalent as well.
8.2. Adjoint of ∇P . To find a formula for the geodesic equation one has to calculate the adjoint of ∇P , see section 6.4. The following calculations at the same time show the existence of the adjoint. It has been shown in section 6.2 that the invariance of the operator P with respect to reparametrizations determines the tangential part of the adjoint: Adj(∇P )(h, k) = grad g g(P h, k) − g(P h, ∇k) + g(∇h, P k) .
It remains to calculate its normal part using the variational formulas from section 5.
In the following calculations there will be terms of the form Tr(g −1 s 1 g −1 s 2 ), where s 1 , s 2 are two-forms on M . When the two-forms are seen as mappings T M → T * M , they can be composed with g −1 : T * M → T M . Thus the expression under the trace is a mapping T M → T M to which the trace can be applied. When one of the two-forms is vector valued, the same tensor components as before are contracted. For example when h ∈ Γ(f * T N ) then s 2 = ∇ 2 h is a two-form on M with values in f * T N . Then in the expression Tr(g −1 .s 1 .g −1 .s 2 ) only T M and T * M components are contracted, whereas the f * T N component remains unaffected.
Using the following symmetry property of the curvature tensor (see [17, 24.4.4]): From this, one can read off the normal part of the adjoint. Thus one gets: Lemma. The adjoint of ∇P defined in section 6.2 for the operator P = 1 + A∆ p is 8.3. Geodesic equations and conserved momentum. The shortest and most convenient formulation of the geodesic equation is in terms of the momentum p = (1 + A∆ p )f t ⊗ vol(g), see sections 6.5 and 6.11.
Theorem. The geodesic equation on Imm(M, N ) for the G P -metric with P = 1 + A∆ p is given by: is constant along any geodesic f in Imm(M, N ).
The horizontal geodesic equation for a general metric on Imm has been derived in section 4.10. In section 6.11 it has been shown that this equation takes a very simple form. Now it is possible to write down this equation specifically for the operator P = 1 + A∆ p : Theorem. The geodesic equation on shape space for the Sobolev-metric G P with P = 1 + A∆ p is equivalent to the set of equations where f is a path of immersions. For the special case of plane curves, this agrees with the geodesic equation calculated in [20, section 4.6].

Surfaces in n-space
This section is about the special case where the ambient space N is R n . The flatness of R n leads to a simplification of the geodesic equation, and the Euclidean motion group acting on R n induces additional conserved quantities. The vector space structure of R n allows to define a Fréchet metric. This metric will be compared to Sobolev metrics. Finally in section 9.5 the space of concentric hyper-spheres in R n is briefly investigated. 9.1. Geodesic equation. The covariant derivative ∇ g on R n is but the usual derivative. Therefore the covariant derivatives ∇ ∂t f t and ∇ ∂t p in the geodesic equation can be replaced by f tt and p t , respectively. (Note that Imm(M, R n ) is an open subset of the Fréchet vector space C ∞ (M, R n ).) Also, the curvature terms disappear because R n is flat. Any of the formulations of the geodesic equation presented so far can thus be adapted to the case N = R n .
We want to show how the geodesic equation simplifies further under the additional assumptions that dim(M ) = dim(N ) − 1 and that M is orientable. Then it is possible define a unit vector field ν to M . The condition that f t is horizontal then simplifies to P f t = a.ν for a ∈ C ∞ (M ). The geodesic equation can then be written as an equation for a. However, the equation is slightly simpler when it is written as an equation for a.vol(g). In practise, vol(g) can be treated as a function on M because one can identify vol(g) with its density with respect to du 1 ∧ . . . ∧ du n−1 , where (u 1 , . . . , u n−1 ) is a chart on M . Thus multiplication by vol(g) does not pose a problem.
Theorem. The geodesic equation for a Sobolev-type metric G P on shape space B i (M, R n ) with dim(M ) = n − 1 is equivalent to the set of equations where f is a path in Imm(M, R n ) and a is a time-dependent function on M .
Proof. Applying g(·, ν) to the geodesic equation 6.11 on shape space in terms of the momentum one gets Let us spell this equation out in even more details for the H 1 -metric. This is the case of interest for the numerical examples in section 11.
Theorem. The geodesic equation on shape space B i (M, R n ) for the Sobolev-metric G P with P = 1 + A∆ is equivalent to the set of equations where f is a path of immersions, a is a time-dependent function on M , s = g(S, ν) ∈ Γ(T 0 2 M ) is the shape operator, L = g −1 s ∈ Γ(T 1 1 M ) is the Weingarten mapping, and Tr(L) is the mean curvature.
Proof. The fastest way to get to this equation is to apply g(·, ν) to the geodesic equation on Imm from section 8.3. This yields Notice that the second order derivatives of f t have canceled out. 9.2. Additional conserved momenta. If P is invariant under the action of the Euclidean motion group R n SO(n), then also the metric G P is in invariant under this group action and one gets additional conserved quantities as described in [20, section 2.5]: Theorem. For an operator P that is invariant under the action of the Euclidean motion group R n SO(n), the linear momentum It is claimed in [15, theorem 13] that d ∞ = dist G ∞ . However, the proof given there only works on the vector space C ∞ (M, R n ) and not on B i (M, R n ). The reason is that convex combinations of immersions are used in the proof, but that the space of immersions is not convex. 9.4. Sobolev versus Fréchet distance. It is a desirable property of any distance on shape space to be stronger than the Fréchet distance. Otherwise, singular points of a shape could move arbitrarily far away without increasing the distance much.
As the following result shows, Sobolev metrics of low order do not have this property. The authors believe that this property is true when the order of the metric is high enough, but were not able to prove this.
Lemma. Let G P be a metric on B i (M, R n ) that is weaker than or at least as weak as a Sobolev H p -metric with p < dim(M ) for all h ∈ T Imm.
Then the Fréchet distance can not be bounded by the G P -distance.
Proof. It is sufficient to prove the claim for P = 1+∆ p . Let f 0 be a fixed immersion of M into R n , and let f 1 be a translation of f 0 by a vector h of length . It will be shown that the H p -distance between π(f 0 ) and π(f 1 ) is bounded by a constant 2L that does not depend on , where π denotes the projection of Imm onto B i . Then it follows that the H p -distance can not be bounded from below by the Fréchet distance, and this proves the claim.
For small r 0 , one calculates the H p -length of the following path of immersions: First scale f 0 by a factor r 0 , then translate it by h, and then scale it again until it has reached f 1 . The following calculation shows that under the assumption p < m/2 + 1 the immersion f 0 can be scaled down to zero in finite H p -pathlength L. Let r : [0, 1] → [0, 1] be a function of time with r(0) = 1 and r(1) = 0.
The last integral converges if m−2p 2 < −1, which holds by assumption. Scaling down to r 0 > 0 needs even less effort. So one sees that the length of the shrinking and growing part of the path is bounded by 2L.
The length of the translation is simply r m 0 Vol(f 0 ) = O(r m/2 ) since the Laplacian of the constant vector field vanishes. Therefore 9.5. Concentric spheres. For a Sobolev type metric G P that is invariant under the action of the SO(n) on R n , the set of hyper-spheres in R n with common center 0 is a totally geodesic subspace of B i (S n−1 , R n ). The reason is that it is the fixed point set of the group SO(n) acting on B i isometrically. (One also needs uniqueness of solutions to the geodesic equation to prove that the concentric spheres are totally geodesic.) This section mainly deals with the case P = 1 + ∆ p .
First we want to determine under what conditions the set of concentric spheres is geodesically complete under the G P -metric.
Lemma. The space of concentric spheres is complete with respect to the G P metric with P = 1 + A∆ p iff p ≥ (n + 1)/2.
Proof. The space is complete if and only if it is impossible to scale a sphere down to zero or up to infinity in finite G P path-length. So let f be a path of concentric spheres. It is uniquely described by its radius r. Its velocity is f t = r t .ν, where ν designates the unit normal vector field. One has g g −1 .S, ν =: L = − 1 r Id T M , Tr(L k ) = (−1) k n−1 r k , Vol = r n−1 nπ n/2 Γ(n/2+1) .
Keep in mind that r and r t are constant functions on the sphere, so that all derivatives of them vanish. Therefore From this it is clear that the path f is horizontal. Therefore its length as a path in B i is the same as its length as a path in Imm. One calculates its length as in the proof of 9.4: n.π n/2 Γ(n/2 + 1) r n−1 dt = n.π n/2 Γ(n/2 + 1) r1 r0 1 + A (n − 1) p r 2p r n−1 dr.
The geodesic equation within the space of concentric spheres reduces to an ODE for the radius that can be read off the geodesic equation in section 8.3: A.(n − 1) p r r 2p + A(n − 1) p .

Diffeomorphism groups
For M = N the space Emb(M, M ) equals the diffeomorphism group of M . An operator P ∈ Γ L(T Emb; T Emb) that is invariant under reparametrizations induces a right-invariant Riemannian metric on this space. Thus one gets the geodesic equation for right-invariant Sobolev metrics on diffeomorphism groups and well-posedness of this equation. To the authors knowledge, well-posedness has so far only been shown for the special case M = S 1 in [7] and for the special case of Sobolev order one metrics in [10]. Theorem 6.6 establishes this result for arbitrary compact M and Sobolev metrics of arbitrary order.
In the decomposition of a vector h ∈ T f Emb into its tangential and normal components h = T f.h +h ⊥ , the normal part h ⊥ vanishes. Also S = ∇T f vanishes. Thus the geodesic equation on Diff(M ) in terms of the momentum p is given by (see 6.5) p = P f t ⊗ vol(g), Note that this equation is not right-trivialized, in contrast to the equation given in [1,18,16], for example. The special case of theorem 6.6 now reads as follows: Theorem. Let p ≥ 1 and k > dim(M ) 2 + 1 and let P satisfy assumptions (1-3) of 6.6.
Then the initial value problem for the geodesic equation has unique local solutions in the Sobolev manifold Diff k+2p of H k+2p -diffeomorphisms. The solutions depend smoothly on t and on the initial conditions f (0, . ) and f t (0, . ). The domain of existence (in t) is uniform in k and thus this also holds in Diff(M ).
Moreover, in each Sobolev completion Diff k+2p , the Riemannian exponential mapping exp P exists and is smooth on a neighborhood of the zero section in the tangent bundle, and (π, exp P ) is a diffeomorphism from a (smaller) neigbourhood of the zero section to a neighborhood of the diagonal in Diff k+2p × Diff k+2p . All these neighborhoods are uniform in k > dim(M )/2 + 1 and can be chosen H k0+2popen, for k 0 > dim(M )/2 + 1. Thus both properties of the exponential mapping continue to hold in Diff(M ).

Numerical results
It is of great interest for shape comparison to solve the boundary value problem for geodesics in shape space. When the boundary value problem can be solved, then any shape can be encoded as the initial momentum of a geodesic starting at a fixed reference shape. Since the initial momenta all lie in the same vector space, this also opens the way to statistics on shape space.
There are two approaches to solving the boundary value problem. In [2] the first approach of minimizing horizontal path energy over the set of curves in Imm connecting two fixed boundary shapes has been taken. This has been done for several almost local metrics. For these metrics it is straightforward to calculate the horizontal energy because the horizontal bundle equals the normal bundle. However, in the case of Sobolev type metrics the horizontal energy involves the inverse of a differential operator (see section 6.8), which makes this approach much harder. The second approach is the method of geodesic shooting. This method is based on iteratively solving the initial value problem while suitably adapting the initial   conditions. The theoretical requirements of existence of solutions to the geodesic equation and smooth dependence on initial conditions are met for Sobolev type metrics, see section 6.6. This makes geodesic shooting a promising approach.
The first step towards this aim is to numerically solve the initial value problem for geodesics, at least for the H 1 -metric and the case of surfaces in R 3 , and that is what will be presented in this work.
The geodesic equation on shape space is equivalent to the horizontal geodesic equation on the space of immersions. For the case of surfaces in R 3 , it takes the form given in section 9.1. This equation can be conveniently set up using the DifferentialGeometry package incorporated in the computer algebra system Maple as demonstrated in figure 6. (The equations that have actually been solved were simplified by multiplying intermediate terms with suitable powers of vol(g), but for the sake of clearness this has not been included in the Maple code in figure 6.) Unfortunately, Maple (as of version 14) is not able to solve PDEs with more than one space variable numerically. Thus the equations were translated into Mathematica. The PDE was solved using the method of lines. Spatial discretization was done using an equidistant grid, and spatial derivatives were replaced by finite differences.
The time-derivative f t appears implicitly in the equation P f (f t ) = a.ν, and this remains so when the operator P f is replaced by finite differences.
The solver that has been used is the Implicit Differential-Algebraic (IDA) solver that is part of the SUNDIALS suite and is integrated into Mathematica. IDA uses backward differentiation of order 1 to 5 with variable step widths. Order 5 is standard and has also been used here. At each time step, the new value of f t is computed using some previous values of f , and then the new value of f is calculated from the equation P f (f t ) = a.ν. The dependence on f in this equation is of course highly nonlinear. A Newton method is used to solve it. This operation is quite costly and has to be done at every step, which is a main disadvantage of backward differentiation algorithms. Explicit methods are probably much better adapted to the problem. The implementation of an explicit solver is ongoing common work of the authors with Martins Bruveris and Colin Cotter.
In the examples that follow, f at time zero is a square [0, π] × [0, π] flatly embedded in R 3 . This is a manifold with boundary, but it can be seen as a part of a bigger closed manifold. Zero boundary conditions are used for both f and a. It remains to specify an initial condition for a. As a first example, let us assume that a at time zero equals sin(x) sin(y), where x, y are the Euclidean coordinates on the square. The resulting geodesic is depicted in figure 1. In the absence of a closed-form solution of the geodesic equation, one way to check if the solution is correct is to see if the energy G P (f t , f t ) is conserved. Figure 2 shows this for the geodesic from figure 1 with various space and time discretizations.
A more complicated example of a geodesic is shown in figure 3 and 4. There, the initial velocity was chosen to be a smoothened version of a black and white picture of the letter A. The initial momentum was computed from it using a discrete Fourier transform.
Finally, it is shown in figure 5 that self-intersections of the surface can occur. This is not due to a numerical error but part of the theory, and can be an advantage or a disadvantage depending on the application.