Anisotropically Weighted and Nonholonomically Constrained Evolutions on Manifolds

We present evolution equations for a family of paths that results from anisotropically weighting curve energies in non-linear statistics of manifold valued data. This situation arises when performing inference on data that have non-trivial covariance and are anisotropic distributed. The family can be interpreted as most probable paths for a driving semi-martingale that through stochastic development is mapped to the manifold. We discuss how the paths are projections of geodesics for a sub-Riemannian metric on the frame bundle of the manifold, and how the curvature of the underlying connection appears in the sub-Riemannian Hamilton-Jacobi equations. Evolution equations for both metric and cometric formulations of the sub-Riemannian metric are derived. We furthermore show how rank-deficient metrics can be mixed with an underlying Riemannian metric, and we relate the paths to geodesics and polynomials in Riemannian geometry. Examples from the family of paths are visualized on embedded surfaces, and we explore computational representations on finite dimensional landmark manifolds with geometry induced from right-invariant metrics on diffeomorphism groups.


Introduction
When manifold valued data have non-trivial covariance (i.e., when anisotropy asserts higher variance in some directions than others), non-zero curvature necessitates special care when generalizing Euclidean space normal distributions to manifold valued distributions: in the Euclidean situation, normal distributions can be seen as transition distributions of diffusion processes, but on the manifold, holonomy makes transport of covariance path-dependent in the presence of curvature, preventing a global notion of a spatially constant covariance matrix. To handle this, in the diffusion principal component analysis (PCA) framework [1], and with the class of anisotropic normal distributions on manifolds defined in [2,3], data on non-linear manifolds are modelled as being distributed according to transition distributions of anisotropic diffusion processes that are mapped from Euclidean space to the manifold by stochastic development (see [4]). The construction is connected to a non-bracket-generating sub-Riemannian metric on the bundle of linear frames of the manifold, the frame bundle, and the requirement that covariance stays covariantly constant gives a nonholonomically constrained system. Velocity vectors and geodesic distances are conventionally used for estimation and statistics in Riemannian manifolds; for example, for estimation of the Frechét mean [5], for Principal Geodesic Analysis [6], and for tangent space statistics [7]. In contrast to this, anisotropy as modelled with anisotropic normal distributions makes a distance for a sub-Riemannian metric the natural vehicle for estimation and statistics. This metric naturally accounts for anisotropy in a similar way as the precision matrix weights the inner product in the negative log-likelihood of a Euclidean normal distribution. The connection between the weighted distance and statistics of manifold valued data was presented in [2], and the underlying sub-Riemannian and fiber-bundle geometry, together with properties of the generated densities, was further explored in [3]. The fundamental idea is to perform statistics on manifolds by maximum likelihood (ML) instead of parametric constructions that use, for example, approximating geodesic subspaces; by defining natural families of probability distributions (in this case using diffusion processes), ML parameter estimates give a coherent way to statistically model non-linear data. The anisotropically weighted distance and the resulting family of extremal paths arises in this situation when the diffusion processes have non-isotropic covariance (i.e., when the distribution is not generated from a standard Brownian motion).
In this paper, we focus on the family of most probable paths for the semi-martingales that drives the stochastic development, and in turn the manifold valued anisotropic stochastic processes. Such paths, as exemplified in Figure 1, extremize the anisotropically weighted action functional. We present derivations of evolution equations for the paths from different viewpoints, and we discuss the role of frames as representing either metrics or cometrics. In the derivation, we explicitly see the influence of the connection and its curvature. We then turn to the relation between the sub-Riemannian metric and the Sasaki-Mok metric on the frame bundle, and we develop a construction that allows the sub-Riemannian metric to be defined as a sum of a rank-deficient generator and an underlying Riemannian metric. Finally, we relate the paths to geodesics and polynomials in Riemannian geometry, and we explore computational representations on different manifolds including a specific case: the finite dimensional manifolds arising in the Large Deformation Diffeomorphic Metric Mapping (LDDMM) [8] landmark matching problem. The paper ends with a discussion concerning statistical implications, open questions, and concluding remarks. A frame u t (black/gray vectors) representing the square root covariance is parallel transported along the curve, enabling the anisotropic weighting with the precision matrix in the action functional. With isotropic covariance, normal MPPs are Riemannian geodesics. In general situations, such as the displayed anisotropic case, the family of MPPs is much larger; (b) The corresponding anti-development in R 2 (red line) of the MPP. Compare with the anti-development of a Riemannian geodesic with same initial velocity (blue dotted line). The frames u t ∈ GL(R 2 , T x t M) provide local frame coordinates for each time t.

Background
Generalizing common statistical tools for performing inference on Euclidean space data to manifold valued data has been the subject of extensive work (e.g., [9]). Perhaps most fundamental is the notion of Frechét or Karcher means [5,10], defined as minimizers of the square Riemannian distance. Generalizations of the Euclidean principal component analysis procedure to manifolds are particularly relevant for data exhibiting anisotropy. Approaches include principal geodesic analysis (PGA, [6]), geodesic PCA (GPCA, [11]), principal nested spheres (PNS, [12]), barycentric subspace analysis (BSA, [13]), and horizontal component analysis (HCA, [14]). Common to these constructions are explicit representations of approximating low-dimensional subspaces. The fundamental challenge here is that the notion of Euclidean linear subspace on which PCA relies has no direct analogue in non-linear spaces.
A different approach taken by diffusion PCA (DPCA, [1,2]) and probabilistic PGA [15] is to base the PCA problem on a maximum likelihood fit of normal distributions to data. In Euclidean space, this approach was first introduced with probabilistic PCA [16]. In DPCA, the process of stochastic development [4] is used to define a class of anisotropic distributions that generalizes the family of Euclidean space normal distributions to the manifold context. DPCA is then a simple maximum likelihood fit in this family of distributions mimicking the Euclidean probabilistic PCA. The approach transfers the geometric complexities of defining subspaces common in the approaches listed above to the problem of defining a geometrically natural notion of normal distributions.
In Euclidean space, squared distances x − x 0 2 between observations x and the mean x 0 are affinely related to the negative log-likelihood of a normal distribution N (x 0 , Id). This makes an ML fit of the mean such as performed in probabilistic PCA equivalent to minimizing squared distances. On a manifold, distances d g (x, x 0 ) 2 coming from a Riemannian metric g are equivalent to tangent space distances Log x 0 x 2 when mapping data from M to T x 0 M using the inverse exponential map Log x 0 . Assuming Log x 0 x are distributed according to a normal distribution in the linear space T x 0 M, this restores the equivalence with a maximum likelihood fit. Let {e 1 , . . . , e d } be the standard basis for R d . If u : R d → T x 0 M is a linear invertible map with ue 1 , . . . , ue d orthonormal with respect to g, the normal distribution in T x 0 M can be defined as uN (0, Id) (see Figure 2).
(a) Normal distributions uN (0, Id) in the tangent space T x 0 M with covariance uu T (blue ellipsis) can be mapped to the manifold by applying the exponential map Exp x 0 to sampled vectors v ∈ T x 0 M (red vectors). This effectively linearises the geometry around x 0 ; (b) The stochastic development map ϕ u maps R d valued paths w t to M by transporting the covariance in each step (blue ellipses) giving a covariance u t along the entire sample path. The approach does not linearise around a single point. Holonomy of the connection implies that the covariance "rotates" around closed loops-an effect which can be illustrated by continuing the transport along the loop created by the dashed path. The anisotropic metric g FM weights step lengths by the transported covariance at each time t.
The map u can be represented as a point in the frame bundle FM of M. When the orthonormal requirement on u is relaxed so that uN (0, Id) is a normal distribution in T x 0 M with anisotropic covariance, the negative log-likelihood in T x 0 M is related to (u −1 Log x 0 x) T (u −1 Log x 0 x) in the same way as the precision matrix Σ −1 is related to the negative log-likelihood (x − x 0 ) T Σ −1 (x − x 0 ) in Euclidean space. The distance is thus weighted by the anisotropy of u, and u can be interpreted as a square root covariance matrix Σ 1/2 . However, the above approach does not specify how u changes when moving away from the base point x 0 . The use of Log x 0 x effectively linearises the geometry around x 0 , but a geometrically natural way to relate u at points nearby to x 0 will be to parallel transport it, equivalently specifying that u when transported does not change as measured from the curved geometry. This constraint is nonholonomic, and it implies that any path from x 0 to x carries with it a parallel transport of u lifting paths from M to paths in the frame bundle FM. It therefore becomes natural to equip FM with a form of metric that encodes the anisotropy represented by u. The result is the sub-Riemannian metric on FM defined below that weights infinitesimal movements on M using the parallel transport of the frame u. Optimal paths for this metric are sub-Riemannian geodesics giving the family of most probable paths for the driving process that this paper concerns. Figure 1 shows one such path for an anisotropic normal distribution with M an ellipsoid embedded in R 3 .

Frame Bundles, Stochastic Development, and Anisotropic Diffusions
Let M be a finite dimensional manifold of dimension d with connection C, and let x 0 be a fixed point in M. When a Riemannian metric is present, and C is its Levi-Civita connection, we denote the metric g R . For a given interval [0, T], we let W(M) denote the Wiener space of continuous paths in M starting at x 0 . Similarly, W(R d ) is the Wiener space of paths in R d . We let H(R d ) denote the subspace of W(R d ) of finite energy paths.
Let now u = (u 1 , . . . , u d ) be a frame for T x M, x ∈ M; i.e., u 1 , . . . , u d is an ordered set of linearly independent vectors in T x M with span{u 1 , . . . , u d } = T x M. We can regard the frame as an isomorphism u : R d → T x M with u(e i ) = u i , where e 1 , . . . , e d denotes the standard basis in R d . Stochastic development (e.g., [4]) provides an invertible map ϕ u from W(R d ) to W(M). Through ϕ u , Euclidean semi-martingales map to stochastic processes on M. When M is Riemannian and u orthonormal, the result is the Eells-Elworthy-Malliavin construction of Brownian motion [17]. We here outline the geometry behind development, stochastic development, the connection, and curvature, focusing in particular on frame bundle geometry.

The Frame Bundle
For each point x ∈ M, let F x M be the set of frames for T x M (i.e., the set of ordered bases for T x M). The set {F x M} x∈M can be given a natural differential structure as a fiber bundle on M called the frame bundle FM. It can equivalently be defined as the principal bundle GL(R d , TM). We let the map π : FM → M denote the canonical projection. The kernel of π * : TFM → TM is the sub-bundle of TFM that consists of vectors tangent to the fibers π −1 (x). It is denoted the vertical subspace VFM. We will often work in a local trivialization u = (x, u 1 , . . . , u d ) ∈ FM, where x = π(u) ∈ M denotes the base point, and for each i = 1, . . . , d, u i ∈ T x M is the ith frame vector. For v ∈ T x M and u ∈ FM with π(u) = x, the vector u −1 v ∈ R d expresses v in components in terms of the frame u. We will denote the vector u −1 v frame coordinates of v.
For a differentiable curve x t in M with x = x 0 , a frame u for T x 0 M can be parallel transported along x t by parallel transporting each vector in the frame, thus giving a path u t ∈ FM. Such paths are called horizontal, and have zero acceleration in the sense C(u i,t ) = 0. For each x ∈ M, their derivatives form a d-dimensional subspace of the d + d 2 -dimensional tangent space T u FM. This horizontal subspace HFM and the vertical subspace VFM together split the tangent bundle of FM (i.e., TFM = HFM ⊕ VFM). The split induces a map π * : HFM → TM, see Figure 3. For fixed u ∈ FM, the restriction π * | H u FM : H u FM → T x M is an isomorphism. Its inverse is called the horizontal lift and is denoted h u in the following. Using h u , horizontal vector fields H e on FM are defined for vectors e ∈ R d by H e (u) = h u (ue). In particular, the standard basis (e 1 , . . . , e d ) on R d gives d globally defined horizontal vector fields H i ∈ HFM, i = 1, . . . , d by H i = H e i . Intuitively, the fields H i (u) model infinitesimal transformations in M of x 0 in direction u i = ue i with corresponding infinitesimal parallel transport of the vectors u 1 , . . . , u d of the frame along the direction u i . A horizontal lift of a differentiable curve x t ∈ M is a curve in FM tangent to HFM that projects to x t . Horizontal lifts are unique up to the choice of initial frame u 0 . Figure 3. Relations between the manifold, frame bundle, the horizontal distribution HFM, the vertical bundle VFM, a Riemannian metric g R , and the sub-Riemannian metric g FM , defined below.
The connection C provides the splitting TFM = HFM ⊕ VFM. The restrictions π * | H u M are invertible maps H u M → T π(u) M with inverse h u , the horizontal lift. Correspondingly, the vertical bundle VFM is isomorphic to the trivial bundle FM × gl(n). The metric g FM : T * FM → TFM has an image in the subspace HFM.

Development and Stochastic Development
Let x t be a differentiable curve on M and u t a horizontal lift. If s t is a curve in R d with components s i t such thatẋ t = H i (u)s i t , x t is said to be a development of s t . Correspondingly, s t is the anti-development of x t . For each t, the vector s t contains frame coordinates ofẋ t as defined above. Similarly, let W t be an R d valued Brownian motion so that sample paths We call the process W t in R d , that through ϕ maps to X t , the driving process of X t . Let ϕ u : W(R d ) → W(M) be the map that for fixed u sends a path in R d to its development on M. Its inverse ϕ −1 u is the anti-development in R d of paths on M given u.
Equivalent to the fact that normal distributions N (0, Σ) in R d can be obtained as the transition distributions of diffusion processes Σ 1/2 W t stopped at time t = 1, a general class of distributions on the manifold M can be defined by stochastic development of processes W t , resulting in M-valued random variables X = X 1 . This family of distributions on M introduced in [2] is denoted anisotropic normal distributions. The stochastic development by construction ensures that U t is horizontal, and the frames are thus parallel transported along the stochastic displacements. The effect is that the frames stay covariantly constant, thus resembling the Euclidean situation where Σ 1/2 is spatially constant and therefore does not change as W t evolves. Thus, as further discussed in Section 3.2, the covariance is kept constant at each of the infinitesimal stochastic displacements. The existence of a smooth density for the target process X t and small time asymptotics are discussed in [3].
Stochastic development gives a map Diff : FM → Prob(M) to the space of probability distributions on M. For each point u ∈ FM, the map sends a Brownian motion in R d to a distribution µ u by stochastic development of the process U t in FM, starting at u and letting µ u be the distribution of X = π(U 1 ). The pair (x, u), x = π(u) is analogous to the parameters (µ, Σ) for a Euclidean normal distribution: the point x ∈ M represents the starting point of the diffusion, and the frame u represents a square root Σ 1/2 of the covariance Σ. In the general situation where µ u has smooth density, the construction can be used to fit the parameters u to data by maximum likelihood. As an example, diffusion PCA fits distributions obtained through Diff by maximum likelihood to observed samples in M; i.e., it optimizes for the most likely parameters u = (x, u 1 , . . . , u d ) for the anisotropic diffusion process, giving a fit to the data of the manifold generalization of the Euclidean normal distribution.

Adapted Coordinates
For concrete expressions of the geometric constructions related to frame bundles, and for computational purposes, it is useful to apply coordinates that are adapted to the horizontal bundle HFM and the vertical bundle VFM together with their duals H * FM and V * FM. The notation below follows the notation used in, for example, [18]. Let z = (u, ξ) be a local trivialization of T * FM, and let To find a basis that is adapted to the horizontal distribution, define the d linearly independent Switching between standard coordinates and the adapted frame and coframes can be expressed in terms of the component matrices A below the frame and coframe induced by the coordinates (x i , u i α ) and the adapted frame D and coframe D * . We have

Connection and Curvature
The TM valued connection C : TM × TM → TM lifts to a principal connection TFM × TFM → VFM on the principal bundle FM. C can then be identified with the gl(n)-valued connection form ω on TFM. The identification occurs by the isomorphism ψ between FM × gl(n) and VFM given by ψ(u, v) = d dt u exp(tv)| t=0 (e.g., [19,20]). The map ψ is equivariant with respect to the GL(n) action g → ug −1 on FM. In order to explicitly see the connection between the usual covariant derivative ∇ : Γ(TM) × Γ(TM) → Γ(TM) on M determined by C and C regarded as a connection on the principal bundle FM, following [19], we let s : M → TM be a local vector field on M; equivalently, s ∈ Γ(TM) is a local section of TM. s determines a map s FM : FM → R d by s FM (u) = u −1 s(π(u)); i.e., it gives the coordinates of s(x) in the frame u at x. The pushforward (s FM ) * : TFM → R d has in its ith component the exterior derivative d(s FM ) i . Let now w(x) be a local section of FM. The composition w • (s FM ) * • h w : TM → TM is identical to the covariant derivative ∇ · s : TM → TM. The construction is independent of the choice of w because of the GL(n)-equivariance of s FM . The connection form ω can be expressed as the matrix The identification becomes particularly simple if the covariant derivative is taken along a curve x t on which w t is the horizontal lift. In this case, we can let i.e., the covariant derivative takes the form of the standard derivative applied to the frame coordinates s i t . The curvature tensor R ∈ T 3 1 (M) gives the gl(n)-valued curvature form Ω : Note that Ω(v u , w u ) = Ω(h u (π * (v u )), h u (π * (w u ))), which we can use to write the curvature R as the gl(n)-valued map R u : In coordinates, the curvature is s be a family of paths in M, and let u t,s ∈ π −1 (x t,s ) be horizontal lifts of x t,s for each fixed s. Writeẋ t,s = ∂ t x t,s andu t,s = ∂ t u t,s . The s-derivative of u t,s can be regarded a pushforward of the horizontal lift and is the curve in TFM This follows from the structure equation dω = −ω ∧ ω + Ω (e.g., [21]). Note that the curve depends on the vertical variation C(∂ s u 0,s ) at only one point along the curve. The remaining terms depend on the horizontal variation or, equivalently, ∂ s x t,s . The t-derivative of ∂ s u t,s is the curve in TTFM satisfying Here, the first and third term in the last expression are identified with elements of T ∂ s u t,s TFM by the natural mapping T u t,s FM → T ∂ s u t,s TFM. When C(∂ s u 0,s ) is zero, the relation reflects the property that the curvature arises when computing brackets between horizontal vector fields. Note that the first term of (3) has values in VFM, while the third term has values in HFM.

The Anisotropically Weighted Metric
For a Euclidean driftless diffusion process with spatially constant stochastic generator Σ, the log-probability of a sample path can formally be written with the norm · Σ given by the inner product v, w , the inner product weighted by the precision matrix Σ −1 . Though only formal, as the sample paths are almost surely nowhere differentiable, the interpretation can be given a precise meaning by taking limits of piecewise linear curves [21]. Turning to the manifold situation with the processes mapped to M by stochastic development, the probability of observing a differentiable path can either be given a precise meaning in the manifold by taking limits of small tubes around the curve, or in R d by considering infinitesimal tubes around the anti-development of the curves. With the former formulation, a scalar curvature correction term must be added to (4), giving the Onsager-Machlup function [22]. The latter formulation corresponds to defining a notion of path density for the driving R d -valued process W t . When M is Riemannian and Σ unitary, taking the maximum of (4) gives geodesics as most probable paths for the driving process. Let now u t be a path in FM, and choose a local trivialization u t = (x t , u 1,t , . . . , u d,t ) such that the matrix [u i α,t ] represents the square root covariance matrix Σ 1/2 at x t . Since u t being a frame defines an invertible map R d → T x t M, the norm · Σ above has a direct analogue in the norm · u t defined by the inner product v, for vectors v, w ∈ T x t M. The transport of the frame along paths in effect defines a transport of inner product along sample paths: the paths carry with them the inner product weighted by the precision matrix, which in turn is a transport of the square root covariance u 0 at x 0 . The inner product can equivalently be defined as a metric g u : T * x M → T x M. Again using that u can be considered a map R d → T x M, g u is defined by ξ → u ( (ξ • u) ), where is the standard identification (R d ) * → R d . The sequence of mappings defining g u is illustrated below: This definition uses the R d inner product in the definition of . Its inverse gives the cometric

Sub-Riemannian Metric on the Horizontal Distribution
We now lift the path-dependent metric defined above to a sub-Riemannian metric on HFM. For any w,w ∈ H u FM, the lift of (5) by π * is the inner product The inner product induces a sub-Riemannian metric g FM : TFM * → HFM ⊂ TFM by w, g FM (ξ) = (ξ|w) , ∀w ∈ H u FM with (ξ|w) denoting the evaluation ξ(w) for the covector ξ ∈ T * FM. The metric g FM gives FM a non-bracket-generating sub-Riemannian structure [23] on FM (see also Figure 3). It is equivalent to the lift of the metric g u above. In frame coordinates, the metric takes the form In terms of the adapted coordinates for TFM described in Section 2.3, with w = w j D j and w =w j D j , we have Define now W kl = δ αβ u k α u l β , so that W ir W rj = δ i j and W ir W rj = δ j i . We can then write the metric g FM directly as One clearly recognizes the dependence on the horizontal H * FM part of T * FM only, and the fact that g FM has image in HFM. The sub-Riemannian energy of an almost everywhere horizontal path u t is i.e., the line element is ds 2 = W ij D i D j in adapted coordinates. The corresponding distance is given by If we wish to express g FM in canonical coordinates on T * FM, we can switch between the adapted frame and the coordinates ( (11), g FM has D, D * components Therefore, g FM has the following components in the coordinates (

Covariance and Nonholonomicity
The metric g FM encodes the anisotropic weighting given the frame u, thus up to an affine transformation measuring the energy of horizontal paths equivalently to the negative log-probability of sample paths of Euclidean anisotropic diffusions as formally given in (4). In addition, the requirement that paths must stay horizontal almost everywhere enforces that C(u t ) = 0 a.e., i.e., that no change of the covariance is measured by the connection. The intuitive effect is that covariance is covariantly constant as seen by the connection. Globally, curvature of C will imply that the covariance changes when transported along closed loops, and torsion will imply that the base point "slips" when travelling along covariantly closed loops on M. However, the zero acceleration requirement implies that the covariance is as close to spatially constant as possible with the given connection. This is enabled by the parallel transport of the frame, and it ensures that the model closely resembles the Euclidean case with spatially constant stochastic generator.
With non-zero curvature of C, the horizontal distribution is non-integrable (i.e., the brackets [H i , H j ] are non-zero for some i, j). This prevents integrability of the horizontal distribution HFM in the sense of the Frobenius theorem. In this case, the horizontal constraint is nonholonomic similarly to nonholonomic constraints appearing in geometric mechanics (e.g., [24]). The requirement of covariantly constant covariance thus results in a nonholonomic system.

Riemannian Metrics on FM
If the horizontality constraint is relaxed, a related Riemannian metric on FM can be defined by pulling back a metric on gl(n) to each fiber using the isomorphism ψ(u, ·) −1 : V u FM → gl(n). Therefore, the metric on HFM can be extended to a Riemannian metric on FM. Such metrics incorporate the anisotropically weighted metric on HFM, however, allowing vertical variations and thus that covariances can change unrestricted.
When M is Riemannian, the metric g FM is in addition related to the Sasaki-Mok metric on FM [18] that extends the Sasaki metric on TM. As for the above Riemannian metric on FM, the Sasaki-Mok metric allows paths in FM to have derivatives in the vertical space VFM. On HFM, the Riemannian metric g R is here lifted to the metric g SM = (v u , w u ) = g R (π * (v u ), π * (w u )) (i.e., the metric is not anisotropically weighted). The line element is in this case ds 2 = g ij dx i dx j + X βα g ij D α i D β j .
Geodesics for g SM are lifts of Riemannian geodesics for g R on M, in contrast to the sub-Riemannian normal geodesics for g FM which we will characterize below. The family of curves arising as projections to M of normal geodesics for g FM includes Riemannian geodesics for g R (and thus projections of geodesics for g SM ), but the family is in general larger than geodesics for g R .

Constrained Evolutions
Extremal paths for (5) can be interpreted as most probable paths for the driving process W t when u 0 defines an anisotropic diffusion. This is captured in the following definition [3]: Definition 1. A most probable path for the driving process (MPP) from x = π(u 0 ) ∈ M to y ∈ M is a smooth path x t : [0, 1] → M with x 0 = x and x 1 = y such that its anti-development ϕ −1 u 0 (x t ) is a most probable path for W t ; i.e., with L R d being the Onsager-Machlup function for the process W t on R d [22].
The definition uses the one-to-one relation between W(R d ) and W(M) provided by ϕ u 0 to characterize the paths using the R d Onsager-Machlup function L R d . When M is Riemannian with metric g R , the Onsager-Machlup function for a g-Brownian motion on M is L(x t ,ẋ t ) = − 1 2 ẋ t 2 g R + 1 12 S g R (x t ) with S g R denoting the scalar curvature. This curvature term vanishes on R d , and therefore L R d (γ t ,γ t ) = − 1 2 γ t 2 for a curve γ t ∈ R d . By pulling x t ∈ M back to R d using ϕ −1 u 0 , the construction removes the 1 12 S g R (x t ) scalar curvature correction term present in the non-Euclidean Onsager-Machlup function. It thereby provides a relation between geodesic energy and most probable paths for the driving process. This is contained in the following characterization of most probable paths for the driving process as extremal paths of the sub-Riemannian distance [3] that follows from the Euclidean space Onsager-Machlup theorem [22].

Theorem 1 ([3]
). Let Q(u 0 ) denote the principal sub-bundle of FM of points z ∈ FM reachable from u 0 ∈ FM by horizontal paths. Suppose the Hörmander condition is satisfied on Q(u 0 ), and that Q(u 0 ) has compact fibers. Then, most probable paths from x 0 to y ∈ M for the driving process of X t exist, and they are projections of sub-Riemannian geodesics in FM minimizing the sub-Riemannian distance from u 0 to π −1 (y).
Below, we will derive evolution equations for the set of such extremal paths that correspond to normal sub-Riemannian geodesics.

Connected to the metric g FM is the Hamiltonian
on the symplectic space T * FM. Lettingπ denote the projection on the bundle T * FM → FM, (8) gives Normal geodesics in sub-Riemannian manifolds satisfy the Hamilton-Jacobi equations [23] with Hamiltonian flowż where Ω here is the canonical symplectic form on T * FM (e.g., [25]). We denote (13) the MPP equations, and we let projections x t = π T * FM (z t ) of minimizing curves satisfying (13) be denoted normal MPPs. The system (13) has 2(d + d 2 ) degrees of freedom, in contrast to the usual 2d degrees of freedom for the classical geodesic equation. Of these, d 2 describes the current frame at time t, while the remaining d 2 allows the curve to "twist" while still being horizontal. We will see this effect visualized in Section 6.
In a local canonical trivialization z = (u, ξ), (13) gives the Hamilton-Jacobi equationṡ Using (3), we have for the second equatioṅ Here ∂u denotes u-derivative in the directionu, equivalently ∂uh u (v) = ∂ t (h u )(v). While the first equation of (14) involves only the horizontal part of ξ, the second equation couples the vertical part of ξ through the evaluation of ξ on the term ψ(u, R u (π * (u), π * (∂ u )). If the connection is curvature-free, which in non-flat cases implies that it carries torsion, this vertical term vanishes. Conversely, when M is Riemannian, C the g R Levi-Civita connection, and u 0 is for all v, w ∈ T π(u t ) M. In this case, a normal MPP π(u t ) will be a Riemannian g R geodesic.

Evolution in Coordinates
In coordinates u = (x i , u i α , ξ i , ξ i α ) for T * FM, we can equivalently writė Combining these expressions, we obtaiṅ

Acceleration and Polynomials for C
We can identify the covariant acceleration ∇ẋ tẋ t of curves satisfying the MPP equations, and hence normal MPPs through their frame coordinates. Let (u t , ξ t ) satisfy (13). Then, u t is a horizontal lift of x t = π(u t ) and hence by (1), (3), (10), and (15), The fact that the covariant derivative vanishes for classical geodesic leads to a definition of higher-order polynomials through the covariant derivative by requiring (∇ẋ t ) kẋ t = 0 for a kth order polynomial (e.g., [26,27]). As discussed above, compared to classical geodesics, curves satisfying the MPP equations have extra d 2 degrees of freedom, allowing the curves to twist and deviate from being geodesic with respect to C while still satisfying the horizontality constraint on FM. This makes it natural to ask if normal MPPs relate to polynomials defined using C. For curves satisfying the MPP equations, using (16) and (15), we have Thus, in general, normal MPPs are not second order polynomials in the sense (∇ẋ t ) 2ẋ t = 0 unless the curvature R u t (u t e i , π * (u t )) is constant in t.
For comparison, in the Riemannian case, a variational formulation placing a cost on covariant acceleration [28,29] leads to cubic splines In (16), the curvature terms appear in the covariant acceleration for normal MPPs, while cubic splines leads to the curvature term appearing in the third order derivative.

Cometric Formulation and Low-Rank Generator
We now investigate a cometric g F k M + λg R , where g R is Riemannian, g F k M is a rank k positive semi-definite inner product arising from k linearly independent tangent vectors, and λ > 0 a weight. We assume that g F k M is chosen so that g F k M + λg R is invertible, even though g F k M is rank-deficient.
The situation corresponds to extracting the first k eigenvectors in Euclidean space PCA. If the eigenvectors are estimated statistically from observed data, this allows the estimation to be restricted to only the first k eigenvectors. In addition, an important practical implication of the construction is that a numerical implementation need not transport a full d × d matrix for the frame, but a potentially much lower dimensional d × k matrix. This point is essential when dealing with high-dimensional data, examples of which are landmark manifolds as discussed in Section 6.
When using the frame bundle to model covariances, the sum formulation is natural to express as a cometric compared to a metric because, with the cometric formulation, g F k M + λg R represents a sum of covariance matrices instead of a sum of precision matrices. Thus, g F k M + λg R can be intuitively thought of as adding isotropic noise of variance λ to the covariance represented by g F k M .
To pursue this, let F k M denote the bundle of rank k linear maps R k → T x M. We define a cometric by ξ,ξ = δ αβ (ξ|h u (u α ))(ξ|h u (u β )) + λ ξ,ξ g R for ξ,ξ ∈ T * u F k M. The sum over α, β is for α, β = 1, . . . , k. The first term is equivalent to the lift (9) of the cometric ξ,ξ = ξ|g u (ξ) given u : R k → T x M. Note that in the definition (6) of g u , the map u is not inverted; thus, the definition of the metric immediately carries over to the rank-deficient case.
Let (x i , u i α ), α = 1, . . . , k be a coordinate system on F k M. The vertical distribution is in this case spanned by the dk vector fields D j β = ∂ u j β . Except for index sums being over k instead of d terms, the situation is thus similar to the full-rank case. Note that (ξ|π −1 * w) = (ξ|w j D j ) = w i ξ i . The cometric in coordinates is We can then write the corresponding sub-Riemannian metric g F k M in terms of the adapted frame D because (ξ|g F k M (ξ)) = ξ,ξ = ξ i W ijξ j . That is, the situation is analogous to (11), except the term λg ij R is added to W ij . The geodesic system is again given by the Hamilton-Jacobi equations. As in the full-rank case, the system is specified by the derivatives of g F k M : hj ,l , Note that the introduction of the Riemannian metric g R implies that W ij are now dependent on the manifold coordinates x i .

Numerical Experiments
We aim at visualizing most probable paths for the driving process and projections of curves satisfying the MPP Equation (13) in two cases: On 2D surfaces embedded in R 3 and on finite dimensional landmark manifolds that arise from equipping a subset of the diffeomorphism group with a right-invariant metric and letting the action descend to the landmarks by a left action. The surface examples are implemented in Python using the Theano [30] framework for symbolic operations, automatic differentiation, and numerical evaluation. The landmark equations are detailed below and implemented in Numpy using Numpy's standard ODE integrators. The code for running the experiments is available at http://bitbucket.com/stefansommer/mpps/.

Embedded Surfaces
We visualize normal MPPs and projections of curves satisfying the MPP Equation (13) 1) With fixed starting point x 0 ∈ M and initial velocityẋ 0 ∈ TM but varying anisotropy represented by changing frame u in the fiber above x 0 ; (2) minimizing normal MPPs with fixed starting point and endpoint x 0 , x 1 ∈ M but changing frame u above x 0 ; (3) fixed starting point x 0 ∈ M and frame u but varying V * FM vertical part of the initial momentum ξ 0 ∈ T * FM. The first and second cases thus show the effect of varying anisotropy, while the third case illustrates the effect of the "twist" that the d 2 degrees in the vertical momentum allows. Note the displayed anti-developed curves in R 2 that for classical C geodesics would always be straight lines.

LDDMM Landmark Equations
We here give a example of the MPP equations using the finite dimensional landmark manifolds that arise from right invariant metrics on subsets of the diffeomorphism group in the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework [8]. The LDDMM metric can be conveniently expressed as a cometric, and, using a rank-deficient inner product g F k M as discussed in Section 5, we can obtain a reduction of the system of equations to 2(2N + 2Nk) compared to 2(2N + (2N) 2 ) with N landmarks in R 2 .
Under sufficient conditions on L, V is reproducing and admits a kernel K inverse to L. K is a Green's kernel when L is a differential operator, or K can be a Gaussian kernel. The Hilbert structure on V gives a Riemannian metric on a subset G V ⊂ Diff(Ω) by setting v 2 ϕ = v • ϕ −1 2 V ; i.e., regarding ·, · V an inner product on T Id G V and extending the metric to G V by right-invariance. This Riemannian metric descends to a Riemannian metric on the landmark space.
Let M be the manifold M = {(p 1 1 , . . . , p d 1 , . . . , The LDDMM metric on the landmark manifold M is directly related to the kernel K when written as a cometric g p (ξ, η) = ∑ N i,j=1 ξ i K(p i , p j )η j . Letting i k denote the index of the kth component of the ith landmark, the cometric is in coordinates g i k j l p = K(p i , p j ) l k . The Christoffel symbols can be written in terms of derivatives of the cometric g ij [31] (recall that δ i j = g ik g kj = g jk g ki ) Γ k ij = 1 2 g ir g kl g rs ,l − g sl g rk ,l − g rl g ks ,l g sj .
This relation comes from the fact that g jm,k = −g jr g rs ,k g sm gives the derivative of the metric. The derivatives of the cometric is simply g i k j l ,r q = (δ i r + δ j r )∂ p q r K(p i , p j ) l k . Using (18), derivatives of the Christoffel symbols can be computed Γ k ij,ξ = 1 2 g ir,ξ g kl g rs ,l − g sl g rk ,l − g rl g ks ,l g sj + 1 2 g ir g kl g rs ,l − g sl g rk ,l − g rl g ks ,l g sj,ξ + 1 2 g ir g kl ,ξ g rs ,l + g kl g rs ,lξ − g sl ,ξ g rk ,l − g sl g rk ,lξ − g rl ,ξ g ks ,l − g rl g ks ,lξ g sj .
This provides the full data for numerical integration of the evolution equations on F k M. In Figure 7 (top row), we plot minimizing normal MPPs on the landmark manifold with two landmarks and varying covariance in the R 2 horizontal and vertical direction. The plot shows the landmark equivalent of the experiment in Figure 5. Note how adding covariance in the horizontal and vertical direction, respectively, allows the minimizing normal MPP to vary more in these directions because the anisotropically-weighted metric penalizes high-covariance directions less. . The action of the corresponding diffeomorphisms on a regular grid is visualized by the deformed grid which is colored by the warp strain. The added covariance allows the paths to have more movement in the horizontal and vertical direction, respectively, because the anisotropically weighted metric penalizes high-covariance directions less. (bottom row, (f-j)) Five landmark trajectories with fixed initial velocity and anisotropic covariance but varying V * FM vertical initial momentum ξ 0 . Changing the vertical momentum "twists" the paths. Figure 7 (bottom row) shows five curves satisfying the MPP equations with varying vertical V * FM initial momentum similarly to the plots in Figure 6. Again, we see how the extra degrees of freedom allows the paths to twist, generating a higher-dimensional family than classical geodesics with respect to C.

Discussion and Concluding Remarks
Incorporating anisotropy in models for data in non-linear spaces via the frame bundle as pursued in this paper leads to a sub-Riemannian structure and metric. A direct implication is that most probable paths to observed data in the sense of sequences of stochastic steps of a driving semi-martingale are not related to geodesics in the classical sense. Instead, a best estimate of the sequence of steps w t ∈ R d that leads to an observation x = ϕ u (w t )| t=1 is an MPP in the sense of Definition 1. As shown in the paper, these paths are generally not geodesics or polynomials with respect to the connection on the manifold. In particular, if M has a Riemannian structure, the MPPs are generally neither Riemannian geodesics nor Riemannian polynomials. Below, we discuss the statistical implications of this result.

Statistical Estimators
Metric distances and Riemannian geodesics have been the traditional vehicle for representing observed data in non-linear spaces. Most fundamentally, the sample Frechét mean of observed data x 1 , . . . , x N ∈ M relies crucially on the Riemannian distance d g R connected to the metric g R . Many PCA constructs (e.g., Principal Geodesics Analysis [6]) use the Riemannian Exp. and Log maps to map between linear tangent spaces and the manifold. These maps are defined from the Riemannian metric and Riemannian geodesics. Distributions modelled as in the random orbit model [32] or Bayesian models [15,33] again rely on geodesics with random initial conditions. Using the frame bundle sub-Riemannian metric g FM , we can define an estimator analogous to the Riemannian Frechét mean estimator. Assuming the covariance is a priori known, the estimator (20) acts correspondingly to the Frechét mean estimator (19). Here s ∈ Γ(FM) is a (local) section of FM that to x ∈ M connects the known covariance represented by s(x) ∈ FM. The distances d FM u, π −1 (x i ) , u = s(x) are realized by MPPs from the mean candidate x to the fibers π −1 (x i ). The Frechét mean problem is thus lifted to the frame bundle with the anisotropic weighting incorporated in the metric g FM . This metric is not related to g R , except for its dependence on the connection C that can be defined as the Levi-Civita connection of g R . The fundamental role of the distance d g R and g R geodesics in (19) is thus removed.
Because covariance is an integral part of the model, sample covariance can also be estimated directly along with the sample mean. In [3], the estimator is suggested. The normalizing term −N log(det g R u) is derived such that the estimator exactly corresponds to the maximum likelihood estimator of mean and covariance for Euclidean Gaussian distributions. The determinant is defined via g R , and the term acts to prevent the covariance from approaching infinity. Maximum likelihood estimators of mean and covariance for normally distributed Euclidean data have unique solutions in the sample mean and sample covariance matrix, respectively.
Uniqueness of the Frechét mean (19) is only ensured for sufficiently concentrated data. For the estimator (21), existence and uniqueness properties are not immediate, and more work is needed in order to find necessary and sufficient conditions.

Priors and Low-Rank Estimation
The low-rank cometric formulation pursued in Section 5 gives a natural restriction of (21) to u ∈ F k M, 1 ≤ k ≤ d. As for Euclidean PCA, most variance is often captured in the span of the first k eigenvectors with k d. Estimates of the remaining eigenvectors are generally ignored, as the variance of the eigenvector estimates increases as the noise captured in the span of the last eigenvectors becomes increasingly uniform. The low-rank cometric restricts the estimation to only the first k eigenvectors, and thus builds the construction directly into the model. In addition, it makes numerical implementation feasible, because a numerical representation need only store and evolve d × k matrices. As a different approach for regularizing the estimator (21), the normalizing term −N log(det g R u) can be extended with other priors (e.g., an L 1 -type penalizing term). Such priors can potentially partly remove existence and uniqueness issues, and result in additional sparsity properties that can benefit numerical implementations. The effects of such priors have yet to be investigated.
In the k = d case, the number of degrees of freedom for the MPPs grows quadratically in the dimension d. This naturally increases the variance of any MPP estimate given only one sample from its trajectory. The low-rank cometric formulation reduces the growth to linear in d. The number of degrees of freedom is however still k times larger than for Riemannian geodesics. With longitudinal data, more samples per trajectory can be obtained, reducing the variance and allowing a better estimate of the MPP. However, for the estimators (20) and (21) above, estimates of the actual optimal MPPs are not needed-only their squared length. It can be hypothesized that the variance of the length estimates is lower than the variance of the estimates of the corresponding MPPs. Further investigation regarding this will be the subject of future work.

Conclusions
The underlying model of anisotropy used in this paper originates from the anisotropic normal distributions formulated in [2] and the diffusion PCA framework [1]. Because many statistical models are defined using normal distributions, this approach to incorporating anisotropy extends to models such as linear regression. We expect that finding most probable paths in other statistical models such as regressions models can be carried out with a program similar to the program presented in this paper.
The difference between MPPs and geodesics shows that the geometric and metric properties of geodesics, zero acceleration, and local distance minimization are not directly related to statistical properties such as maximizing path probability. Whereas the concrete application and model determines if metric or statistical properties are fundamental, most statistical models are formulated without referring to metric properties of the underlying space. It can therefore be argued that the direct incorporation of anisotropy and the resulting MPPs are natural in the context of many models of data variation in non-liner spaces.