Anisotropically Weighted and Nonholonomically Constrained Evolutions on Manifolds ^†

Abstract

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.

1. 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.

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 $\parallel x-x_0{\parallel }^2$ between observations x and the mean $x_0$ are affinely related to the negative log-likelihood of a normal distribution $\mathcal{N}(x_0,\mathrm{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 $\parallel {\mathrm{Log}}_{x_0}{x\parallel }^2$ when mapping data from M to $T_{x_0}M$ using the inverse exponential map ${\mathrm{Log}}_{x_0}$. Assuming ${\mathrm{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,\dots ,e_d\}$ be the standard basis for ${\mathbb{R}}^d$. If $u:{\mathbb{R}}^d\to T_{x_0}M$ is a linear invertible map with $u{e}_1,\dots ,u{e}_d$ orthonormal with respect to g, the normal distribution in $T_{x_0}M$ can be defined as $u\mathcal{N}(0,\mathrm{Id})$ (see Figure 2).

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 $u\mathcal{N}(0,\mathrm{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 ${\left(u^{-1}{\mathrm{Log}}_{x_0}x\right)}^T\left(u^{-1}{\mathrm{Log}}_{x_0}x\right)$ in the same way as the precision matrix ${\Sigma }^{-1}$ is related to the negative log-likelihood ${(x-x_0)}^T{\Sigma }^{-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 ${\Sigma }^{1/2}$.

However, the above approach does not specify how u changes when moving away from the base point $x_0$. The use of ${\mathrm{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 ${\mathbb{R}}^3$.

2. Frame Bundles, Stochastic Development, and Anisotropic Diffusions

Let M be a finite dimensional manifold of dimension d with connection $\mathcal{C}$, and let $x_0$ be a fixed point in M. When a Riemannian metric is present, and $\mathcal{C}$ is its Levi–Civita connection, we denote the metric $g_R$. For a given interval $[0,T]$, we let $W\left(M\right)$ denote the Wiener space of continuous paths in M starting at $x_0$. Similarly, $W\left({\mathbb{R}}^d\right)$ is the Wiener space of paths in ${\mathbb{R}}^d$. We let $H\left({\mathbb{R}}^d\right)$ denote the subspace of $W\left({\mathbb{R}}^d\right)$ of finite energy paths.

Let now $u=(u_1,\dots ,u_d)$ be a frame for $T_{x}M$, $x\in M$; i.e., $u_1,\dots ,u_d$ is an ordered set of linearly independent vectors in $T_{x}M$ with $\mathrm{span}\{u_1,\dots ,u_d\}=T_{x}M$. We can regard the frame as an isomorphism $u:{\mathbb{R}}^d\to T_{x}M$ with $u\left(e_i\right)=u_i$, where $e_1,\dots ,e_d$ denotes the standard basis in ${\mathbb{R}}^d$. Stochastic development (e.g., [4]) provides an invertible map ${\phi }_u$ from $W\left({\mathbb{R}}^d\right)$ to $W\left(M\right)$. Through ${\phi }_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.

2.1. The Frame Bundle

For each point $x\in 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 ${\left\{F_{x}M\right\}}_{x\in 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 $\mathrm{GL}({\mathbb{R}}^d,TM)$. We let the map $\pi :FM\to M$ denote the canonical projection. The kernel of ${\pi }_{*}:TFM\to TM$ is the sub-bundle of $TFM$ that consists of vectors tangent to the fibers ${\pi }^{-1}\left(x\right)$. It is denoted the vertical subspace $VFM$. We will often work in a local trivialization $u=(x,u_1,\dots ,u_d)\in FM$, where $x=\pi \left(u\right)\in M$ denotes the base point, and for each $i=1,\dots ,d$, $u_i\in T_{x}M$ is the ith frame vector. For $v\in T_{x}M$ and $u\in FM$ with $\pi \left(u\right)=x$, the vector $u^{-1}v\in {\mathbb{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\in FM$. Such paths are called horizontal, and have zero acceleration in the sense $\mathcal{C}\left({\dot{u}}_{i,t}\right)=0$. For each $x\in 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\oplus VFM$). The split induces a map ${\pi }_{*}:HFM\to TM$, see Figure 3. For fixed $u\in FM$, the restriction ${\pi }_{*}{|}_{H_{u}FM}:H_{u}FM\to 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\in {\mathbb{R}}^d$ by $H_e\left(u\right)=h_u\left(ue\right)$. In particular, the standard basis $(e_1,\dots ,e_d)$ on ${\mathbb{R}}^d$ gives d globally defined horizontal vector fields $H_i\in HFM$, $i=1,\dots ,d$ by $H_i=H_{e_i}$. Intuitively, the fields $H_i\left(u\right)$ model infinitesimal transformations in M of $x_0$ in direction $u_i=u{e}_i$ with corresponding infinitesimal parallel transport of the vectors $u_1,\dots ,u_d$ of the frame along the direction $u_i$. A horizontal lift of a differentiable curve $x_t\in 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$.

2.2. 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 ${\mathbb{R}}^d$ with components $s_t^i$ such that ${\dot{x}}_t=H_i\left(u\right)s_t^i$, $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 ${\dot{x}}_t$ as defined above. Similarly, let $W_t$ be an ${\mathbb{R}}^d$ valued Brownian motion so that sample paths $W_t\left(\omega \right)\in W\left({\mathbb{R}}^d\right)$. A solution to the stochastic differential equation $d{U}_t={\sum }_{i=1}^d{H}_i\left(U_t\right)\circ d{W}_t^i$ in $FM$ is called a stochastic development of $W_t$ in $FM$. The solution projects to a stochastic development $X_t=\pi \left(U_t\right)$ in M. We call the process $W_t$ in ${\mathbb{R}}^d$, that through φ maps to $X_t$, the driving process of $X_t$. Let ${\phi }_u:W\left({\mathbb{R}}^d\right)\to W\left(M\right)$ be the map that for fixed u sends a path in ${\mathbb{R}}^d$ to its development on M. Its inverse ${\phi }_u^{-1}$ is the anti-development in ${\mathbb{R}}^d$ of paths on M given u.

Equivalent to the fact that normal distributions $\mathcal{N}(0,\Sigma )$ in ${\mathbb{R}}^d$ can be obtained as the transition distributions of diffusion processes ${\Sigma }^{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 ${\Sigma }^{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 ${\int }_{\mathrm{Diff}}:FM\to \mathrm{Prob}\left(M\right)$ to the space of probability distributions on M. For each point $u\in FM$, the map sends a Brownian motion in ${\mathbb{R}}^d$ to a distribution ${\mu }_u$ by stochastic development of the process $U_t$ in $FM$, starting at u and letting ${\mu }_u$ be the distribution of $X=\pi \left(U_1\right)$. The pair $(x,u)$, $x=\pi \left(u\right)$ is analogous to the parameters $(\mu ,\Sigma )$ for a Euclidean normal distribution: the point $x\in M$ represents the starting point of the diffusion, and the frame u represents a square root ${\Sigma }^{1/2}$ of the covariance Σ. In the general situation where ${\mu }_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 ${\int }_{\mathrm{Diff}}$ by maximum likelihood to observed samples in M; i.e., it optimizes for the most likely parameters $u=(x,u_1,\dots ,u_d)$ for the anisotropic diffusion process, giving a fit to the data of the manifold generalization of the Euclidean normal distribution.

2.3. 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,\xi )$ be a local trivialization of $T^{*}FM$, and let $(x^i,u_{\alpha }^i)$ be coordinates on $FM$ with $u_{\alpha }^i$ satisfying $u_{\alpha }=u_{\alpha }^i{\partial }_{x^i}$ for each $\alpha =1,\dots ,d$.

To find a basis that is adapted to the horizontal distribution, define the d linearly independent vector fields $D_j={\partial }_{x^j}-{\mathsf{\Gamma }}_j^{h_{\gamma }}{\partial }_{u_{\gamma }^h}$ where ${\mathsf{\Gamma }}_j^{h_{\gamma }}={{\mathsf{\Gamma }}^h}_{ij}{u}_{\gamma }^i$ is the contraction of the Christoffel symbols ${{\mathsf{\Gamma }}^h}_{ij}$ for the connection $\mathcal{C}$ with $u_{\alpha }^i$. We denote this adapted frame D. The vertical distribution is correspondingly spanned by $D_{j_{\beta }}={\partial }_{u_{\beta }^j}$. The vectors $D^h=d{x}^h,$ and $D^{h_{\gamma }}={\mathsf{\Gamma }}_j^{h_{\gamma }}d{x}^j+d{u}_{\gamma }^h$ constitutes a dual coframe $D^{*}$. The map ${\pi }_{*}:HFM\to TM$ is in coordinates of the adapted frame ${\pi }_{*}\left(w^j{D}_j\right)=w^j{\partial }_{x^j}$. Correspondingly, the horizontal lift $h_u$ is $h_u\left(w^j{\partial }_{x^j}\right)=w^j{D}_j$. The map $u:{\mathbb{R}}^d\to T_{x}M$ is given by the matrix $\left[u_{\alpha }^i\right]$ so that $uv=u_{\alpha }^i{v}^{\alpha }{\partial }_{x^i}=u_{\alpha }{v}^{\alpha }$.

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_{\alpha }^i)$ and the adapted frame D and coframe $D^{*}$. We have

$${}_{({\partial }_{x^i},{\partial }_{u_{\alpha }^i})}{A}_D=\left[\begin{array}{cc}I& 0\\ -\mathsf{\Gamma }& I\end{array}\right]withinverse{}_D{A}_{({\partial }_{x^i},{\partial }_{u_{\alpha }^i})}=\left[\begin{array}{cc}I& 0\\ \mathsf{\Gamma }& I\end{array}\right]$$

writing Γ for the matrix $\left[{\mathsf{\Gamma }}_j^{h_{\gamma }}\right]$. Similarly, the component matrices of the dual frame $D^{*}$ are

$${}_{{({\partial }_{x^i},{\partial }_{u_{\alpha }^i})}^{*}}{A}_{D^{*}}=\left[\begin{array}{cc}I& {\mathsf{\Gamma }}^T\\ 0& I\end{array}\right]and{}_{D^{*}}{A}_{{({\partial }_{x^i},{\partial }_{u_{\alpha }^i})}^{*}}=\left[\begin{array}{cc}I& -{\mathsf{\Gamma }}^T\\ 0& I\end{array}\right].$$

2.4. Connection and Curvature

The $TM$ valued connection $\mathcal{C}:TM\times TM\to TM$ lifts to a principal connection $TFM\times TFM\to VFM$ on the principal bundle $FM$. $\mathcal{C}$ can then be identified with the $\mathfrak{gl}\left(n\right)$-valued connection form ω on $TFM$. The identification occurs by the isomorphism ψ between $FM\times \mathfrak{gl}\left(n\right)$ and $VFM$ given by $\psi (u,v)=\frac{d}{dt}{uexp\left(tv\right)|}_{t=0}$ (e.g., [19,20]).

The map ψ is equivariant with respect to the $\mathrm{GL}\left(n\right)$ action $g\mapsto u{g}^{-1}$ on $FM$. In order to explicitly see the connection between the usual covariant derivative $\nabla :\mathsf{\Gamma }\left(TM\right)\times \mathsf{\Gamma }\left(TM\right)\to \mathsf{\Gamma }\left(TM\right)$ on M determined by $\mathcal{C}$ and $\mathcal{C}$ regarded as a connection on the principal bundle $FM$, following [19], we let $s:M\to TM$ be a local vector field on M; equivalently, $s\in \mathsf{\Gamma }\left(TM\right)$ is a local section of $TM$. s determines a map $s^{FM}:FM\to {\mathbb{R}}^d$ by $s^{FM}\left(u\right)=u^{-1}s\left(\pi \left(u\right)\right)$; i.e., it gives the coordinates of $s\left(x\right)$ in the frame u at x. The pushforward ${\left(s^{FM}\right)}_{*}:TFM\to {\mathbb{R}}^d$ has in its ith component the exterior derivative $d{\left(s^{FM}\right)}^i$. Let now $w\left(x\right)$ be a local section of $FM$. The composition $w\circ {\left(s^{FM}\right)}_{*}\circ h_w:TM\to TM$ is identical to the covariant derivative ${\nabla }_{·}s:TM\to TM$. The construction is independent of the choice of w because of the $\mathrm{GL}\left(n\right)$-equivariance of $s^{FM}$. The connection form ω can be expressed as the matrix $(s_1^{FM}\circ h_w,\dots ,s_d^{FM}\circ h_w)$ when letting $s_i^{FM}\left(u\right)=e_i$.

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 $s_t=w_{t,i}{s}_t^i$. Then, $s^{FM}\left(w_t\right)={(s_t^1,\dots ,s_t^d)}^T$, and

$$w_t^{-1}{\nabla }_{{\dot{x}}_t}s={\left(s^{F}M\right)}_{*}\left(h_{w_t}\left({\dot{x}}_t\right)\right)=\frac{d}{dt}{(s_t^1,\dots ,s_t^d)}^T;$$

(1)

i.e., the covariant derivative takes the form of the standard derivative applied to the frame coordinates $s_t^i$.

The curvature tensor $R\in {\mathcal{T}}_1^3\left(M\right)$ gives the $\mathfrak{gl}\left(n\right)$-valued curvature form $\Omega :TFM\times TFM\to \mathfrak{gl}\left(n\right)$ on $TFM$ by

$$\Omega (v_u,w_u)=u^{-1}R({\pi }_{*}\left(v_u\right),{\pi }_{*}\left(w_u\right))u,v_u,w_v\in TFM.$$

Note that $\Omega (v_u,w_u)=\Omega (h_u\left({\pi }_{*}\left(v_u\right)\right),h_u\left({\pi }_{*}\left(w_u\right)\right))$, which we can use to write the curvature R as the $\mathfrak{gl}\left(n\right)$-valued map $R_u:T^2\left(T_{\pi \left(u\right)}M\right)\to \mathfrak{gl}\left(n\right)$, $(v,w)\mapsto \Omega (h_u\left({\pi }_{*}\left(v_u\right)\right),h_u\left({\pi }_{*}\left(w_u\right)\right))$ for fixed u. In coordinates, the curvature is

$$R{}_{ijk}^s=\mathsf{\Gamma }{}_{ik}^l\mathsf{\Gamma }{}_{jl}^s-\mathsf{\Gamma }{}_{jk}^l\mathsf{\Gamma }{}_{il}^s+\mathsf{\Gamma }{}_{ik;j}^s-\mathsf{\Gamma }{}_{jk;i}^s$$

where ${{\mathsf{\Gamma }}^s}_{ik;j}={\partial }_{x^j}{{\mathsf{\Gamma }}^s}_{ik}$.

Let $x_{t,s}$ be a family of paths in M, and let $u_{t,s}\in {\pi }^{-1}\left(x_{t,s}\right)$ be horizontal lifts of $x_{t,s}$ for each fixed s. Write ${\dot{x}}_{t,s}={\partial }_t{x}_{t,s}$ and ${\dot{u}}_{t,s}={\partial }_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$

$$\begin{array}{cc}\hfill {\partial }_s{u}_{t,s}& =\psi \left(u_{t,s},{\psi }_{u_{0,s}}^{-1}\left(\mathcal{C}\left({\partial }_s{u}_{0,s}\right)\right)+{\int }_0^s\Omega ({\dot{u}}_{r,s},{\partial }_s{u}_{r,s})dr\right)+h_{u_{t,s}}\left({\partial }_s{x}_{t,s}\right)\hfill \\ & =\psi \left(u_{t,s},{\psi }_{0,s}^{-1}\left(\mathcal{C}\left({\partial }_s{u}_{0,s}\right)\right)+{\int }_0^s{R}_{u_{r,s}}({\dot{x}}_{r,s},{\partial }_s{x}_{r,s})dr\right)+h_{u_{t,s}}\left({\partial }_s{x}_{t,s}\right).\hfill \end{array}$$

(2)

This follows from the structure equation $d\omega =-\omega \wedge \omega +\Omega$ (e.g., [21]). Note that the curve depends on the vertical variation $\mathcal{C}\left({\partial }_s{u}_{0,s}\right)$ at only one point along the curve. The remaining terms depend on the horizontal variation or, equivalently, ${\partial }_s{x}_{t,s}$. The t-derivative of ${\partial }_s{u}_{t,s}$ is the curve in $TTFM$ satisfying

$$\begin{array}{cc}\hfill {\partial }_s{h}_{u_{t,s}}\left({\dot{x}}_{t,s}\right)& =\psi \left(u_{t,s},R_{u_{t,s}}({\dot{x}}_{t,s},{\partial }_s{x}_{t,s})\right)+{\partial }_t\psi \left(u_{t,s},{\psi }_{0,s}^{-1}\left(\mathcal{C}\left({\partial }_s{u}_{0,s}\right)\right)\right)+{\partial }_t\left(h_{u_{t,s}}\left({\partial }_s{x}_{t,s}\right)\right)\hfill \\ & =\psi \left(u_{t,s},R_{u_{t,s}}({\dot{x}}_{t,s},{\partial }_s{x}_{t,s})\right)+{\partial }_t\psi \left(u_{t,s},{\psi }_{0,s}^{-1}\left(\mathcal{C}\left({\partial }_s{u}_{0,s}\right)\right)\right)\hfill \\ & +h_{u_{t,s}}\left({\partial }_t{\partial }_s{x}_{t,s}\right)+\left({\partial }_t{h}_{u_{t,s}}\right)\left({\partial }_s{x}_{t,s}\right).\hfill \end{array}$$

(3)

Here, the first and third term in the last expression are identified with elements of $T_{{\partial }_s{u}_{t,s}}TFM$ by the natural mapping $T_{u_{t,s}}FM\to T_{{\partial }_s{u}_{t,s}}TFM$. When $\mathcal{C}\left({\partial }_s{u}_{0,s}\right)$ 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$.

3. 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

$$ln{\tilde{p}}_{\Sigma }\left(x_t\right)\propto -{\int }_0^1{\parallel {\dot{x}}_t\parallel }_{\Sigma }^{2}dt+c_{\Sigma }$$

(4)

with the norm ${\parallel ·\parallel }_{\Sigma }$ given by the inner product ${〈v,w〉}_{\Sigma }=〈{\Sigma }^{-1/2}v,{\Sigma }^{-1/2}w〉=v{\Sigma }^{-1}w$; i.e., the inner product weighted by the precision matrix ${\Sigma }^{-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 ${\mathbb{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 ${\mathbb{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},\dots ,u_{d,t})$ such that the matrix $\left[u_{\alpha ,t}^i\right]$ represents the square root covariance matrix ${\Sigma }^{1/2}$ at $x_t$. Since $u_t$ being a frame defines an invertible map ${\mathbb{R}}^d\to T_{x_t}M$, the norm ${\parallel ·\parallel }_{\Sigma }$ above has a direct analogue in the norm ${\parallel ·\parallel }_{u_t}$ defined by the inner product

$${〈v,w〉}_{u_t}={〈u_t^{-1}v,u_t^{-1}w〉}_{{\mathbb{R}}^d}$$

(5)

for vectors $v,w\in 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\to T_{x}M$. Again using that u can be considered a map ${\mathbb{R}}^d\to T_{x}M$, $g_u$ is defined by $\xi \mapsto u{(\xi \circ u)}^{♯}$, where ♯ is the standard identification ${\left({\mathbb{R}}^d\right)}^{*}\to {\mathbb{R}}^d$. The sequence of mappings defining $g_u$ is illustrated below:

$$\begin{array}{ccccccc}{T}_x^{*}M& \to & {\left({\mathbb{R}}^d\right)}^{*}& \to & {\mathbb{R}}^d& \to & T_{x}M\\ \xi & \mapsto & \xi \circ u& \mapsto & {(\xi \circ u)}^{♯}& \mapsto & u{(\xi \circ u)}^{♯}.\end{array}$$

(6)

This definition uses the ${\mathbb{R}}^d$ inner product in the definition of ♯. Its inverse gives the cometric $g_u^{-1}:T_{x}M\to T_x^{*}M$; i.e., $v\mapsto {\left(u^{-1}v\right)}^{♭}\circ u^{-1}$.

$$\begin{array}{ccccccc}{T}_{x}M& \to & {\mathbb{R}}^d& \to & {\left({\mathbb{R}}^d\right)}^{*}& \to & T_x^{*}M\\ v& \mapsto & u^{-1}v& \mapsto & {\left(u^{-1}\right)}^{♭}& \mapsto & {\left(u^{-1}\right)}^{♭}\circ u^{-1}.\end{array}$$

(7)

3.1. 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,\tilde{w}\in H_{u}FM$, the lift of (5) by ${\pi }_{*}$ is the inner product

$$〈w,\tilde{w}〉={〈u^{-1}{\pi }_{*}w,u^{-1}{\pi }_{*}\tilde{w}〉}_{{\mathbb{R}}^d}.$$

The inner product induces a sub-Riemannian metric $g_{FM}:TF{M}^{*}\to HFM\subset TFM$ by

$$〈w,g_{FM}\left(\xi \right)〉=\left(\xi \right|w),\forall w\in H_{u}FM$$

(8)

with $\left(\xi \right|w)$ denoting the evaluation $\xi \left(w\right)$ for the covector $\xi \in 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

$$\xi \mapsto h_u\left(g_u(\xi \circ h_u)\right),\xi \in T_{u}FM$$

(9)

of the metric $g_u$ above. In frame coordinates, the metric takes the form

$$u^{-1}{\pi }_{*}{g}_{FM}\left(\xi \right)=\left(\begin{array}{c}\xi \left(H_1\left(u\right)\right)\\ ⋮\\ \xi \left(H_d\left(u\right)\right)\end{array}\right).$$

(10)

In terms of the adapted coordinates for $TFM$ described in Section 2.3, with $w=w^j{D}_j$ and $\tilde{w}={\tilde{w}}^j{D}_j$, we have

$$\begin{array}{cc}\hfill 〈w,\tilde{w}〉& =〈w^i{D}_i,{\tilde{w}}^j{D}_j〉=〈u^{-1}{w}^i{\partial }_{x^i},u^{-1}{\tilde{w}}^j{\partial }_{x^j}〉\hfill \\ & ={〈w^i{u}_i^{\alpha },{\tilde{w}}^j{u}_j^{\alpha }〉}_{{\mathbb{R}}^d}={\delta }_{\alpha \beta }{w}^i{u}_i^{\alpha }{\tilde{w}}^j{u}_j^{\beta }=W_{ij}{w}^i{\tilde{w}}^j\hfill \end{array}$$

where $\left[u_i^{\alpha }\right]$ is the inverse of $\left[u_{\alpha }^i\right]$ and $W_{ij}={\delta }_{\alpha \beta }{u}_i^{\alpha }{u}_j^{\beta }$. Define now $W^{kl}={\delta }^{\alpha \beta }{u}_{\alpha }^k{u}_{\beta }^l$, so that $W^{ir}{W}_{rj}={\delta }_j^i$ and $W_{ir}{W}^{rj}={\delta }_i^j$. We can then write the metric $g_{FM}$ directly as

$$g_{FM}({\xi }_h{D}^h+{\xi }_{h_{\gamma }}{D}^{h_{\gamma }})=W^{ih}{\xi }_h{D}_i,$$

(11)

because $〈w,g_{FM}\left(\xi \right)〉=〈w,W^{jh}{\xi }_h{D}_j〉=W_{ij}{w}^i{W}^{jh}{\xi }_h=w^i{\xi }_i={\xi }_h{D}^h\left(w^j{D}_j\right)=\xi \left(w\right)$. 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

$$l_{FM}\left(u_t\right)=\int g_{FM}({\dot{u}}_t,{\dot{u}}_t)dt;$$

i.e., the line element is $d{s}^2=W_{ij}{D}^i{D}^j$ in adapted coordinates. The corresponding distance is given by

$$d_{FM}(u_1,u_2)=inf\{l_{FM}\left(\gamma \right)\mid \gamma \left(0\right)=u_1,\gamma \left(1\right)=u_2\}.$$

If we wish to express $g_{FM}$ in canonical coordinates on $T^{*}FM$, we can switch between the adapted frame and the coordinates $(x^i,u_{\alpha }^i,{\xi }^i,{\xi }_{\alpha }^i)$. From (11), $g_{FM}$ has $D,D^{*}$ components

$${}_D{g}_{FM,D^{*}}=\left[\begin{array}{cc}{W}^{-1}& 0\\ 0& 0\end{array}\right].$$

Therefore, $g_{FM}$ has the following components in the coordinates $(x^i,u_{\alpha }^i,{\xi }_h,{\xi }_{h_{\gamma }})$

$${}_{({\partial }_{x^i},{\partial }_{u_{\alpha }^i})}{g}_{FM,{({\partial }_x,{\partial }_{u_{\alpha }^i})}^{*}}={}_{({\partial }_{x^i},{\partial }_{u_{\alpha }^i})}{A}_D{}_D{g}_{FM,D^{*}}{}_{D^{*}}{A}_{{({\partial }_{x^i},{\partial }_{u_{\alpha }^i})}^{*}}=\left[\begin{array}{cc}{W}^{-1}& -W^{-1}{\mathsf{\Gamma }}^T\\ -\mathsf{\Gamma }{W}^{-1}& \mathsf{\Gamma }{W}^{-1}{\mathsf{\Gamma }}^T\end{array}\right]$$

or $g_{FM}^{ij}=W^{ij}$, $g_{FM}^{i{j}_{\beta }}=-W^{ih}{\mathsf{\Gamma }}_h^{j_{\beta }}$, $g_{FM}^{i_{\alpha }j}=-{\mathsf{\Gamma }}_h^{i_{\alpha }}{W}^{hj}$, and $g_{FM}^{i_{\alpha }{j}_{\beta }}={\mathsf{\Gamma }}_k^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_h^{j_{\beta }}$.

3.2. 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 $\mathcal{C}\left({\dot{u}}_t\right)=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 $\mathcal{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 $\mathcal{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.

3.3. 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 $\mathfrak{gl}\left(n\right)$ to each fiber using the isomorphism $\psi {(u,·)}^{-1}:V_{u}FM\to \mathfrak{gl}\left(n\right)$. 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({\pi }_{*}\left(v_u\right),{\pi }_{*}\left(w_u\right))$ (i.e., the metric is not anisotropically weighted). The line element is in this case $d{s}^2=g_{ij}d{x}^{i}d{x}^j+X_{\beta \alpha }{g}_{ij}{D}^{{\alpha }_i}{D}^{{\beta }_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$.

4. 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=\pi \left(u_0\right)\in M$ to $y\in M$ is a smooth path $x_t:[0,1]\to M$ with $x_0=x$ and $x_1=y$ such that its anti-development ${\phi }_{u_0}^{-1}\left(x_t\right)$ is a most probable path for $W_t$; i.e.,

$$x_t\in {\mathrm{argmin}}_{\sigma ,{\sigma }_0=x,{\sigma }_1=y}{\int }_0^1-L_{{\mathbb{R}}^d}({\phi }_{u_0}^{-1}\left({\sigma }_t\right),\frac{d}{dt}{\phi }_{u_0}^{-1}\left({\sigma }_t\right))dt$$

with $L_{{\mathbb{R}}^d}$ being the Onsager–Machlup function for the process $W_t$ on ${\mathbb{R}}^d$ [22].

The definition uses the one-to-one relation between $W\left({\mathbb{R}}^d\right)$ and $W\left(M\right)$ provided by ${\phi }_{u_0}$ to characterize the paths using the ${\mathbb{R}}^d$ Onsager–Machlup function $L_{{\mathbb{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,{\dot{x}}_t)=-\frac{1}{2}{\parallel {\dot{x}}_t\parallel }_{g_R}^2+\frac{1}{12}{S}_{g_R}\left(x_t\right)$ with $S_{g_R}$ denoting the scalar curvature. This curvature term vanishes on ${\mathbb{R}}^d$, and therefore $L_{{\mathbb{R}}^d}({\gamma }_t,{\dot{\gamma }}_t)=-\frac{1}{2}{\parallel {\dot{\gamma }}_t\parallel }^2$ for a curve ${\gamma }_t\in {\mathbb{R}}^d$.

By pulling $x_t\in M$ back to ${\mathbb{R}}^d$ using ${\phi }_{u_0}^{-1}$, the construction removes the $\frac{1}{12}{S}_{g_R}\left(x_t\right)$ 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\left(u_0\right)$ denote the principal sub-bundle of $FM$ of points $z\in FM$ reachable from $u_0\in FM$ by horizontal paths. Suppose the Hörmander condition is satisfied on $Q\left(u_0\right)$, and that $Q\left(u_0\right)$ has compact fibers. Then, most probable paths from $x_0$ to $y\in 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 ${\pi }^{-1}\left(y\right)$.

Below, we will derive evolution equations for the set of such extremal paths that correspond to normal sub-Riemannian geodesics.

4.1. Normal Geodesics for $g_{FM}$

Connected to the metric $g_{FM}$ is the Hamiltonian

$$H\left(z\right)=\frac{1}{2}\left(z\right|g_{FM}\left(z\right))$$

(12)

on the symplectic space $T^{*}FM$. Letting $\widehat{\pi }$ denote the projection on the bundle $T^{*}FM\to FM$, (8) gives

$$H\left(z\right)=\frac{1}{2}〈g_{FM}\left(z\right)|g_{FM}\left(z\right)〉=\frac{1}{2}{\parallel z\circ h_{\widehat{\pi }\left(z\right)}\circ \widehat{\pi }\left(z\right)\parallel }_{{\left({\mathbb{R}}^d\right)}^{*}}^2=\frac{1}{2}\sum _{i=1}^d\xi {\left(H_i\left(u\right)\right)}^2.$$

Normal geodesics in sub-Riemannian manifolds satisfy the Hamilton–Jacobi equations [23] with Hamiltonian flow

$${\dot{z}}_t=X_H={\Omega }^{\#}dH\left(z\right)$$

(13)

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={\pi }_{T^{*}FM}\left(z_t\right)$ 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,\xi )$, (13) gives the Hamilton–Jacobi equations

$$\begin{array}{cc}\hfill \dot{u}& ={\partial }_{\xi }H(u,\xi )=g_{FM}(u,\xi )=h_u\left(u{(\xi \left(H_1\left(u\right)\right),\dots ,\xi \left(H_d\left(u\right)\right))}^T\right)\hfill \\ \hfill \dot{\xi }& =-{\partial }_{u}H(u,\xi )=-{\partial }_u\frac{1}{2}{\parallel \xi \circ h_u\circ u\parallel }_{{\left({\mathbb{R}}^d\right)}^{*}}^2=-{\partial }_u\frac{1}{2}\sum _{i=1}^d\xi {\left(H_i\left(u\right)\right)}^2.\hfill \end{array}$$

(14)

Using (3), we have for the second equation

$$\begin{array}{cc}\hfill \dot{\xi }& =-\sum _{i=1}^d\xi \left(H_i\left(u\right)\right)\xi \left({\partial }_u{h}_u\left(u{e}_i\right)\right)\hfill \\ & =-\sum _{i=1}^d\xi \left(H_i\left(u\right)\right)\xi \left(\psi (u,R_u(u{e}_i,{\pi }_{*}\left({\partial }_u\right)))+{\partial }_{h_u\left(u{e}_i\right)}\psi \left(u,{\psi }^{-1}\left(\mathcal{C}\left({\partial }_u\right)\right)\right)+{\partial }_{h_u\left(u{e}_i\right)}{h}_u\left({\pi }_{*}\left({\partial }_u\right)\right)\right)\hfill \\ & =-\xi \left(\psi (u,R_u({\pi }_{*}\left(\dot{u}\right),{\pi }_{*}\left({\partial }_u\right)))+{\partial }_{\dot{u}}\psi \left(u,{\psi }^{-1}\left(\mathcal{C}\left({\partial }_u\right)\right)\right)+{\partial }_{\dot{u}}{h}_u\left({\pi }_{*}\left({\partial }_u\right)\right)\right).\hfill \end{array}$$

(15)

Here ${\partial }_{\dot{u}}$ denotes u-derivative in the direction $\dot{u}$, equivalently ${\partial }_{\dot{u}}{h}_u\left(v\right)={\partial }_t\left(h_u\right)\left(v\right)$. 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 $\psi (u,R_u({\pi }_{*}\left(\dot{u}\right),{\pi }_{*}\left({\partial }_u\right))$. 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, $\mathcal{C}$ the $g_R$ Levi–Civita connection, and $u_0$ is $g_R$ orthonormal, $g_{FM}(h_u\left(v\right),h_u\left(w\right))=g_R(v,w)$ for all $v,w\in T_{\pi \left(u_t\right)}M$. In this case, a normal MPP $\pi \left(u_t\right)$ will be a Riemannian $g_R$ geodesic.

4.2. Evolution in Coordinates

In coordinates $u=(x^i,u_{\alpha }^i,{\xi }_i,{\xi }_{i_{\alpha }})$ for $T^{*}FM$, we can equivalently write

$$\begin{array}{cc}& {\dot{x}}^i=g^{ij}{\xi }_j+g^{i{j}_{\beta }}{\xi }_{j_{\beta }}=W^{ij}{\xi }_j-W^{ih}{\mathsf{\Gamma }}_h^{j_{\beta }}{\xi }_{j_{\beta }}\hfill \\ & {\dot{X}}_{\alpha }^i=g^{i_{\alpha }j}{\xi }_j+g^{i_{\alpha }{j}_{\beta }}{\xi }_{j_{\beta }}=-{\mathsf{\Gamma }}_h^{i_{\alpha }}{W}^{hj}{\xi }_j+{\mathsf{\Gamma }}_k^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_h^{j_{\beta }}{\xi }_{j_{\beta }}\hfill \\ & {\dot{\xi }}_i=-\frac{1}{2}\left({\partial }_{y^i}{g}_y^{hk}{\xi }_h{\xi }_k+{\partial }_{y^i}{g}_y^{h{k}_{\delta }}{\xi }_h{\xi }_{k_{\delta }}+{\partial }_{y^i}{g}_y^{h_{\gamma }k}{\xi }_{h_{\gamma }}{\xi }_k+{\partial }_{y^i}{g}_y^{h_{\gamma }{k}_{\delta }}{\xi }_{h_{\gamma }}{\xi }_{k_{\delta }}\right)\hfill \\ & {\dot{\xi }}_{i_{\alpha }}=-\frac{1}{2}\left({\partial }_{y^{i_{\alpha }}}{g}_y^{hk}{\xi }_h{\xi }_k+{\partial }_{y^{i_{\alpha }}}{g}_y^{h{k}_{\delta }}{\xi }_h{\xi }_{k_{\delta }}+{\partial }_{y^{i_{\alpha }}}{g}_y^{h_{\gamma }k}{\xi }_{h_{\gamma }}{\xi }_k+{\partial }_{y^{i_{\alpha }}}{g}_y^{h_{\gamma }{k}_{\delta }}{\xi }_{h_{\gamma }}{\xi }_{k_{\delta }}\right)\hfill \end{array}$$

with ${\mathsf{\Gamma }}_{k,i}^{h_{\gamma }}$ for ${\partial }_{y^i}{\mathsf{\Gamma }}_k^{h_{\gamma }}$, and where

$$\begin{array}{cc}& {\partial }_{y^l}{g}^{ij}=0,{\partial }_{y^l}{g}^{i{j}_{\beta }}=-W^{ih}{\mathsf{\Gamma }}_{h,l}^{j_{\beta }},{\partial }_{y^l}{g}^{i_{\alpha }j}=-{\mathsf{\Gamma }}_{h,l}^{i_{\alpha }}{W}^{hj},{\partial }_{y^l}{g}^{i_{\alpha }{j}_{\beta }}={\mathsf{\Gamma }}_{k,l}^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_h^{j_{\beta }}+{\mathsf{\Gamma }}_k^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_{h,l}^{j_{\beta }},\hfill \end{array}$$

$$\begin{array}{cc}& {\partial }_{y^{l_{\zeta }}}{g}^{ij}={W^{ij}}_{,l\zeta },{\partial }_{y^{l_{\zeta }}}{g}^{i{j}_{\beta }}=-{W^{ih}}_{,l\zeta }{\mathsf{\Gamma }}_h^{j_{\beta }}-W^{ih}{\mathsf{\Gamma }}_{h,l_{\zeta }}^{j_{\beta }},{\partial }_{y^{l_{\zeta }}}{g}^{i_{\alpha }j}=-{\mathsf{\Gamma }}_{h,l_{\zeta }}^{i_{\alpha }}{W}^{hj}-{\mathsf{\Gamma }}_h^{i_{\alpha }}{W^{hj}}_{,l\zeta },\hfill \\ & {\partial }_{y^{l_{\zeta }}}{g}^{i_{\alpha }{j}_{\beta }}={\mathsf{\Gamma }}_{k,l_{\zeta }}^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_h^{j_{\beta }}+{\mathsf{\Gamma }}_k^{i_{\alpha }}{W^{kh}}_{,l\zeta }{\mathsf{\Gamma }}_h^{j_{\beta }}+{\mathsf{\Gamma }}_k^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_{h,l_{\zeta }}^{j_{\beta }},\hfill \\ & {\mathsf{\Gamma }}_{h,l_{\zeta }}^{i_{\alpha }}={\partial }_{y^{l_{\zeta }}}\left(\mathsf{\Gamma }{}_{hk}^i{u}_{\alpha }^k\right)={\delta }^{\zeta \alpha }\mathsf{\Gamma }{}_{hl}^i,{W^{ij}}_{,l\zeta }={\delta }^{il}{u}_{\zeta }^j+{\delta }^{jl}{u}_{\zeta }^i.\hfill \end{array}$$

Combining these expressions, we obtain

$$\begin{array}{cc}& {\dot{x}}^i=W^{ij}{\xi }_j-W^{ih}{\mathsf{\Gamma }}_h^{j_{\beta }}{\xi }_{j_{\beta }},{\dot{X}}_{\alpha }^i=-{\mathsf{\Gamma }}_h^{i_{\alpha }}{W}^{hj}{\xi }_j+{\mathsf{\Gamma }}_k^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_h^{j_{\beta }}{\xi }_{j_{\beta }}\hfill \\ & {\dot{\xi }}_i=W^{hl}{\mathsf{\Gamma }}_{l,i}^{k_{\delta }}{\xi }_h{\xi }_{k_{\delta }}-\frac{1}{2}\left({\mathsf{\Gamma }}_{k,i}^{h_{\gamma }}{W}^{kh}{\mathsf{\Gamma }}_h^{k_{\delta }}+{\mathsf{\Gamma }}_k^{h_{\gamma }}{W}^{kh}{\mathsf{\Gamma }}_{h,i}^{k_{\delta }}\right){\xi }_{h_{\gamma }}{\xi }_{k_{\delta }}\hfill \\ & {\dot{\xi }}_{i_{\alpha }}={\mathsf{\Gamma }}_{k,i_{\alpha }}^{h_{\delta }}{W}^{kh}{\mathsf{\Gamma }}_h^{k_{\delta }}{\xi }_{h_{\gamma }}{\xi }_{k_{\delta }}-\left({W^{hl}}_{,i\alpha }{\mathsf{\Gamma }}_l^{k_{\delta }}+W^{hl}{\mathsf{\Gamma }}_{l,i_{\alpha }}^{k_{\delta }}\right){\xi }_h{\xi }_{k_{\delta }}-\frac{1}{2}\left({W^{hk}}_{,i\alpha }{\xi }_h{\xi }_k+{\mathsf{\Gamma }}_k^{h_{\delta }}{W^{kh}}_{,i\alpha }{\mathsf{\Gamma }}_h^{k_{\delta }}{\xi }_{h_{\gamma }}{\xi }_{k_{\delta }}\right).\hfill \end{array}$$

4.3. Acceleration and Polynomials for $\mathcal{C}$

We can identify the covariant acceleration ${\nabla }_{{\dot{x}}_t}{\dot{x}}_t$ of curves satisfying the MPP equations, and hence normal MPPs through their frame coordinates. Let $(u_t,{\xi }_t)$ satisfy (13). Then, $u_t$ is a horizontal lift of $x_t=\pi \left(u_t\right)$ and hence by (1), (3), (10), and (15),

$$\begin{array}{cc}\hfill u_t^{-1}{\nabla }_{{\dot{x}}_t}{\dot{x}}_t& =\frac{d}{dt}\left(\begin{array}{c}\xi \left(h_{u_t}\left(u_t{e}_1\right)\right)\\ ⋮\\ \xi \left(h_{u_t}\left(u_t{e}_d\right)\right)\end{array}\right)=\left(\begin{array}{c}\dot{\xi }\left(h_{u_t}\left(u_t{e}_1\right)\right)\\ ⋮\\ \dot{\xi }\left(h_{u_t}\left(u_t{e}_d\right)\right)\end{array}\right)+\left(\begin{array}{c}\xi \left({\partial }_t{h}_{u_t}\left(u_t{e}_1\right)\right)\\ ⋮\\ \xi \left({\partial }_t{h}_{u_t}\left(u_t{e}_d\right)\right)\end{array}\right)\hfill \\ & =-\left(\begin{array}{c}\xi ({\partial }_{h_{u_t}\left(u_t{e}_1\right)}{h}_{u_t}\left({\pi }_{*}\left({\dot{u}}_t\right)\right)\\ ⋮\\ \xi ({\partial }_{h_{u_t}\left(u_t{e}_d\right)}{h}_{u_t}\left({\pi }_{*}\left({\dot{u}}_t\right)\right)\end{array}\right)+\left(\begin{array}{c}\xi \left({\partial }_{h_{u_t}\left({\pi }_{*}\left({\dot{u}}_t\right)\right)}{h}_{u_t}\left(u_t{e}_1\right)\right)\\ ⋮\\ \xi \left({\partial }_{h_{u_t}\left({\pi }_{*}\left({\dot{u}}_t\right)\right)}{h}_{u_t}\left(u_t{e}_d\right)\right)\end{array}\right)\hfill \\ & =\left(\begin{array}{c}\xi \left(\psi (u_t,R_{u_t}(u_t{e}_1,{\pi }_{*}\left({\dot{u}}_t\right)))\right)\\ ⋮\\ \xi \left(\psi (u_t,R_{u_t}(u_t{e}_d,{\pi }_{*}\left({\dot{u}}_t\right)))\right)\end{array}\right).\hfill \end{array}$$

(16)

The fact that the covariant derivative vanishes for classical geodesic leads to a definition of higher-order polynomials through the covariant derivative by requiring ${\left({\nabla }_{{\dot{x}}_t}\right)}^k{\dot{x}}_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 $\mathcal{C}$ while still satisfying the horizontality constraint on $FM$. This makes it natural to ask if normal MPPs relate to polynomials defined using $\mathcal{C}$. For curves satisfying the MPP equations, using (16) and (15), we have

$$\begin{array}{cc}\hfill u_t^{-1}{\left({\nabla }_{{\dot{x}}_t}\right)}^2{\dot{x}}_t& =\frac{d}{dt}\left(\begin{array}{c}\xi \left(\psi (u_t,R_{u_t}(u_t{e}_1,{\pi }_{*}\left({\dot{u}}_t\right)))\right)\\ ⋮\\ \xi \left(\psi (u_t,R_{u_t}(u_t{e}_d,{\pi }_{*}\left({\dot{u}}_t\right)))\right)\end{array}\right)=\left(\begin{array}{c}\xi \left(\psi (u_t,\frac{d}{dt}{R}_{u_t}(u_t{e}_1,{\pi }_{*}\left({\dot{u}}_t\right)))\right)\\ ⋮\\ \xi \left(\psi (u_t,\frac{d}{dt}{R}_{u_t}(u_t{e}_d,{\pi }_{*}\left({\dot{u}}_t\right)))\right)\end{array}\right).\hfill \end{array}$$

Thus, in general, normal MPPs are not second order polynomials in the sense ${\left({\nabla }_{{\dot{x}}_t}\right)}^2{\dot{x}}_t=0$ unless the curvature $R_{u_t}(u_t{e}_i,{\pi }_{*}\left({\dot{u}}_t\right))$ 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

$${\left({\nabla }_{{\dot{x}}_t}\right)}^2{\dot{x}}_t=-R({\nabla }_{{\dot{x}}_t}{\dot{x}}_t,x_t,){\dot{x}}_t.$$

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.

5. Cometric Formulation and Low-Rank Generator

We now investigate a cometric $g_{F^{k}M}+\lambda 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 $\lambda >0$ a weight. We assume that $g_{F^{k}M}$ is chosen so that $g_{F^{k}M}+\lambda 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\times d$ matrix for the frame, but a potentially much lower dimensional $d\times 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}+\lambda g_R$ represents a sum of covariance matrices instead of a sum of precision matrices. Thus, $g_{F^{k}M}+\lambda 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 ${\mathbb{R}}^k\to T_{x}M$. We define a cometric by

$$〈\xi ,\tilde{\xi }〉={\delta }^{\alpha \beta }\left(\xi \right|h_u\left(u_{\alpha }\right))\left(\tilde{\xi }\right|h_u\left(u_{\beta }\right))+\lambda {〈\xi ,\tilde{\xi }〉}_{g_R}$$

for $\xi ,\tilde{\xi }\in T_u^{*}{F}^{k}M$. The sum over $\alpha ,\beta$ is for $\alpha ,\beta =1,\dots ,k$. The first term is equivalent to the lift (9) of the cometric $〈\xi ,\tilde{\xi }〉=\left(\xi |g_u\left(\widehat{\xi }\right)\right)$ given $u:{\mathbb{R}}^k\to 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_{\alpha }^i)$, $\alpha =1,\dots ,k$ be a coordinate system on $F^{k}M$. The vertical distribution is in this case spanned by the $dk$ vector fields $D_{j_{\beta }}={\partial }_{u_{\beta }^j}$. Except for index sums being over k instead of d terms, the situation is thus similar to the full-rank case. Note that $\left(\xi \right|{\pi }_{*}^{-1}w)=\left(\xi \right|w^j{D}_j)=w^i{\xi }_i$. The cometric in coordinates is

$$〈\xi ,\tilde{\xi }〉={\delta }^{\alpha \beta }{u}_{\alpha }^i{\xi }_i{u}_{\beta }^j{\tilde{\xi }}_j+\lambda g_R^{ij}{\xi }_i{\tilde{\xi }}_j={\xi }_i\left({\delta }^{\alpha \beta }{u}_{\alpha }^i{u}_{\beta }^j+\lambda g_R^{ij}\right){\tilde{\xi }}_j={\xi }_i{W}^{ij}{\tilde{\xi }}_j$$

with $W^{ij}={\delta }^{\alpha \beta }{u}_{\alpha }^i{u}_{\beta }^j+\lambda g_R^{ij}$. We can then write the corresponding sub-Riemannian metric $g_{F^{k}M}$ in terms of the adapted frame D

$$g_{F^{k}M}({\xi }_h{D}^h+{\xi }_{h_{\gamma }}{D}^{h_{\gamma }})=W^{ih}{\xi }_h{D}_i$$

(17)

because $\left(\xi \right|g_{F^{k}M}\left(\tilde{\xi }\right))=〈\xi ,\tilde{\xi }〉={\xi }_i{W}^{ij}{\tilde{\xi }}_j$. That is, the situation is analogous to (11), except the term $\lambda g_R^{ij}$ 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}$:

$$\begin{array}{cc}& {\partial }_{y^l}{g}_{F^{k}M}^{ij}={W^{ij}}_{,l},{\partial }_{y^l}{g}_{F^{k}M}^{i{j}_{\beta }}=-{W^{ih}}_{,l}{\mathsf{\Gamma }}_h^{j_{\beta }}-W^{ih}{\mathsf{\Gamma }}_{h,l}^{j_{\beta }},{\partial }_{y^l}{g}_{F^{k}M}^{i_{\alpha }j}=-{\mathsf{\Gamma }}_{h,l}^{i_{\alpha }}{W}^{hj}-{\mathsf{\Gamma }}_h^{i_{\alpha }}{W^{hj}}_{,l},\hfill \\ & {\partial }_{y^l}{g}_{F^{k}M}^{i_{\alpha }{j}_{\beta }}={\mathsf{\Gamma }}_{k,l}^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_h^{j_{\beta }}+{\mathsf{\Gamma }}_k^{i_{\alpha }}{W^{kh}}_{,l}{\mathsf{\Gamma }}_h^{j_{\beta }}+{\mathsf{\Gamma }}_k^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_{h,l}^{j_{\beta }},\hfill \\ & {\partial }_{y^{l_{\zeta }}}{g}_{F^{k}M}^{ij}={W^{ij}}_{,l\zeta },{\partial }_{y^{l_{\zeta }}}{g}_{F^{k}M}^{i{j}_{\beta }}=-{W^{ih}}_{,l\zeta }{\mathsf{\Gamma }}_h^{j_{\beta }}-W^{ih}{\mathsf{\Gamma }}_{h,l_{\zeta }}^{j_{\beta }},{\partial }_{y^{l_{\zeta }}}{g}_{F^{k}M}^{i_{\alpha }j}=-{\mathsf{\Gamma }}_h^{i_{\alpha }}{W^{hj}}_{,l\zeta }-{\mathsf{\Gamma }}_{h,l_{\zeta }}^{i_{\alpha }}{W}^{hj},\hfill \\ & {\partial }_{y^{l_{\zeta }}}{g}_{F^{k}M}^{i_{\alpha }{j}_{\beta }}={\mathsf{\Gamma }}_{k,l_{\zeta }}^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_h^{j_{\beta }}+{\mathsf{\Gamma }}_k^{i_{\alpha }}W{}_{,l_{\zeta }}^{kh}{\mathsf{\Gamma }}_h^{j_{\beta }}+{\mathsf{\Gamma }}_k^{i_{\alpha }}{W}^{kh}{\mathsf{\Gamma }}_{h,l_{\zeta }}^{j_{\beta }},\hfill \\ & {\mathsf{\Gamma }}_{h,l_{\zeta }}^{i_{\alpha }}={\partial }_{y^{l_{\zeta }}}\left(\mathsf{\Gamma }{}_{hk}^i{u}_{\alpha }^k\right)={\delta }^{\zeta \alpha }\mathsf{\Gamma }{}_{hl}^i,{W^{ij}}_{,l}=\lambda g{}_R^{ij},l,{W^{ij}}_{,l\zeta }={\delta }^{il}{u}_{\zeta }^j+{\delta }^{jl}{u}_{\zeta }^i.\hfill \end{array}$$

Note that the introduction of the Riemannian metric $g_R$ implies that $W^{ij}$ are now dependent on the manifold coordinates $x^i$.

6. 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 ${\mathbb{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/.

6.1. Embedded Surfaces

We visualize normal MPPs and projections of curves satisfying the MPP Equation (13) on surfaces embedded in ${\mathbb{R}}^3$ in three cases: The sphere ${\mathbb{S}}^2$, on an ellipsoid, and on a hyperbolic surface. The surfaces are chosen in order to have both positive and negative curvature, and to have varying degree of symmetry. In all cases, an open subset of the surfaces are represented in a single chart by a map $F:{\mathbb{R}}^2\to {\mathbb{R}}^3$. For the sphere and ellipsoid, this gives a representation of the surface, except for the south pole. The metric and Christoffel symbols are calculated using the symbolic differentiation features of Theano. The integration are performed by a simple Euler integrator.

Figure 4, Figure 5, Figure 6 show families of curves satisfying the MPP equations in three cases: (1) With fixed starting point $x_0\in M$ and initial velocity ${\dot{x}}_0\in 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\in M$ but changing frame u above $x_0$; (3) fixed starting point $x_0\in M$ and frame u but varying $V^{*}FM$ vertical part of the initial momentum ${\xi }_0\in 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 ${\mathbb{R}}^2$ that for classical $\mathcal{C}$ geodesics would always be straight lines.

6.2. 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+{\left(2N\right)}^2)$ with N landmarks in ${\mathbb{R}}^2$.

Let $\{p_1,\dots ,p_N\}$ be landmarks in a subset $\Omega \subset {\mathbb{R}}^d$. The diffeomorphism group $\mathrm{Diff}(\Omega )$ acts on the left on landmarks with the action $\phi .\{p_1,\dots ,p_N\}=\{\phi \left(p_1\right),\dots ,\phi \left(p_N\right)\}$. In LDDMM, a Hilbert space structure is imposed on a linear subspace V of $L^2(\Omega ,{\mathbb{R}}^d)$ using a self-adjoint operator $L:V\to V^{*}\subset L^2(\Omega ,{\mathbb{R}}^d)$ and defining the inner product ${〈·,·〉}_V$ by

$${〈v,w〉}_V={〈Lv,w〉}_{L^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\subset \mathrm{Diff}(\Omega )$ by setting ${\parallel v\parallel }_{\phi }^2={\parallel v\circ {\phi }^{-1}\parallel }_V^2$; i.e., regarding ${〈·,·〉}_V$ an inner product on $T_{\mathrm{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=\left\{(p_1^1,\dots ,p_1^d,\dots ,p_N^1,\dots ,p_N^d)\right|(p_i^1,\dots ,p_i^d)\in {\mathbb{R}}^d\}$. The LDDMM metric on the landmark manifold M is directly related to the kernel K when written as a cometric $g_p(\xi ,\eta )={\sum }_{i,j=1}^N{\xi }^{i}K(p_i,p_j){\eta }^j$. Letting $i^k$ denote the index of the kth component of the ith landmark, the cometric is in coordinates $g_p^{i^k{j}^l}=K{(p_i,p_j)}_k^l$. The Christoffel symbols can be written in terms of derivatives of the cometric $g^{ij}$[31] (recall that ${\delta }_j^i=g^{ik}{g}_{kj}=g_{jk}{g}^{ki}$)

$${{\mathsf{\Gamma }}^k}_{ij}=\frac{1}{2}{g}_{ir}\left(g^{kl}{g^{rs}}_{,l}-g^{sl}{g^{rk}}_{,l}-g^{rl}{g^{ks}}_{,l}\right)g_{sj}.$$

(18)

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}=({\delta }_r^i+{\delta }_r^j){\partial }_{p_r^q}K{(p_i,p_j)}_k^l$. Using (18), derivatives of the Christoffel symbols can be computed

$$\begin{array}{cc}\hfill {{\mathsf{\Gamma }}^k}_{ij,\xi }& =\frac{1}{2}{g}_{ir,\xi }\left(g^{kl}{g^{rs}}_{,l}-g^{sl}{g^{rk}}_{,l}-g^{rl}{g^{ks}}_{,l}\right)g_{sj}+\frac{1}{2}{g}_{ir}\left(g^{kl}{g^{rs}}_{,l}-g^{sl}{g^{rk}}_{,l}-g^{rl}{g^{ks}}_{,l}\right)g_{sj,\xi }\hfill \\ & +\frac{1}{2}{g}_{ir}\left({g^{kl}}_{\xi }{g^{rs}}_{,l}+g^{kl}{g^{rs}}_{,l\xi }-{g^{sl}}_{,\xi }{g^{rk}}_{,l}-g^{sl}{g^{rk}}_{,l\xi }-{g^{rl}}_{,\xi }{g^{ks}}_{,l}-g^{rl}{g^{ks}}_{,l\xi }\right)g_{sj}.\hfill \end{array}$$

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 ${\mathbb{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.

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 $\mathcal{C}$.

7. 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\in {\mathbb{R}}^d$ that leads to an observation $x={\phi }_u\left(w_t\right){|}_{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.

7.1. 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

$$\widehat{x}={\mathrm{argmin}}_{x\in M}\sum _{i=1}^N{d}_{g_R}{\left(x,x_i\right)}^2$$

(19)

of observed data $x_1,\dots ,x_N\in 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 $\mathrm{Exp}.$ and $\mathrm{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

$$\widehat{x}={\mathrm{argmin}}_{u\in s\left(M\right)}\sum _{i=1}^N{d}_{FM}{\left(u,{\pi }^{-1}\left(x_i\right)\right)}^2$$

(20)

acts correspondingly to the Frechét mean estimator (19). Here $s\in \mathsf{\Gamma }\left(FM\right)$ is a (local) section of $FM$ that to $x\in M$ connects the known covariance represented by $s\left(x\right)\in FM$. The distances $d_{FM}\left(u,{\pi }^{-1}\left(x_i\right)\right)$, $u=s\left(x\right)$ are realized by MPPs from the mean candidate x to the fibers ${\pi }^{-1}\left(x_i\right)$. 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 $\mathcal{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

$$\widehat{u}={\mathrm{argmin}}_{u\in FM}\sum _{i=1}^N{d}_{FM}{\left(u,{\pi }^{-1}\left(x_i\right)\right)}^2-Nlog\left({\det }_{g_R}u\right)$$

(21)

is suggested. The normalizing term $-Nlog\left({\det }_{g_R}u\right)$ 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.

7.2. Priors and Low-Rank Estimation

The low-rank cometric formulation pursued in Section 5 gives a natural restriction of (21) to $u\in F^{k}M$, $1\le k\le d$. As for Euclidean PCA, most variance is often captured in the span of the first k eigenvectors with $k\ll 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\times k$ matrices. As a different approach for regularizing the estimator (21), the normalizing term $-Nlog\left({\det }_{g_R}u\right)$ 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.

7.3. 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.