A Geometric Framework for Stochastic Shape Analysis

TL;DR

This paper develops a stochastic model within the LDDMM framework to analyze shape uncertainties, providing methods to infer model parameters from landmark data through moment matching and likelihood optimization.

Contribution

It introduces a novel stochastic diffeomorphism model for shape analysis and proposes two methods for parameter inference from landmark data.

Findings

The stochastic model captures shape deformation uncertainties.

Two inference methods effectively estimate noise parameters.

Numerical tests validate the model's ability to recover spatial correlation length.

Abstract

We introduce a stochastic model of diffeomorphisms, whose action on a variety of data types descends to stochastic evolution of shapes, images and landmarks. The stochasticity is introduced in the vector field which transports the data in the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework for shape analysis and image registration. The stochasticity thereby models errors or uncertainties of the flow in following the prescribed deformation velocity. The approach is illustrated in the example of finite dimensional landmark manifolds, whose stochastic evolution is studied both via the Fokker-Planck equation and by numerical simulations. We derive two approaches for inferring parameters of the stochastic model from landmark configurations observed at discrete time points. The first of the two approaches matches moments of the Fokker-Planck equation to sample moments of the data, while the second approach employs an Expectation-Maximisation based algorithm using a Monte Carlo bridge sampling scheme to optimise the data likelihood. We derive and numerically test the ability of the two approaches to infer the spatial correlation length of the underlying noise.

Equations

$$\begin{array}{rl}d{\mathbf{q}}_i& =\frac{\partial h}{\partial {\mathbf{p}}_i}dt+\sum _l{\sigma }_l({\mathbf{q}}_i)\circ d{W}_t^l,\\ d{\mathbf{p}}_i& =-\frac{\partial h}{\partial {\mathbf{q}}_i}dt-\sum _l\frac{\partial }{\partial {\mathbf{q}}_i}\left({\mathbf{p}}_i\cdot {\sigma }_l({\mathbf{q}}_i)\right)\circ d{W}_t^l.\end{array}$$

$$d{\mathbf{q}}_i^{\alpha }$$

$$E(u_t)={\int }_0^{1}l(u_t)dt+\frac{1}{2{\lambda }^2}\Vert g_1.I_0-I_1{\Vert }^2,$$

$${\partial }_t{g}_t=u_t{g}_t,$$

$$\frac{d}{dt}\frac{\delta l}{\delta u}+{\mathrm{ad}}_{u_t}^{\ast }\frac{\delta l}{\delta u}=0,$$

$$m(t)=-\frac{1}{{\lambda }^2}{J}_t^0\diamond (g_{t,1}(J_1^0-J_1^1{)}^{♭}),$$

$$⟨I\diamond f,u{⟩}_{\mathfrak{g}\times {\mathfrak{g}}^{\ast }}:=⟨f,u.I{⟩}_{V\times V^{\ast }}.$$

$$d{g}_t=u_t{g}_{t}dt+\sum _{l=1}^J{\sigma }_l{g}_t\circ d{W}_t^l,$$

$$d\frac{\delta l}{\delta u}+{\mathrm{ad}}_{u_t}^{\ast }\frac{\delta l}{\delta u}dt+\sum _{l=1}^J{\mathrm{ad}}_{{\sigma }_l}^{\ast }\frac{\delta l}{\delta u}\circ d{W}_t^l=0.$$

$$\delta d{g}_t=\delta ugdt+u\delta gdt+\sum _{l=1}^J{\sigma }_l\delta g\circ d{W}_t^l,$$

$$d(g^{-1}\delta g)$$

$$=-g^{-1}(udt+\sum _{l=1}^J{\sigma }_l\circ d{W}_t^l)\delta g+g^{-1}(\delta ugdt+u\delta gdt+\sum _{l=1}^J{\sigma }_l\delta g\circ d{W}_t^l)$$

$$=g^{-1}\delta ugdt$$

$$:={\mathrm{Ad}}_g\delta udt.$$

$$\frac{\delta l}{\delta u}={\mathrm{Ad}}_{g^{-1}}^{\ast }(I_0\diamond (g_1^{-1}\pi )),\mathrm{where}\pi =\frac{1}{{\lambda }^2}(g_i{I}_0-I_1{)}^{♭}.$$

$$d{\mathrm{Ad}}_{g^{-1}}^{\ast }\eta =-{\mathrm{ad}}_{dg{g}^{-1}}^{\ast }{\mathrm{Ad}}_{g^{-1}}^{\ast }\eta ,$$

$$d\frac{\delta l}{\delta u}$$

$$={\mathrm{ad}}_{dg{g}^{-1}}^{\ast }{\mathrm{Ad}}_{g^{-1}}^{\ast }(I_0\diamond (g_1^{-1}\pi ))$$

$$={\mathrm{ad}}_{dg{g}^{-1}}^{\ast }\frac{\delta l}{\delta u},$$

$$dg{g}^{-1}=udt+\sum _{l=1}^N{\sigma }_l\circ d{W}_l^t.$$

$$\mathbf{m}(x,t)=\sum _{i=0}^N{\mathbf{p}}_i(t)\delta (\mathbf{x}-{\mathbf{q}}_i(t)),$$

$$S(\mathbf{u},\mathbf{q},\mathbf{p})=\iint l(\mathbf{u})d\mathbf{x}dt+\sum _i\int {\mathbf{p}}_i\cdot \left(\circ d{\mathbf{q}}_i-\mathbf{u}({\mathbf{q}}_i)dt+\sum _l{\sigma }_l({\mathbf{q}}_i)\circ d{W}_t^l\right).$$

$$2l(\mathbf{u})=\int \mathbf{m}(\mathbf{q},\mathbf{p})(\mathbf{x})\cdot \mathbf{u}(\mathbf{x})d\mathbf{x}=\sum _i{\mathbf{p}}_i\cdot \mathbf{u}({\mathbf{q}}_i).$$

$$S(\mathbf{q},\mathbf{p})$$

$$h(\mathbf{q},\mathbf{p})=\frac{1}{2}\sum _{i,j=1}^N({\mathbf{p}}_i\cdot {\mathbf{p}}_j)K({\mathbf{q}}_i-{\mathbf{q}}_j).$$

$${\varphi }_l(\mathbf{q},\mathbf{p}):=\sum _i{\mathbf{p}}_i\cdot {\sigma }_l({\mathbf{q}}_i).$$

$$\begin{array}{rl}d{\mathbf{q}}_i& =\frac{\partial h}{\partial {\mathbf{p}}_i}dt+\sum _l\frac{\partial {\varphi }_l}{\partial {\mathbf{p}}_i}\circ d{W}_t^l,\\ d{\mathbf{p}}_i& =-\frac{\partial h}{\partial {\mathbf{q}}_i}dt-\sum _l\frac{\partial {\varphi }_l}{\partial {\mathbf{q}}_i}\circ d{W}_t^l.\end{array}$$

$$\begin{array}{rl}d{\mathbf{q}}_i& =\sum _j{\mathbf{p}}_{j}K({\mathbf{q}}_i-{\mathbf{q}}_j)dt+\sum _l{\sigma }_l({\mathbf{q}}_i)\circ d{W}_t^l,\\ d{\mathbf{p}}_i& =-\sum _j{\mathbf{p}}_i\cdot {\mathbf{p}}_j\frac{\partial }{\partial {\mathbf{q}}_i}K({\mathbf{q}}_i-{\mathbf{q}}_j)dt-\sum _l\frac{\partial }{\partial {\mathbf{q}}_i}\left({\mathbf{p}}_i\cdot {\sigma }_l({\mathbf{q}}_i)\right)\circ d{W}_t^l.\end{array}$$

$$\begin{array}{rl}d{q}_i^{\alpha }& =\frac{\partial h}{\partial p_i^{\alpha }}dt+\sum _l{\sigma }_l^{\alpha }({\mathbf{q}}_i)\circ d{W}_t^l,\\ d{p}_i^{\alpha }& =-\frac{\partial h}{\partial q_i^{\alpha }}dt-\sum _{l,\beta }\frac{\partial {\sigma }_l^{\beta }({\mathbf{q}}_i)}{\partial q_i^{\alpha }}{p}_i^{\beta }\circ d{W}_t^l,\end{array}$$

$$\frac{d}{dt}\mathbb{P}(\mathbf{q},\mathbf{p},t)=\{h,\mathbb{P}{\}}_{\mathrm{can}}+\frac{1}{2}\sum _l\{{\varphi }_l,\{{\varphi }_l,\mathbb{P}{\}}_{\mathrm{can}}{\}}_{\mathrm{can}}:={\mathcal{L}}^{\ast }\mathbb{P},$$

$$df(\mathbf{p},\mathbf{q})=\{f,h{\}}_{\mathrm{can}}dt+\sum _l\{f,{\varphi }_l{\}}_{\mathrm{can}}\circ d{W}_t^l,$$

$$\int \mathbb{P}(\mathbf{q},\mathbf{p},t)\mathcal{L}f(\mathbf{q},\mathbf{p})d\mathbf{q}d\mathbf{p}=\int {\mathcal{L}}^{\ast }\mathbb{P}(\mathbf{q},\mathbf{p},t)f(\mathbf{q},\mathbf{p})d\mathbf{q}d\mathbf{p},$$

$$d\mathbf{x}=b(\mathbf{x},t)dt+\Sigma (\mathbf{x},t)dW,{\mathbf{x}}_0=\mathbf{u}$$

Full text

1. Introduction

In this work, we aim at modelling variability of shapes using a theory of stochastic perturbations consistent with the action of the diffeomorphism group underlying the Large Deformation Diffeomorphic Metric Mapping framework (LDDMM, see [You10]). In applications, such variability arise and can be observed, for example, when human organs are influenced by disease processes, as analysed in computational anatomy [YAM09]. Spatially independent white noise contains insufficient information to describe these large-scale variabilities of shapes. In addition, the coupling of the spatial correlations of the noise must be adapted to a variety of transformation properties of the shape spaces. The theory developed here addresses this problem by introducing spatially correlated transport noise which respects the geometric structure of the data. This method provides a new way of characterising stochastic variability of shapes using spatially correlated noise in the context of the standard LDDMM framework.

We will show that this specific type of noise can be used for all the data structures to which the LDDMM framework applies. The LDDMM theory was initiated by [Tro95, CRM96, DGM98, MTY02, BMTY05] based on the pattern theory of [Gre94]. LDDMM models the dynamics of shapes by the action of diffeomorphisms (smooth invertible transformations) on shape spaces. It gives a unified approach to shape modelling and shape analysis that is valid for a range of structures such as landmarks, curves, surfaces, images, densities or even tensor-valued images. For any such data structure, the optimal shape deformations are described via the Euler-Poincaré equation of the diffeomorphism group, usually referred to as the EPDiff equation [HMR98, HM05, YAM09]. In this work, we will show how to obtain a stochastic EPDiff equation valid for any data structure, and in particular for the finite dimensional spaces of landmarks. For this, we will follow the LDDMM derivation in [BGBHR11] based on geometric mechanics [MR99, Hol11]. This view is based on the existence of momentum maps, which are characterized by the transformation properties of the data structures for images and shapes. These momentum maps persist in the process of introducing noise into the EPDiff equation, and they thereby preserve most of the technology developed for shape analysis in the deterministic context and in computational anatomy.

This work is not the first to consider stochastic evolutions in LDDMM. Indeed, [TV12, Via13] and more recently [MS17] have already investigated the possibility of stochastic perturbations of landmark dynamics. In these works, the noise is introduced into the momentum equation, as though it was an external random force acting on each landmark independently. In [MS17], an extra dissipative force was added to balance the energy input from the noise and to make the dynamics correspond to a certain type of heat bath used in statistical physics. [SY17, SAKJ17] considered evolutions on the landmark manifold with stochastic parts being Brownian motion with respect to a Riemannian metric and estimated parameters of the models from observed data. Here, we will introduce Eulerian noise directly into the reconstruction relation used to find the deformation flows from the velocity fields, which are solutions of the EPDiff equation [HM05, You10]. As we will see, this derivation of stochastic models is compatible with variational principles, preserves the momentum map structure and yields a stochastic EPDiff equation with a novel type of multiplicative noise, depending on the gradient of the solution, as well as its magnitude. This model is based on the previous works [Hol15, ACH18], where, respectively, stochastic perturbations of infinite and finite dimensional mechanical systems were considered. The Eulerian nature of the noise discussed here implies that the noise correlation depends on the image position and not, as for example in [TV12, MS17], on the landmarks themselves. Consequently, the present method for the introduction of noise is compatible with any data structure, for any choice of its spatial correlation. We also mention the conference paper [AHPS17] in which the basic theory underlying the present work was applied to shape transformations of the corpus callosum. We discuss possibilities for including Lagrangian noise advected with the flow in contrast to the present Eulerian case, and possibilities for including non-stationary correlation statistics that responds to the evolution of advected quantities, in the conclusion of the paper.

To illustrate this framework and give an immediate demonstration of stochastic landmark dynamics, we display in Figure 1 three experiments which compare the proposed model with a stochastic forcing model, of the type studied in [TV12]. The proposed model introduces the following stochastic Hamiltonian system for the positions of the landmarks, $\mathbf{q}_{i}$, and their canonically conjugate momenta, $\mathbf{p}_{i}$,

[TABLE]

In (1.1), the $\sigma_{i}$ are prescribed functions of space which represent the spatial correlations of the noise. In Figure 1, the $\sigma_{i}$ fields are Gaussians whose variance is equal to twice their separation distance and locations are indicated by black dots. We compare this model with the system,

[TABLE]

where $\sigma$ is a constant. In this case, the noise corresponds to a stochastic force acting on the landmarks, whose corresponding Brownian motion is different for each landmark. We show on the first panel of Figure 1 that for a small number of landmarks and a large range of spatial correlations of the noise, both types of stochastic deformations in (1.1) and (1.2) visually coincide. This is shown for a simple experiment in translating a circle (from the black circle to the black dashed circle). By doubling the number of landmarks (middle panel), the dynamics of (1.2) results in small-scale noise correlation (magenta), whereas the proposed model (blue) remains equivalent to the first experiment. This figure illustrates shape evolution when the noise is Eulerian and independent of the data structure. Indeed, the limit of a large number of landmarks corresponds to a certain continuum limit, in this case corresponding to curve dynamics. Finally, in the right-most panel, we reduce the range of the spatial correlation of the noise by adding more noise fields. This arrangement allows us to qualitatively reproduce the dynamics of the equation (1.2) with the same number of landmarks as the amount of noise and its spatial correlation is similar in both cases. Indeed, the spatial correlations are dictated by the Eulerian functions $\sigma_{l}$ defined in fixed space for our model, and by the density of landmarks in the stochastically forced landmark model.

Modelling large-scale shape variability with noise is of interest for applications in computational anatomy, in which sources of variability include natural ageing, the influence of diseases such as Alzheimer’s disease, and intra-subject population scale variations. In the LDDMM context, these effects are sometimes modelled using the random orbit model [MBC*+*97]. The random orbit approach models variability in the observed data by using an ensemble of initial velocities in matching a template to a set of observations via geodesic flows, see [VMYT04]. The randomness is confined to the initial velocity as opposed to the evolving stochastic processes used in the present work. A prior can be defined by assuming a distribution of the initial velocities, and Bayesian approaches can then be used for inference of the template shape as well as additional unknown parameters [AAT07, MMTY08, ZSF13]. The stochastic model developed here can also be applied to model random warps and to generate distributions used in Bayesian shape modelling, and for coupling warps and functional variations such as those in [RSM14, KSPR17]. Indeed, because the proposed probabilistic approach assigns a likelihood to random deformations, the model can be used for general likelihood-based inference tasks.

In the present model, the observed shape variability indicates the required spatial correlation of the noise, which must be specified or inferred for each application. As this correlation is generally unknown, estimating the parameters of the correlation structure becomes an important part of the framework. We will address the problem of inferring the noise parameters by considering two different methods in the context of the representation of shapes by landmarks: The first method is based on estimating the time evolution of the probability distribution of each landmark. We will derive a set of differential equations approximating the time evolution of the complete distribution via its first moments. We can then solve the inverse problem of estimating the noise correlation from known initial and final distribution of landmarks by minimization of a certain cost function, solved using a genetic algorithm. The second method is based on an Expectation-Maximisation (EM) algorithm which can infer unknown parameters for a parametric statistical model from observed data. In this context, since only initial and final landmarks positions are observed, the full stochastic trajectories are regarded as missing information. For this algorithm, we need to estimate the likelihood of stochastic paths connecting sets of observed landmarks. We achieve this by adapting the theory of diffusion bridges to the stochastic landmark equation. As discussed in the concluding remarks, inference methods for other data structures, in particular for infinite dimensional shape representations, are not treated in this paper and left as outstanding problems for future work.

Finally, we wish to mention that multiple additional approaches for shapes analysis exist outside the LDDMM context, particularly exemplified by the Kendall shape spaces [Ken84], see also [DM16]. We focus this paper on the LDDMM framework leaving possibilities for extending the presented methods to include stochastic dynamics and noise inference in other shape analysis approaches to future work.

Plan of this work

We begin by developing a general theory of stochastic perturbations for inexact matching in section 2. We then focus on exact landmark matching in section 3, which is the simplest example of this theory. In particular, we derive the Fokker-Planck equation in section 3.2 and diffusion bridge simulation in section 3.3. In section 4, we describe the two methods we use for estimating parameters of the noise from observations. The Fokker-Planck based method is discussed in section 4.2 and the Expectation-Maximisation algorithm is treated in section 4.3. We end the paper with numerical examples in section 5, in which we investigate the effect of the noise on landmark dynamics and compare the two methods for estimating the noise amplitude.

2. Stochastic Large Deformation Matching

In this section, we will first review the geometrical framework of LDDMM, following [BGBHR11], and then introduce noise following [Hol15] to preserve the geometrical structure of LDDMM. The key ingredient for both topics is the momentum map, which we will use as the main tool for reducing the infinite dimensional equation on the diffeomorphism group to equations on shape spaces.

2.1. The Deterministic LDDMM Model

Here, we will briefly review the theory of reduction by symmetry, as applied to the LDDMM context, following the presentation of [BGBHR11]. We detail the proof of the formulas below in the next section when we include noise. Define an energy functional $E$ by

[TABLE]

where $I_{0},I_{1}\in V$ are shapes represented in a vector space $V$ on which the diffeomorphism group $\mathrm{Diff}(\mathbb{R}^{d})$ acts, $u_{t}$ is a time-dependent vector field, and $\lambda$ is a weight, or tolerance, which allows the matching to be inexact. The flow $g_{t}\in\mathrm{Diff}(\mathbb{R}^{d})$ corresponding to $u_{t}$ is found by solving the reconstruction relation

[TABLE]

and $I_{0}$ is matched against $I_{1}$ through the action $g_{1}.I_{0}$ of $g_{1}$ on $s_{0}$. The vector field $u_{t}$ can be considered an element of the Lie algebra $\mathfrak{X}(\mathbb{R}^{d})$. In the case of $I_{0},I_{1}$ being images $I:\mathbb{R}^{d}\to\mathbb{R}$, the action is by push-forward, $g.I=I\circ g^{-1}$, and when $I$ represents $N$ landmarks with positions $\mathbf{q}_{i}\in\mathbb{R}^{d}$, the action is by evaluation $g.\mathbf{q}=(g(\mathbf{q}_{1}),\ldots,g(\mathbf{q}_{N}))$ (see [BGBHR11] for more details). The group elements can act on various additional shape structures such as tensor fields.

Remark 2.1** (Nonlinear shape structures).**

This framework can be extended to structures that are not represented by a vector space $V$, such as curves or surfaces. We leave this extension for future work.

Using the calculus of variations for the functional (2.1) results in the equation of motion for $u_{t}$ of the form

[TABLE]

which is called the Euler-Poincaré equation. The operation $\mathrm{ad}^{*}$ is the coadjoint action of the Lie algebra of vector fields associated to the diffeomorphism group. The operation $\mathrm{ad}^{*}$ acts on the variations ${\delta l}/{\delta u}$, which are 1-form densities, in the dual of the Lie algebra of vector fields, under the $L^{2}$ pairing. When $l(u)$ is a norm, this equation is the geodesic equation for that norm, in the case that $\lambda=\infty$; that is, with exact matching. We will focus on this case later in section 3 when discussing landmark dynamics. Here, the inexact matching term constrains the form of the momentum $m=\frac{\partial l}{\partial u}$ to depend on the geodesic path. Following the notation of [BGBHR11], the momentum map is defined as

[TABLE]

where $g_{t,s}$ is the solution of (2.2) at time $t$ with initial conditions at time $s$, while $J_{t}^{0}=g_{t,0}I_{0}$ and $J_{t}^{1}=g_{t,1}I_{1}$. The value $J_{1}^{0}$ corresponds to the initial shape, pushed forward to time $t=1$, and $J_{1}^{1}=I_{1}$ is the target shape.

The operations $\diamond$ and $\flat$ in the momentum map formula (2.4) are defined, as follows. The Lagrangian $l$ in (2.1) may be taken as kinetic energy, which defines a scalar product and norm $l(u)=\langle u,Lu\rangle_{L^{2}}=\|u\|^{2}_{L^{2}}$ on the space of vector fields $\mathfrak{X}(\mathbb{R}^{d})$. The quantity $Lu={\delta l}/{\delta u}$ may then be regarded as the momentum conjugate to the velocity $u$. Similarly, for the image data space $V$, we define the dual space $V^{*}$ with the $L^{2}$ pairing $\langle f,I\rangle=\int_{\Omega}f(x)I(x)dx$, where $f\in V^{*}$ and $\Omega$ is the image domain $\Omega\in\mathbb{R}^{d}$. This identification defines the $\flat$ operator as $\flat:V\to V^{*}$. When an element $g_{t}$ of the diffeomorphism group acts on $V$ by push-forward, $I_{t}=g_{t}.I_{0}=(g_{t})_{*}I_{0}$, the corresponding infinitesimal action of the velocity $u$ in the Lie algebra of vector fields $u\in\mathfrak{X}(\mathbb{R}^{d})$ is given by $u.I:=[g_{t}^{*}\frac{d}{dt}(g_{t})_{*}I_{0}]_{t=0}$. In terms of this infinitesimal action, we can then define the operation $\diamond:V\times V^{*}\to\mathfrak{g}^{*}$ as

[TABLE]

A detailed derivation of this formula for the momentum map can be found in [BGBHR11].

Remark 2.2** (Solving this equation).**

We will just add here the important remark that the relation (2.4) introduces nonlocality into the problem, as the momentum implicitly depends on the value of the group at later times. This is exactly what is needed in order to solve the boundary value problem coming from the matching of images $I_{1}$ and $I_{0}$. The optimal vector field can be found with a shooting method or a gradient descent algorithm on the energy functional (2.1), see [BMTY05]. For more information about the relation of the momentum map approach of [BGBHR11] to the LDDMM approach of [BMTY05], see [BH15].

2.2. Stochastic Reduction Theory

The aim here is to introduce noise in the Euler-Poincaré equation (2.3) while preserving the momentum map (2.4); so that the noise descends to the shape spaces. Following [Hol15], we introduce noise in the reconstruction relation (2.2) and proceed with the theory of reduction by symmetry. We will focus on a noise described by a set of $J$ real-valued independent Wiener processes $W^{i}_{t}$ together with $J$ associated vector fields $\sigma_{i}\in\mathfrak{X}(\mathbb{R}^{d})$ on the data domain. We will later discuss particular forms of these fields and methods for estimating unknown parameters of the fields in the context of landmark matching.

Remark 2.3** (Dimension of the noise).**

We proceed here with a finite number of $J$ associated vector fields and finite dimensional noise while leaving possible extension to infinite dimensional noise such as done by [Via13] for later works.

We replace the reconstruction relation (2.2) by the following stochastic process

[TABLE]

where $\circ$ denotes Stratonovich integration. That is, the Lie group trajectory $g_{t}$ is now a stochastic process. With this noise construction, the previous derivations of (2.3) and (2.4) in [BGBHR11] still apply and we obtain the following result for the stochastic vector field, $u_{t}$.

Proposition 2.4**.**

Under stochastic perturbations of the form (2.6), the momentum map (2.4) persists, and the Euler-Poincaré equation takes the form

[TABLE]

Proof.

We first show that the momentum map formula (2.4) persists in the presence of noise. The key step in its computation is to prove the formula in the lemma $2.5$ of [BGBHR11] which is given by $\partial_{t}(g^{-1}\delta g)=\mathrm{Ad}_{g}\delta u$, where $\mathrm{Ad}$ is the adjoint action on the diffeomorphism group on its Lie algebra. We first compute the variations of (2.6)

[TABLE]

and then prove this formula by a direct computation

[TABLE]

This key formula is the same as in [BMTY05] and [BGBHR11] for the deterministic case. In particular, it does not explicitly depend on the Wiener processes $W_{t}^{l}$. This ensures that the momentum map formula (2.4) remains the same as in the deterministic case. The last step of the proof is to derive the stochastic Euler-Poincaré equation (2.7). This is done by computing the stochastic evolution of the momentum, given by

[TABLE]

The only time dependence is in the coadjoint action, and, by the standard formula

[TABLE]

we obtain the result

[TABLE]

where we have used the stochastic reconstruction relation (2.6) in the form

[TABLE]

∎

In summary, this stochastic perturbation of the LDDMM framework preserves the form of momentum map (2.4), although it does affect the reconstruction relation (2.6) and the Euler-Poincaré equation (2.7). As shown in [BGBHR11], various data structures fit into this framework including landmarks, images, shapes, and tensor fields. In practice, for inexact matching, a gradient descent algorithm can be used to minimise the energy functional (2.1). The noise will only appear in the evaluation of the matching cost via the reconstruction relation. The algorithm of [BMTY05] then directly applies, provided the stochastic reconstruction relation can be integrated with enough accuracy. We will not treat the full inexact matching problem here. Instead, we will study the simpler case of exact matching, where the energy functional consists only of the kinetic term.

The exact matching problem in computational anatomy possesses many parallels with the geometric approach to classical mechanics and ideal fluid dynamics. In particular, Poincaré’s fundamental paper in 1901, which started the field of geometric mechanics in finite dimensions, has recently been generalised to the stochastic case [CHR17]. In addition, the fundamental analytical properties of Euler’s fluid equations have been shown to extend to the stochastic case in [CFH17].

We expect that these advances in the analysis of SPDEs occurring in fluid dynamics and other parallel fields will inform computational anatomy, and eventually will apply to infinite-dimensional representations of shape. One reason for our optimism is that the fundamental analytical properties of incompressible Euler fluid dynamics in three dimensions have already been found in [CFH17] to persist under the introduction of the present type of stochasticity. Namely, the properties of local-in-time existence and uniqueness, as well as the Beal-Kato-Majda criterion for blow-up for the deterministic 3D Euler fluid motion equations, all persist in detail for stochastic Euler fluid motion, under the introduction of the type of stochastic Lie transport by cylindrical Stratonovich noise that we have proposed here for stochastic shape analysis.

The persistence of deterministic analytical properties in passing to the SPDEs governing stochastic 3D incompressible continuum fluid dynamics is a type of infinite-dimensional result that has not yet been proven for the evolution of shapes. The corresponding results in the analysis of SPDEs for embeddings, immersions and curves representing data structures for shape evolution, for example, have not yet been discovered, and they remain now as outstanding open problems. However, we believe that the prospects for successfully performing the necessary analysis are hopeful because the type of noise we propose here preserves the fundamental properties of diffeomorphic flow for both continuum fluids and shapes. For example, the momentum maps for the deterministic and stochastic evolution of shapes of any data structure are identical. Thus, the only difference in the present approach from the deterministic case is that the diffeomorphic time evolution of the various shape momentum maps proceeds by the action of Lie derivative by a stochastic vector field, instead of a deterministic one. Since the stochastic part of the vector field is as smooth as we wish, we are hopeful that the analytical properties for the deterministic evolution of a large class of infinite-dimensional representations of shape (such as smooth embeddings) will also persist under the introduction of the type of stochastic transport proposed here. For the remainder of the paper, we restrict ourselves to the treatment of stochastic landmark dynamics.

3. Exact stochastic Landmark Matching

In this section, we apply the previous ideas of stochastic deformation of LDDMM to exact matching with landmark dynamics. This is the simplest data structure in the LDDMM framework, and it will serve to give interesting insights into the effect of the noise in this context. Since exact matching means that the energy functional contains only a kinetic energy, the optimal vector field is found from a boundary value problem with the Euler-Poincaré equation (2.3). For exact matching, the singular momentum map for landmarks takes the simple familiar form for the reduction of the EPDiff equation (see [CH93, HM05])

[TABLE]

for $N$ landmarks with momenta $\mathbf{p}_{i}$ and positions $\mathbf{q}_{i}$, with $i=1,2,\dots,N$. A direct substitution of $u=K*m$ into the stochastic Euler-Poincaré equation (2.7) gives the stochastic landmark equations in (3.6). Here, $K$ is a given kernel corresponding to the Green’s function of the differential operator $L$ used to construct the Lagrangian. Below, we take a different approach and proceed from a variational principle in which the stochastic landmark dynamics is constrained. We refer the interested reader to for example [JS14] for a detailed exposition of this derivation in the deterministic context.

3.1. Stochastic Landmarks dynamics

Recall that for $N$ landmarks in $\mathbb{R}^{d}$, the diffeomorphism group elements $g$ act on the landmarks by evaluation of their position $g.\mathbf{q}=(g(q_{1}),\ldots,g(q_{N}))$, and the associated momentum map is (3.1). The original action functional (2.1) can be equivalently written as a constrained variational principle where the $\mathbf{p}_{i}$ play the role of Lagrange multipliers enforcing the stochastic reconstruction relation (2.6). This procedure is based on the Clebsch action principle, which for landmark dynamics has been studied for one dimensional motion of landmarks on the real line in [HT16b]

[TABLE]

Notice that only the Lagrangian depends on the spatial (Eulerian) variable $\mathbf{x}$ on the image domain. We now use the singular momentum map (3.1) which provides us with the relation

[TABLE]

This relation reduces the action functional (3.2) to the finite dimensional space of landmarks. We arrive at the action integral

[TABLE]

where the Hamiltonian only depends on the landmark variables, as

[TABLE]

The action integral in (3.3) involves the phase space Lagrangian (3.4) and the stochastic potential, given by

[TABLE]

Taking free variations of (3.3) gives the stochastic Hamilton equations in the form

[TABLE]

Explicitly, we have

[TABLE]

In coordinates, the stochastic equations (3.6) become

[TABLE]

where $\alpha,\beta$ run through the domain directions, $\alpha,\beta=1,\ldots,d$.

In order to have a unique strong solution of this stochastic differential equation, we need the drift and volatility to be Lipschitz functions with a linear growth condition after converting to Itô form, and for the volatility to be uniformly bounded, see [KS91]. This requirement is achieved when the functions $\sigma_{l}$ are twice continuously differentiable and uniformly bounded in the position variable. The latter property will hold with these functions being $C^{2}$ kernel functions. The particular form of the stochastic potential in (3.5) arises from the Legendre transformation of (3.2). The solutions of (3.8) represent the singular solutions of the stochastic EPDiff equation, corresponding to a stochastic path in the diffeomorphism group. In previous works such as [TV12, Via13, MS17], noise has been introduced additively and only in the momentum equation, corresponding to a stochastic force. Also, the noise has typically been taken to be different for each landmark, and one can interpret it having been attached to each landmark. In the present case, the noise is not additive and the Wiener processes are not related to the landmarks, but to the domain of the image. Nearby landmarks will thus be affected by a similar noise, controlled by the spatial correlations of the noise. We refer to Figure 1 in the Introduction for a numerical experiment demonstrating this effect.

Remark 3.1** (Geometric noise).**

The geometric origin of the Hamiltonian stochastic equations in (3.6) deserves a bit more explanation. In the position equation (3.6), the noise arises as the infinitesimal transformation by the action of the stochastic vector field in (2.6), namely $dgg^{-1}=udt+\sum_{l}\sigma_{l}\circ dW_{l}^{t}$, on the manifold of positions of the landmarks, which is generated by the $J$ stochastic potentials, $\Phi_{l}(\mathbf{q}_{i},\mathbf{p}_{i}):=\mathbf{p}_{i}\cdot\sigma_{l}(\mathbf{q}_{i}))$. Since this stochastic Hamiltonian is linear in the canonical momenta, the noise perturbing the evolution of the landmark positions is independent of the landmark momenta. On the other hand, the noise in the momentum equations arises as the cotangent lift of the action of the stochastic vector field $dgg^{-1}$ on the positions of the landmarks. This cotangent lift determines the action on the momentum fibres attached to the perturbed position of each of the landmarks in phase space. The cotangent lift transformation is given explicitly by the product of the momentum and the gradient of the spatial fields $\sigma_{l}$ with respect to the position $\mathbf{q}_{i}$ of the $i$-th landmark. This difference increases the effect of the noise in regions where the $\sigma_{l}$ fields have large spatial gradients, provided the landmarks are moving rapidly enough for their momenta to be non-negligible. We will see in the example that in certain cases this balance in the product of the momentum and the spatial gradient of the noise parameters can significantly affect the dynamics of the landmarks.

3.2. The Fokker-Planck Equation

In this section, we study the evolution of the probability density function of the stochastic landmarks by using the Fokker-Planck equation. This study is possible in the case of landmarks because the associated phase space is finite dimensional.

We will denote the probability density by $\mathbb{P}(\mathbf{q},\mathbf{p},t)$, on the phase space $\mathbb{R}^{2dN}$ at time $t$. The Fokker-Planck equation can be computed using standard procedures and is given in the following proposition.

Proposition 3.2**.**

The Fokker-Planck equation associated to the stochastic process (3.6) for the probability distribution $\mathbb{P}:\mathbb{R}^{2dN}\times\mathbb{R}\to\mathbb{R}$ is given by

[TABLE]

where $\{F,G\}_{\mathrm{can}}=\nabla F^{T}\mathbb{J}\nabla G$ is the canonical Poisson bracket with $\mathbb{J}=\begin{pmatrix}0&1\\ -1&0\end{pmatrix}$ and $\phi_{l}(\mathbf{q},\mathbf{p})=\sum_{i}\mathbf{p}_{i}\cdot\sigma_{l}(\mathbf{q}_{i})$ are the stochastic potentials. This formula also defines the forward Kolmogorov operator, $\mathscr{L}^{*}$.

Proof.

The proof follows the standard derivation of the Fokker-Planck equation, by taking into account the geometrical structure of the stochastic process (3.6). The time evolution of an arbitrary function $f:\mathbb{R}^{2dN}\to\mathbb{R}$ can be written as

[TABLE]

as both drift and volatility have the same Hamiltonian form in the Stratonovich formulation. We then compute the Itô correction of this stochastic process, which is can be written as a double Poisson bracket form; namely, $\frac{1}{2}\sum_{l}\{\{f,\phi_{l}\}_{\mathrm{can}},\phi_{l}\}_{\mathrm{can}}dt$. The Itô correction is the quadratic variation of the Stratonovich term in the stochastic differential equation, which equals the non-stochastic part of one half of the time derivative of the volatility (where a square Brownian motion becomes $dt$). We refer to [ACH18, CHR17] for a more detailed derivation of this formula in a general setting. Taking the expectation of the Itô process then removes the noise term and defines the forward Kolmogorov operator such that $\dot{f}=\mathscr{L}f$. By pairing this formula with the density function $\mathbb{P}(\mathbf{q},\mathbf{p},t)$ over the phase space $(\mathbf{q},\mathbf{p})$ by using the usual $L^{2}$ pairing, as

[TABLE]

we obtain the Fokker-Planck equation $\dot{\mathbb{P}}=\mathscr{L}^{*}\mathbb{P}$, which is explicitly given by (3.9) as the double bracket term is self-adjoint and the advection term anti-self-adjoint. Notice that here we have used a special property of the Poisson bracket; namely, that the Poisson bracket is also a symplectic $2$-form, which is exact and whose integral over the whole phase space vanishes, provided we choose suitable boundary conditions. We again refer to [ACH18, CHR17] for more details about this derivation. ∎

Of course, the direct study of this equation is not possible, even numerically, because of its high dimensionality. The main use here of the Fokker-Planck equation will be to understand the time evolution of uncertainties around each landmark. Indeed, for each landmark $\mathbf{q}_{i}$, the corresponding marginal distribution (integrating $\mathbb{P}$ over all the other variables) will represent the time evolution of the error on the mean trajectory of this landmark. We will show in the next section how to approximate the Fokker-Planck equation with a finite set of ordinary differential equations which describe the dynamics of the first moments of the distribution $\mathbb{P}$. This will then be used to estimate parameters of the noise fields $\sigma_{l}$ for given sets of initial and final landmarks.

Remark 3.3** (On ergodicity).**

The question of ergodicity of the process (3.6) is not relevant here, as we will only consider this process for finite times, usually between $t=0$ and $t=1$. The existence of stationary measures of the Fokker-Planck equation via Hörmander’s theorem is thus not needed. Nevertheless, we will rely on a notion of reachability in the landmark position in the next section, where we will show how to sample diffusion bridges for landmarks with fixed initial and final positions. This ensures that there exists a noise realisation which can bring any set of landmarks to any other set of landmarks. This property is weaker than the Hörmander condition and was introduced in [Sus73].

3.3. Diffusion Bridges

The transition probability and solution to the Fokker-Planck equation $\mathbb{P}(\mathbf{q},\mathbf{p},t)$ can also be estimated by Monte Carlo sampling of diffusion bridges. This approach will, in particular, be natural for maximum likelihood estimation of parameters of landmark processes using the Expectation-Maximisation (EM) algorithm that will involve expectation over unobserved landmark trajectories, or for direct optimization of the data likelihood. The EM estimation approach will be used in section 4.3. Here, we develop a theory of conditioned bridge processes for landmark dynamics which we will employ in the estimation. The approach is based on the method of [DH06] with two main modifications. The scheme and its modifications will be detailed after a short description of the general situation. Alternative methods for simulating conditioned diffusion bridges can be found in e.g. [PR12, BFS16, SMZ17].

In [DH06], a Girsanov formula [Gir60], generalized to account for unbounded drifts, is used to show that when the diffusion field $\Sigma(\mathbf{x},t)$ of an $\mathbb{R}^{d}$-valued diffusion process

[TABLE]

is uniformly invertible, the corresponding process conditioned on hitting a point $\mathbf{v}\in\mathbb{R}^{d}$ at time $T>0$ is absolutely continuous with respect to an explicitly constructed unconditioned process $\hat{\mathbf{x}}$ that will hit $\mathbf{v}$ at time $T$ a.s.. The modified process $\hat{\mathbf{x}}$ is constructed by adding an additional drift term that forces the process towards the target $\mathbf{v}$. In [DH06], this process is constructed as a modification of (3.10)

[TABLE]

Letting $P_{\mathbf{x}|\mathbf{v}}$ denote the law of $\mathbf{x}$ conditioned on hitting $\mathbf{v}$ with corresponding expectation $\mathbb{E}_{\mathbf{x}|\mathbf{v}}$, the Cameron-Martin-Girsanov theorem implies that $P_{\mathbf{x}|\mathbf{v}}$ is absolutely continuous with respect to $P_{\hat{\mathbf{x}}}$, see for example [Øk03] and the discussion in [PR12]. An explicit expression for the Radon-Nikodym derivative $dP_{\mathbf{x}|\mathbf{v}}/dP_{\hat{\mathbf{x}}}$ can be computed, and this derivative is central for using simulations of the process $\hat{\mathbf{x}}$ to compute expectations over the conditioned process $\mathbf{x}|\mathbf{v}$. In particular, as shown in [DH06], the conditioned process $\mathbf{x}|\mathbf{v}$ and the modified process $\hat{\mathbf{x}}$ are related by

[TABLE]

where $\varphi(\hat{\mathbf{x}})$ is a correction factor applied to each stochastic bridge $\hat{\mathbf{x}}$. Notice here that $f$ is a real-valued function of the stochastic path from $t=0$ to $t=T$.

Returning to landmark evolutions in the phase space $\mathbb{R}^{2dN}$, the process (3.6) has two vector variables $(\mathbf{q},\mathbf{p})$ that typically will be conditioned on hitting a fixed set of landmark positions $\mathbf{v}$ at time $T$. The conditioning thus happens only in the $\mathbf{q}$ variables by requiring $\mathbf{q}_{T}=\mathbf{v}$. To construct bridges with an approach similar to the scheme of [DH06], we need to find an appropriate extra drift term and handle the fact that the diffusion field may not be invertible in general. Recall first that the landmark process (3.6) has diffusion field

[TABLE]

where $\sigma_{j}(\mathbf{q})$ denotes the vector $(\sigma_{j}(q_{1}),\ldots,\sigma_{j}(q_{N}))^{T}$. Notice that this matrix is not square and has dimension $2dN\times J$ so that $\Sigma(\mathbf{q},\mathbf{p})\circ dW_{t}$ with $dW_{t}$ a $J$-vector corresponds to the stochastic term of (3.6). Though $\Sigma(\mathbf{q},\mathbf{p})$ couples the $\mathbf{q}$ and $\mathbf{p}$ equation, when the number of noise fields $J$ is sufficiently large, the $\mathbf{q}$ part $\Sigma_{\mathbf{q}}(\mathbf{q})$ will often be surjective as a linear map $\mathbb{R}^{J}\rightarrow\mathbb{R}^{dN}$. In this situation, by letting $\Sigma_{\mathbf{q}}(\mathbf{q})^{\dagger}$ denote the Moore-Penrose pseudo-inverse of $\Sigma_{\mathbf{q}}(\mathbf{q})$, we can construct a guiding drift term as

[TABLE]

This term, when added to the process (3.6) and when measures are taken to control the unbounded drift of (3.6), will ensure that the modified process hits $\mathbf{q}_{T}$ a.s. at time $T$. The drift term (3.14) is a direct generalization of the term added in (3.11). If $\Sigma$ had been invertible then $\Sigma\Sigma^{\dagger}={\rm Id}$ resulting in the guiding term of [DH06] used in equation (3.11). In the current non-invertible case, $\Sigma\Sigma_{\mathbf{q}}^{\dagger}(\mathbf{q}-\mathbf{v})$ uses the difference $\mathbf{q}-\mathbf{v}$ which only involves the landmark position but affects both the position and the momentum equations. We stress here the fact that the introduction of noise in the $\mathbf{q}$ equation in (3.6) is essential in our present approach. When conditioning on the $\mathbf{q}$ variable, a guided process could not directly be constructed in this way, if the noise had been introduced only in the $\mathbf{p}$ equation, as in [TV12, Via13, MS17]. The fact that this term is weighted by $\Sigma\Sigma^{\dagger}$ is also important as it allows the guiding term to be more efficient in the noisy regions of the image, where there is more freedom to deviate from the deterministic path. The guiding term can be interpreted as originating from a time-rescaled gradient flow, and with the guiding term added, the diffusion process can be seen as a stochastically perturbed gradient flow, see [AHPS17].

The guiding term (3.14) is, in practice, not always appropriate for landmarks. Because the correction is dependent only on the difference to the target in the position equation, a phenomenon of over-shooting is often observed. In such cases, the landmarks travel too fast initially due to a large momentum, strengthened by the guiding term that forces the landmarks towards $\mathbf{v}$. The high initial speed is only corrected when the time approaches $T$ and the guiding term brings the landmark back to their final position. This effect is illustrated in Figure 4 in section 5.2 and results in low values of the correction factor $\varphi(\mathbf{q},\mathbf{p})$ used to compute the expectation in (3.12). This effect tends to produce inefficient samples when approximating (3.12) by Monte Carlo sampling. As an alternative, upon letting $b(\mathbf{q},\mathbf{p})$ denote the drift term of (3.6), we employ a guided diffusion process of the form

[TABLE]

for some appropriately chosen function $\phi_{t,T}:\mathbb{R}^{2dN}\rightarrow\mathbb{R}^{dN}$ that gives an estimate of the value of $\hat{\mathbf{q}}_{T}$ using the value of the modified stochastic process $(\hat{\mathbf{q}}_{t},\hat{\mathbf{p}}_{t})$ at time $t$. The hat denotes the solution of the process (3.15), which is different from the original dynamics of the process (3.6) written without the hats. The choice $\phi_{t,T}(\hat{\mathbf{q}},\hat{\mathbf{p}}):=\hat{\mathbf{q}}$ recovers the guiding term (3.14). It would be natural to define $\phi_{t,T}(\hat{\mathbf{q}},\hat{\mathbf{p}}):=\mathbb{E}_{(\mathbf{q},\mathbf{p})}(\mathbf{q}_{T}|(\mathbf{q}_{t},\mathbf{p}_{t})=(\hat{\mathbf{q}},\hat{\mathbf{p}}))$. The resulting guiding term will only be driven by the expected amount needed at the endpoint, not from the value at time $t$. A similar choice but easier to handle is to let $\phi_{t,T}(\hat{\mathbf{q}},\hat{\mathbf{p}})$ be the solution at time $T$ of the original deterministic landmark dynamics (2.3), obtained from the initial conditions $(\hat{\mathbf{q}}_{t},\hat{\mathbf{p}}_{t})=(\hat{\mathbf{q}},\hat{\mathbf{p}})$. We will use this latter choice with a modification to ensure its time derivative is bounded. The approach is visualised in Figure 4. To ensure convergence of $\hat{\mathbf{q}}_{t}\to\mathbf{v}$ for $t\rightarrow T$, a bounded approximation $\tilde{b}$ will be chosen to replace the original unbounded drift $b$ in (3.15). As it turns out, this choice has little influence in practice.

The matrix $\Sigma(\hat{\mathbf{q}},\hat{\mathbf{p}})\Sigma_{\mathbf{q}}(\hat{\mathbf{q}})^{\dagger}$ in (3.15) only accounts for the $\mathbf{q}$ dynamics in the pseudo-inverse $\Sigma_{\mathbf{q}}(\hat{\mathbf{q}})^{\dagger}$. When the momentum is high and the noise fields $\sigma_{j}$ have high gradients, this fact can again lead to improbable sample paths. In such cases, the scheme can be further generalised by using a guiding term of the form

[TABLE]

The matrix D_{h}\big{(}\phi_{t,T}(\Sigma(\hat{\mathbf{q}},\hat{\mathbf{p}})\mathbf{h})\big{)}|_{\mathbf{h}=\mathbf{0}} is a linear approximation of the expected endpoint dynamics as a function of the noise vector $\mathbf{h}\in\mathbb{R}^{J}$. Again, with $\phi_{t,T}(\hat{\mathbf{q}},\hat{\mathbf{p}}):=\hat{\mathbf{q}}$, the original guiding term (3.14) is recovered, and the term is close to the guiding term of (3.15) when the momentum or gradients of $\sigma_{j}$ are low. We use this term for the experiments in section 5.2 involving high momentum dynamics, e.g. Figure 6.

The following result is an extension of [DH06, Theorem 5] and [Mar11, Theorem 3] to the modified guided SDE (3.15). It is the basis for the EM approach for estimating the parameters of the landmark processes developed in section 4.3. Please note that the Girsanov theorem [DH06, Thm. 1] which relates the modified and original process, does not assume that $\Sigma$ is invertible. The main analytic consequence of the non-invertibility is that the process is semi-elliptic and the transition density, therefore, cannot be bounded by Aronson’s estimation [Aro67]. Instead, we here assume continuity and boundedness of the density of $\mathbf{q}$ in small intervals of $(0,T]$ in the sense of the assumption below. We write $\mathbb{P}(\mathbf{q}_{0},\mathbf{p}_{0};\mathbf{q},\mathbf{p},t)$ for the transition density at time $t$ of a solution $(\mathbf{q},\mathbf{p})$ to (3.6) started at $(\mathbf{q}_{0},\mathbf{p}_{0})$. Similarly, when conditioning only on $\mathbf{q}$, we write $\mathbb{P}(\mathbf{q}_{0},\mathbf{p}_{0};\mathbf{q},t)=\int_{\mathbb{R}^{dN}}\mathbb{P}(\mathbf{q}_{0},\mathbf{p}_{0};\mathbf{q},\mathbf{p},t)d\mathbf{p}$.

Assumption 1**.**

For any $(\mathbf{q}_{0},\mathbf{p}_{0})$ and $(\mathbf{q},\mathbf{p})$, the process $(\mathbf{q}_{t},\mathbf{p}_{t})$ has a density $\mathbb{P}(\mathbf{q}_{0},\mathbf{p}_{0};\mathbf{q},\mathbf{p},t)$ and the map $(\mathbf{q},t)\mapsto\int_{\mathbb{R}^{dN}}g_{0}(\mathbf{q}_{0},\mathbf{p}_{0})\mathbb{P}(\mathbf{q}_{0},\mathbf{p}_{0};\mathbf{q},t)d(\mathbf{q}_{0},\mathbf{p}_{0})$ is continuous in $t$ and $\mathbf{q}$ and bounded on sets $\{(\mathbf{q},t)|s-\epsilon\leq t\leq s\}$ for $s\in(0,T]$, sufficiently small $\epsilon>0$, and any integrable function $g_{0}$.

The interpretation of Assumption 1 is that, given any distribution of initial conditions $(\mathbf{q}_{0},\mathbf{p}_{0})$ with density $g_{0}$, the resulting $\mathbf{q}$-transition density of the process is continuous and bounded in $\mathbf{q}$ and $t$. As shown in Lemma A.11, Assumption 1 can be slightly weakened if Theorem 3.4 is only used to approximate the transition density at time $T$ as opposed to expectations $\mathbb{E}[f(\mathbf{q},\mathbf{p})|\mathbf{q}_{T}=\mathbf{v}]$ for arbitrary measurable functions $f$.

We let $W(\mathbb{R}^{2dN})$ denote the Wiener space of continuous paths $[0,T]\rightarrow\mathbb{R}^{2dN}$.

Theorem 3.4**.**

Assume $\Sigma_{\mathbf{q}}(\mathbf{q}):\mathbb{R}^{J}\rightarrow\mathbb{R}^{dN}$ is surjective for all $\mathbf{q}$ with $\Sigma_{\mathbf{q}}(\mathbf{q})^{\dagger}$ bounded, and that $\Sigma$ is $C^{1,2}$, bounded, and with bounded derivatives. Let $\tilde{b}_{\mathbf{q}}$ be a bounded approximation of the $\mathbf{q}$-part of the drift $b$, and set $\tilde{b}=b+\Sigma(\mathbf{q},\mathbf{p})\Sigma_{\mathbf{q}}(\mathbf{q})^{\dagger}(\tilde{b}_{\mathbf{q}}-b_{\mathbf{q}})$. Let $v\in\mathbb{R}^{dN}$ be a point with $\mathbb{P}(\mathbf{q}_{0},\mathbf{p}_{0};\mathbf{v},t)$ positive, and let $P_{(\mathbf{q},\mathbf{p})|\mathbf{v}}$ be the law of $(\mathbf{q},\mathbf{p})\,|\,\mathbf{q}_{T}=\mathbf{v}$. Let $(\hat{\mathbf{q}},\hat{\mathbf{p}})$ be solution to (3.15), $(\hat{\mathbf{q}}_{0},\hat{\mathbf{p}}_{0})=(\mathbf{q}_{0},\mathbf{p}_{0})$ with $\varphi_{t,T}:\mathbb{R}^{2dN}\rightarrow\mathbb{R}^{dN}$ a map with $\frac{\varphi_{t,T}(\mathbf{q},\mathbf{p})-\mathbf{q}}{T-t}$ bounded on $[0,T)$. Then, for positive measurable $f:W(\mathbb{R}^{2dN})\rightarrow\mathbb{R}$,

[TABLE]

with

[TABLE]

*where A({\mathbf{q}})=\big{(}\Sigma_{\mathbf{q}}({\mathbf{q}})\Sigma_{\mathbf{q}}({\mathbf{q}})^{T}\big{)}^{-1}.

In addition,*

[TABLE]

In the Theorem, $[\cdot,\cdot]$ denotes the quadratic variation of semimartingales. As mentioned above, a bounded approximation $\tilde{b}$ must be used to replace the original drift term $b$ in (3.15). The last integrals in the expression for $\log\varphi(\mathbf{q},\mathbf{p},t)$ are results of this approximation, and the use of the map $\varphi_{t,T}$.

The result is proved in Appendix A. If $\Sigma$ had been invertible and if the guidance scheme (3.11) was used, the result of [DH06] would imply that the right hand side limit of (3.17) would equal

[TABLE]

Extending the convergence argument to the present non-invertible case is non-trivial, and we postpone investigating this to future work. For numerical computations, $\varphi(\hat{\mathbf{q}},\hat{\mathbf{p}},t)$ can be approximated by finite differences. As described later in the paper, we do this using a framework that allows symbolic evaluation of gradients and thus subsequent optimization for parameters of the processes.

4. Estimating the Spatial Correlation of the Noise

We now assume a set of $n$ observed landmark configurations $\mathbf{q}^{1},\ldots,\mathbf{q}^{n}$ at time $T$, i.e. the observations are considered realisations of the stochastic process at some positive time $T$. From this data, we aim at inferring parameters of the model. This can be both parameters of the noise fields $\sigma_{l}$ and parameters for the initial configuration $(\mathbf{q}(0),\mathbf{p}(0))$. The initial configuration can be deterministic with fixed known or unknown parameters, or it can be randomly chosen from a distribution with known or unknown parameters. We develop two different strategies for performing the inference. The first inference method in section 4.2 is a shooting method based on solving the evolution of the first moments of the probability distribution of the landmark positions while the second method in section 4.3 is based on the Expectation-Maximisation (EM) algorithm. The discussion is here in the context of landmarks, although these ideas may also apply in the more general context of section 2.

4.1. The Noise Fields

We start by discussing the form of the unknown $J$ noise fields $\sigma_{l}$. To estimate them from a finite amount of observed data, we are forced to require the fields to be specified by a finite number of parameters. A possible choice for a family of noise fields is to select $J$ linearly independent elements $\sigma_{1},\ldots,\sigma_{J}$ from a dense subset of $C^{1}(\mathbb{R}^{d},\mathbb{R}^{d})$. We here use a kernel $k$ with length-scale $r_{l}$ and a noise amplitude $\lambda_{l}\in\mathbb{R}^{d}$, that is

[TABLE]

where $\delta_{l}$ denotes the kernel positions. Possible choices for the kernel include Gaussians $k_{r_{l}}(x)=e^{-x^{2}/(2r_{l}^{2})}$, or cubic B-splines $k_{r_{l}}(x)=S_{3}(x/r_{l})$. The Gaussian kernel has the advantage of simplifying calculations of the moment equations whereas the B-spline representation is compactly supported and gives a partition of unity when used in a regular grid. Other interesting choices may include a cosine or a polynomial basis of the image domain.

In principle, the methods below allow all parameters of the noise fields to be estimated given sufficient amount of data. However, for simplicity, we will fix the length-scale and the position of the kernels. The unknown parameters for the noise can then be specified in a single vector variable $\theta=(\lambda_{1},\ldots,\lambda_{K})$. The aim of the next sections will be to estimate this vector, possibly in addition to the initial configuration $(\mathbf{q}(0),\mathbf{p}(0))$, from data using the method of moments in 4.2 and EM in 4.3, respectively.

Remark 4.1**.**

For the bridge simulation scheme, we required $\Sigma_{\mathbf{q}}(\mathbf{q})$ to be surjective as a linear map $\mathbb{R}^{J}\rightarrow\mathbb{R}^{dN}$. This assumption can be satisfied when the number of landmarks is low relative to the number of noise fields having non-zero support in the area where the landmarks reside. On the other hand, if the number of landmarks is increased while the number of noise fields is fixed, the assumption eventually cannot be satisfied. Intuitively, in such cases, the extra drift added to the bridge SDE must guide through a nonlinear submanifold of the phase space to ensure the landmarks will hit the target point $\mathbf{v}$ exactly. This limitation can be handled in three ways: (1) The method of moments as described below avoids matching individual point configurations, and it can, therefore, be used in situations where the surjectivity condition is not satisfied. (2) As discussed in remark 2.3, the noise can be made infinite dimensional. This can be done while keeping correlation structure similarly to the case with finite $J$. See also [AHPS17] for a discussion of noise in the form of a Gaussian process. (3) The bridge matching can be made inexact mimicking the inexact matching pursued in deterministic LDDMM. This could potentially relax the requirements on the extra drift term to only ensure convergence towards a given distance of the target. Inexact observations of stochastic processes are for example treated in [vdMS17].

4.2. Method of moments

We describe here our first method for estimating the parameters $\theta$ by solving a shooting problem on the space of first and second order moments. Given an estimate of the endpoint distributions $\mathbb{P}(\mathbf{q},\mathbf{p},T)$, we will solve the inverse problem which consists in using the Fokker-Planck equation (3.9) to find the values of $\theta$ such that we can reproduce the observed final distribution. Solving the Fokker-Planck equation directly is infeasible due to its high dimensionality. Instead, we will derive a set of finite dimensional equations approximating the solution of the Fokker-Planck equation (3.9) for the probability distribution $\mathbb{P}$ in terms of its first moments. This approach has been developed in the field of plasma physics for the Liouville equation, which is similar to the Fokker-Planck equation (3.9).

Remark 4.2** (Geometric moment equation).**

As the Fokker-Planck (3.9) is written in term of the canonical bracket, we could expect to be able to apply a geometrical version of the method of moments such as the one developed by [HPT07]. Although this method seems to fit the present geometric derivation of the stochastic equations, we will not use it as it is not in our case practically useful. Indeed, it requires the expansions of the Hamiltonian functions in term of the moments, but we cannot obtain here a valid expansion with a finite number of terms. This is due to the fact that the LDDMM kernel and the noise kernels cannot generally be globally approximated by finite polynomials with bounded approximation error for large distances. This would, in turn, produce spurious strong interactions between distant landmarks.

The method for approximating the Fokker-Planck that we will use here is the following. We first define the moments

[TABLE]

where we have written only two possible moments, although any combinations of $p$ and $q$ at any order are possible. In this work we will only consider moments up to the second order, that is the moments $\braket{q_{i}^{\alpha}},\braket{p_{i}^{\alpha}},\braket{q_{i}^{\alpha}q_{j}^{\beta}},\braket{q_{i}^{\alpha}p_{j}^{\beta}}$ and $\braket{p_{i}^{\alpha}p_{j}^{\beta}}$. Notice that the first moment are $(1,1)$-tensors, and the second moments are $(2,2)$-tensors, although we will only use index notation here.

We illustrate this method with the first moment $\braket{q_{i}^{\alpha}}$, which represents the mean position of the landmarks. We compute its time derivative and use the property of the Kolmogorov operator $\mathscr{L}$ defined in (3.9) to obtain

[TABLE]

We thus first need to apply the Kolmogorov operator $\mathscr{L}$ to $q_{i}^{\alpha}$ to obtain

[TABLE]

which corresponds to the $q$ part of the drift of the stochastic process with Itô correction. Similarly, for the momentum evolution, we obtain

[TABLE]

Most of the terms on the right hand side of (4.5) and (4.6) are nonlinear; so their expected value cannot be written in terms of only the first moments. This is the usual closure problem of moment equations, such as the BBGKY problem arising in many-body problems in quantum mechanics. The solution to this problem is to truncate the hierarchy of moments for a given order and consider the system of ODEs as an approximation of the complete Fokker-Planck solution. Here we will apply the so-called cluster expansion method described in [KK11]. We refer to the Appendix B.1 for more details about this method.

Apart from the first approximation $\braket{q_{i}^{\alpha}q_{j}^{\beta}}\to\braket{q_{i}^{\alpha}}\braket{q_{j}^{\beta}}$, the next order of approximation is to keep track of the correlations

[TABLE]

This quantity is also called a centred second moment as for $i=j$ it corresponds to the covariance of the probability distribution for the landmark $i$. In general, it corresponds to the correlation between the positions of two landmarks. The dynamical equation for this correlation is found from the equation of the second moment, which gives

[TABLE]

where $(i\leftrightarrow j)$ stands for the same term but with $i$ and $j$ exchanged. This equation is interesting to study in more detail, as it already gives us information about the nature of the dynamics for the spatial covariance of landmarks. Indeed, we have three types of terms with the following effects.

1. (1)

The $\sigma_{l}$-dependent terms. This first term is quadratic in the $\sigma$’s, not proportional to any linear or quadratic polynomial in $q$ or $p$. This term is a direct contribution from the noise in the $q$ equation and will have the effect of almost linearly increasing the centred covariance, wherever a $\sigma_{l}>0$. 2. (2)

The $h$-dependent terms. From the form of this term, we expect it to be proportional to a correlation. It will thus have an exponential effect on the dynamics, triggered by the linear contribution of the first term. Notice that this term only depends on the Hamiltonian, and, thus, on the interaction between landmarks. If two landmarks interact, we expect their covariance to be averaged. This term will capture their averaged covariance. 3. (3)

The $\nabla_{q}\sigma_{l}$-dependent terms. These terms are related to the noise in the $p$ equation and will account for the effect on the landmark position of the interaction of the momentum of the landmark with the gradients of the noise.

Notice that the last two types of terms describe second order effects with respect to the spatial covariance of the landmarks, as they depend linearly on the correlations. In the expansion of these nonlinear terms, the other correlations involving $p$ will appear. This means that all of the possible second order correlations must be computed. This computation is done in Appendix B, where we also approximate the expected value of the kernels as $\braket{K(\mathbf{q})}\approx K(\braket{\mathbf{q}})$. As we will see in the numerical examples in section 5, these approximations can give a reliable estimate of the landmark covariance, but this should be rigorously justified to obtain a precise estimate of the errors. Such a study is beyond the scope of this work and is left open.

Given the equations for the moment evolution, we can estimate the parameters $\theta$ by minimising the cost function

[TABLE]

where $\gamma_{1}$ and $\gamma_{2}$ are weights. We denote by $\braket{\mathbf{q}}$ and $\Delta_{2}\braket{\mathbf{qq}}$ the target first and second moments and by $\braket{\mathbf{q}}(T)$ and $\Delta_{2}\braket{\mathbf{qq}}(T)$ the estimated moments which implicitly depend on the parameters of the noise and the initial momentum. The choice of the norm is free here, and we chose a norm which only considers $i=j$ and normalises each term to $1$ so that all the covariance of the landmarks contribute equally to the cost. Other choices could be made, depending on applications. Also, the cost function may depend on other parameters, but this would make its minimisation more difficult.

To minimise the cost (4.8), we can use gradient based methods such as the BFGS algorithm. Such methods require the evaluation of the Jacobian of $C$ with respect to all of its arguments. Usually, for the estimation of the initial momentum, a linear adjoint equation is used. However, the derivative with respect to the parameters of the noise cannot be computed in this way. We will evaluate the gradients symbolically by using the Theano library in Python [The16]. To improve the efficiency of the algorithm, we first match the mean final position, by only updating the initial momentum. Then, with this initial condition, we match for both first and second moments and update the initial momentum as well as the parameters $\lambda_{l}$. As we will see in the numerical experiments in section 5, gradient-based methods are not optimal, and genetic algorithms, such as the differential evolution algorithm of [SP97] designed for global minimizations, turn out to perform better.

4.3. Maximum Likelihood and Expectation-Maximization

We now describe how to estimate the unknown parameters collected in the variable $\theta$ by a maximum likelihood estimation based on the expectation-maximisation (EM) algorithm of [DLR77]. The likelihood of $n$ independent observations $(\mathbf{q}^{1},\ldots,\mathbf{q}^{n})$ at time $T$ given parameters $\theta$ takes the form

[TABLE]

The parameters $\theta$ can be estimated by maximizing the likelihood, that is by letting

[TABLE]

For this, the likelihood could be directly computed by numerical approximation of $\mathbb{P}_{\theta}(\mathbf{q}_{i},T)$ using an approximation of the Fokker-Planck equation (3.9). Alternatively, the fact that the stochastic process is only sampled at time $T$ suggests a missing data approach that marginalises out the unobserved trajectories up to time $T$. Let $(\mathbf{q},\mathbf{p};\theta)$ denote the stochastic landmark process with parameters $\theta$, and let $P(\mathbf{q},\mathbf{p};\theta)$ denote its law. Let $\mathcal{L}(\mathbf{q},\mathbf{p};\theta)$ denote the likelihood of the entire stochastic path for a given realisation of the noise, and computed with respect to the parameter $\theta$. Notice that this likelihood is only defined for finite time discretizations of the process and there is no notion of path density for the infinite dimensional process. We thus proceed formally, while noting that the approach can be justified rigorously, see e.g. [DS08]. An alternative approach is to optimize the likelihood (4.9) directly using (3.18). This is pursued in [SAKJ17].

The EM algorithm finds a sequence of parameter estimates $\{\theta_{k}\}$ converging to a $\hat{\theta}$ by iterating over the following two steps:

1. (1)

Expectation: Compute the expected value of the log-likelihood given the previous parameter estimate $\theta_{k-1}$:

[TABLE]

The expectation (4.10) over the process conditioned on the observations $\mathbf{q}_{i}$ integrates the likelihood over all sample paths reaching $\mathbf{q}_{i}$. For this, we employ the bridge simulation approach developed in section 3.3. For each $\mathbf{q}^{i}$, we thus exchange $(\mathbf{q}_{t},\mathbf{p}_{t};\theta)$ with a guided process $(\hat{\mathbf{q}},\hat{\mathbf{p}};\theta,\mathbf{q}^{i})$ and use the equality (3.17) from Theorem 3.4. The expectation on the right-hand side of (3.17) can be approximated by drawing samples from the guided process. Note that the correction factor $\varphi(\mathbf{q},\mathbf{p}|\theta_{k-1},\mathbf{q}_{i})$ makes the approach equal to importance sampling of the conditioned process with the guided process as proposal distribution. 2. (2)

Maximisation: Find the new parameter estimate

[TABLE]

The maximisation step can be approximated by updating $\theta_{k}$ such that it increases $Q(\theta|\theta_{k-1})$ instead of maximising it. This is the approach of the generalised EM algorithm [NH98]. The update of $\theta$ is thus computed by taking a gradient step

[TABLE]

where $\epsilon>0$. The gradient which is evaluated for each of the sampled paths can be computed symbolically using the Theano library [The16]. Theano allows the entire computational chain from the definition of the Hamiltonian and noise fields to the time-discrete stochastic integration to be specified symbolically. The framework can therefore automatically derive gradients symbolically before the expressions are compiled to efficient numerical code. See also [KAS17] for more details on the use of Theano for differential geometric and stochastic computations.

The resulting estimation algorithm is listed in Algorithm 1. For each $\mathbf{q}^{i}$, the expectation $\mathbb{E}_{\theta_{k-1}}[\log\mathbb{P}_{\theta}(\mathbf{q},\mathbf{p}|\mathbf{q}_{i})]$ is estimated by sampling $N_{\mathrm{bridges}}$ bridges. The algorithm can perform a fixed number $K$ of updates to the estimate $\theta_{k}$ or stop at convergence.

5. Numerical examples

We now present several numerical tests of the stochastic perturbation of the landmark dynamics. In particular, we want to illustrate aspects of the effect of the noise on the landmarks and test the algorithms for estimation of the spatial correlation of the noise. We will focus here on synthetic examples and refer to [AHPS17] for an application of the methods on a dataset of Corpora Callosa shapes represented by 77 landmarks. The numerical simulations of this work have been done in Python, using the symbolic computation framework Theano [The16]. The code is available from the public repository https://bitbucket.org/stefansommer/stochlandyn. See also [KAS17] for additional details.

5.1. Solution of the Fokker-Planck equation

We first consider a simple experiment with a single landmark, subjected to a square array of noise fields with Gaussian noise kernel. To a first order approximation, the mean trajectory of the landmark is a straight line with constant momenta as the Hamiltonian is a pure kinetic energy.

This experiment is displayed in Figure 2(a) where we used two arrays of four by four noise fields with either $\lambda_{l}=(0.08,0)$ or $\lambda_{l}=(0,0.08)$ and three values of the noise radius $r_{l}=0.5,0.05,0.03$. For large values of $r_{l}$, the noise is mostly uniform and the gradients of the $\sigma_{l}$ are negligible. The only term contributing to the final covariance of the landmark is therefore $\Braket{\sigma_{l}^{\alpha}(\mathbf{q}_{i})\sigma_{l}^{\beta}(\mathbf{q}_{j})}$. Notice that because there is only one landmark, thus a linear drift, the deterministic part does not affect the covariance. This term only has a linear effect on the covariance which is thus an ellipse proportional to the noise fields. Here the noise has equal strength in both the $x$ and $y$ coordinate thus we observe a circle. For smaller values of $r_{l}$, the gradient of $\sigma_{l}$ is large enough for the other term in the momentum equation which couples the momentum and the gradient of $\sigma$ to affect the moment dynamics. This effect is shown in Figure 2(a) where the covariance has a larger value in the direction of the gradient of $\sigma_{l}$ than in the other directions. This is explained by the fact that this coupling is of the form $\frac{\partial}{\partial q_{i}}(\sigma_{l}(\mathbf{q}_{i})\cdot\mathbf{p}_{i})$, thus the ellipse is in the direction of the gradient, not the momenta. Notice that there should be some noise in the direction of the momenta for this term to have an effect.

Using the same experiment, we compared the estimation of the covariance from the moment equation with a direct sampling obtained by solving the stochastic landmarks equations. We did this experiment for $r_{l}=0.5$ and $r_{l}=0.03$ in Figure 2(b),2(c). The left panel with $r_{l}=0.5$ shows an excellent agreement between the two methods but the right panel with $r_{l}=0.03$ shows differences. This type of error in the estimation of the covariance is explained by the fact that the final distribution has a large skewness. This effect is not captured by the moment equations as we neglected the effects of order higher than $2$, and the skewness is a third order effect described by terms such as $\Delta_{3}\braket{q_{i}^{\alpha}q_{j}^{\beta}q_{k}^{\gamma}}$. Nevertheless, the final covariance is close enough to the correct one to be able to use it in the estimation of the noise fields. This demonstrates that even in rather extreme cases, which are not realistic for applications, the second order approximation used to derive the moment equation still produces reliable results.

We did a similar experiment but with $5$ interacting landmarks arranged in an ellipse configuration and with initial conditions obtained from the deterministic shooting method such that the endpoint of the deterministic landmark equations match another ellipse. We display these experiments in Figure 3 with the same noise as in the previous tests and with $r_{l}=0.2$. We modified here the landmark interaction length scale $\alpha$ from $\alpha=0.02$ (no-interactions) to $\alpha=0.2$ (neighbours interactions) to see the effect of the noise with the landmark interactions. Due to the different length scales, the trajectories to the target ellipse are slightly different so the landmarks will be subject to different noise. The larger length scale has the effect of reducing the differences between the covariances of interacting landmarks.

5.2. Bridge Sampling

Here, we aim at visualising the effect of the constructed bridge sampling scheme. In Figure 4, the effect of the guiding term is visualized on a sample path. At $t=T/2$, the predicted endpoint $\phi_{t,T}(\hat{\mathbf{q}}(t),\hat{\mathbf{p}}(t))$ is calculated and the difference $\phi_{t,T}(\hat{\mathbf{q}},\hat{\mathbf{p}})-\mathbf{v}$ is used to guide the evolution of the path towards the target $\mathbf{v}$. The guiding term ensures that $\hat{\mathbf{q}}$ will hit $\mathbf{v}$ almost surely at time $T$. Notice that the difference $\phi_{t,T}(\hat{\mathbf{q}},\hat{\mathbf{p}})-\mathbf{v}$ is generally much smaller than the difference $\hat{\mathbf{q}}-\mathbf{v}$. The introduction of $\phi_{t,T}$ therefore implies that the process is modified less giving more likely bridges. Without $\phi_{t,T}$, the process is generally attracted too quickly towards the target as can be seen by the landmarks at $t=0.5$ being almost at their final positions in Figure 4 (b). The path thus overshoots the target. This effect is not present when using $\phi_{t,T}$ in Figure 4 (a).

5.3. Estimating the noise amplitudes

We here aim at estimating the noise amplitude from sampled data using both the method of moments and maximum likelihood.

We first use the genetic algorithm of [SP97] called differential evolution algorithm to minimise the cost function $C$ in (4.8). This algorithm has in experiments proven successful in avoiding local minima during the optimisation. We compared it with the standard BFGS gradient descent algorithm with a single landmark in Figure 5. This algorithm relies on the Jacobian of the cost functional computed symbolically using the Theano package of [The16]. It is able to estimate the noise amplitude along the trajectory of the landmark where the signal from the gradient of $C$ is the strongest. For the other regions of the image, the algorithm cannot detect any signal to update the noise fields. The genetic algorithm can overcome this issue as it is based on evolving a population of solutions which uniformly cover the entire parameter space. In this way, the solution obtained is a better approximation of the global minimum of $C$. It is interesting that even if the final moment of figure 5 is well matched with the genetic algorithm, the noise amplitude is not perfectly recovered. This illustrates the expected degeneracy of this model for a low number of landmarks. When more landmarks are added, the noise amplitude estimation is closer to the expected one, see figure 7 below. In these experiments, we set the initial variance of the momentum, and the position/momentum correlation to [math] for simplicity, and because we used these values to generate the final variance. In practice, one may expect to have other prior for the initial variance of the momenta or can try to find it as an unknown parameter of the problem. Having them as unknown may results in a large parameter space, thus simplifications such as all the landmarks have the same initial variance in the momentum could be used. We leave such investigation for later when applied to real data, with a possible meaningful prior.

In Figure 6, the same experiment is performed with MLE and the bridge sampling scheme. The noise kernels are in this experiment cubic B-splines placed in a grid providing a partition of unity. In the optimisation, $\lambda_{l}$ are fixed to be equal for all $l=1,\ldots,J$ implying that the total noise variance will be uniform at each point of the domain. The figure shows the experiment performed with low momentum (Figure 6 (a)) and high momentum (Figure 6 (b)). In the low momentum case, the noise parameters are estimated correctly and the sample covariance with the estimated parameters matches the covariance of the original samples. The SDE (3.15) is here used for the bridge sampling scheme. In contrast to the previous method, the algorithm is now optimising for the maximum likelihood of the samples and not directly for matching the final covariance. A higher difference in the endpoint covariance is, therefore, to be expected.

With higher initial momentum, the coupling between the guidance and noise makes the scheme (3.15) overestimate the variance. Instead, the guidance term (3.16) is used. Notice that even though the sample covariance with the estimated parameters matches the covariance of the original samples, the estimated $\lambda_{l}$ are different than the original values. This indicates that the maximum likelihood estimate of the parameters may not match the original setting in the highly nonlinear case occurring when the coupling between noise and momentum is high. Because of the nonlinearity, the noise is able to generate horizontal variation in the position of the final the landmark even though the variation with the estimated parameters are mainly vertical along the trajectory.

Figure 7 and Figure 8 show the result of noise estimation using different configurations of the ellipse and both the method of moments and MLE. The noise parameters $\lambda_{l}$ are allowed to vary with $l$ in both cases giving spatially non-uniform noise amplitude. The algorithms find the correct noise parameters in the areas covered by the landmark trajectories.

6. Discussion and Outlook

As the first topic of this work, we raised the issue of how to include stochasticity and uncertainty in the framework of large deformation matching in a systematic and geometrically consistent way. In section 2, we exposed a general theory of stochastic deformations in the LDDMM framework, based on the momentum map representation of images in [BGBHR11], by introducing spatially correlated time-dependent noise in the reconstruction relation that is used to compute the deformation map from its velocity field. By taking this approach, we have preserved most of the advantages of the theory of reduction by symmetry. In particular, we have preserved the capability of applying this stochastic model to general data structures. The dynamical equation is the stochastic EPDiff equation, in which the noise appears in a certain multiplicative form with spatial correlations encoded in a set of spatially dependent functions $\sigma_{l}$. The key feature of this noise is that the structure of the original equation provided by the theory of reduction by symmetry still remains. In particular, the persistence of the momentum map allows for both exact and inexact matching in this stochastic context.

The question of local in time existence and uniqueness of this equation is important, but it is not treated in this work. We refer to [BFM16] for such a study for the 2D Euler equation and to [CFH17] for the 3D case. Another possible extension would be to consider an infinite number of $\sigma_{l}$ fields with an infinite dimensional Wiener process for the stochastic EPDiff equation as investigated in [Via13], also in the context of stochastic shape analysis. We considered have time-independent $\sigma_{l}$ fields. However, there are several approaches for making these fields time dependent besides simply prescribing them as functions of time. Some of these other approaches were derived by [GBH18] in the context of stochastic fluid dynamics. In particular, the idea of having the noise fields being carried by the deformation could be of interest in this context as well. Yet another possibility would be to introduce two different types of noise fields, one modelling small-scale noise correlations and the other one for larger scale noise correlations. In this case, it would make sense for the small scale variability to be advected by the large-scale deformation, as in the multi-scale model of [HT12].

After defining the general model in section 2, we applied it to exact landmark matching in section 3, which is the simplest non-trivial application of the LDDMM framework. This approach allowed investigation of the effects of the noise on large deformation matching in a finite dimensional model. Introducing the noise in both the momentum and the position equations of the landmarks made the landmark trajectories rougher than they would have been, otherwise, had the noise been only in the momentum equation. The noise in the position equation also increased the flexibility for controlling the landmark trajectories. This flexibility was used to derive a scheme for simulating diffusion bridges with corresponding sampling correction factor that allowed evaluation of expectations with respect to the original conditioned landmark dynamics. In addition, we used the finite dimensionality of the system to derive the Fokker-Planck equation and apply it to the dynamics of moments of the probability distribution function.

Some modifications to the standard theory of diffusion bridges were made to accommodate the case of landmark dynamics and to improve the speed and accuracy of the estimation of expectations over conditioned landmarks trajectories. The landmarks represent the simplest cases for numerical shape analysis, especially in the context of stochastic systems. We used a simple Heun method to solve the stochastic landmark equations. Higher order integration schemes could have been used, such as the stochastic variational integrators of [HT16a]. The next step in extending the landmark example is to allow for inexact matching and to study the trade-off between the effect of noise and the tolerance of the matching.

Several issues regarding ergodicity and other properties of the Kolmogorov operator were left open in this paper, whose future treatments could add to the theoretical understanding of the model. Finally, the stochastic LDDMM framework can be applied to other types of data structures, in particular to images with inexact matching as originally done in [BMTY05]. Studying the effects of the stochastic model on other nonlinear data structures such as curves or surfaces would also be of great interest for future works.

As a second topic, we raised the issue of determining the noise correlation from data sets which would allow the theory of stochastic deformations to be used with observed data. We developed two independent methods which we implemented and applied to several test examples. First, the moment equation allows matching of the sample moments. It is deterministic, making optimisation of the noise parameters stable and efficient, and it does not require special conditions on the noise fields. Its accuracy depends on the approximation order in the moment equation. Scaling the moment equation to a large number of landmarks or continuous shapes such as curves may be challenging as well as optimising for a high number of unknown parameters. In the landmark experiments we presented above, this approach allowed us to reliably estimate the underlying noise, but an extension of this method to infinite dimensional representations of shapes is not possible unless a discretized version of the equations is used. For this method, we also made two approximations that could possibly be improved elsewhere. One of them is the truncation to retain only second order terms in the moment equation, and the other is to approximate the expectation of a kernel function as the kernel of the expected values. Both approximations were shown to work well for cases with small enough noise, which would be the case in most applications. Finally, it is also important to notice that we did not use the freedom of the initial value of the variance of the momentum and the position/momentum correlation. These parameters could either be inferred using this scheme (with a larger parameter space) or be obtained by using other information about the data.

The second method is the MLE optimisation, a Monte Carlo method which evaluates expectations over conditioned stochastic trajectories. The bridge sampling scheme we used requires the noise fields to span the entire $\mathbf{q}$-space to allow guiding the landmarks towards their target. With high nonlinearity as may happen with large initial momentum and high gradients of the noise fields, guiding the trajectories towards their target with high-probability bridges can be challenging. In general, the stochastic nature of the algorithm makes it harder to control than the matching provided by the moment equation. The bridge sampling scheme can be interpreted as a gradient flow, as discussed in [AHPS17] when applied to images. It allows the likelihood of observed images to be evaluated without a prior image registration step. The method may thus be applicable to image analysis problems, and more generally for inexact matching of shapes in which case the requirement of the noise to span the $\mathbf{q}$-space may be relaxed.

The inference of noise parameters treated here can be extended to more general statistical inference problems on shape spaces. Inferring the initial $\mathbf{q}_{0}$ positions can be regarded as estimating a most-likely mean, thereby drawing similarities to the Frechét mean [Fré48] and to means defined by the maximum likelihood of probability distributions in nonlinear spaces [Som15]. When generalised to images, the approach can be used for simultaneous estimation of template images [JDJG04], possible time-dependent transformations in the momentum as caused for example by disease processes [MF12], and population variation in the spatial noise correlation.

It is possible to generalise the stochastic equations we have introduced here to allow for time-dependent noise amplitude as done in [GBH18] for fluid dynamics. In this case, the noise fields could be advected by the diffeomorphism and only the initial condition of the noise field would have to be inferred. This requires the choice of a meaningful advection scheme. By construction of its metric LDDMM is right-invariant, and the flow energy is therefore measured in Eulerian coordinates. This leads us to define stochastic flows that are compatible with this right-invariant geometry thus giving noise in Eulerian coordinates. In the deterministic setting, left-invariant metrics [SRV15] provide a Lagrangian view of the metric that thus, in a medical context, follows the advected anatomy. We leave it as an open and very relevant problem to consider advected, or left-invariant Lagrangian noise.

Extending the inference methods presented here to other data structures, in particular to infinite dimensional shapes spaces, would again constitute an interesting future direction. As discussed in detail at the end of section 2, we believe that the methods presented here with suitable modifications can be applied also for infinite-dimensional representations of shapes, and that additional methods could be introduced, such as stochastic filtering for further data assimilation of the results in infinite-dimensional cases, see e.g. [BC].

Acknowledgements

We are grateful to M. Bruveris, M. Bauer, A Pai, N. Ganaba, C. Tronci, M.R. Schauer, T. Tyranowski and F. Van Der Meulen for helpful discussions of this material and to the anonymous referees for thoughtful comments which improved the exposition of this paper. AA and DH were partially supported by the European Research Council Advanced Grant 267382 FCCA held by DH. DH is also grateful for support from EPSRC Grant EP/N023781/1. AA acknowledges partial support from an Imperial College London Roth Award and the EPSRC through award EP/N014529/1 funding the EPSRC Centre for Mathematics of Precision Healthcare. SS is partially supported by the CSGB Centre for Stochastic Geometry and Advanced Bioimaging funded by a grant from the Villum Foundation. The research was partially completed while the authors were visiting the Institute for Mathematical Sciences, National University of Singapore in 2016.

Appendix A Bridge Sampling

We here follow [DH06] and the later paper [Mar11] to argue for the almost sure hitting of a target $\mathbf{v}$ for the guided process (3.15) and to find the correction term $\varphi(\mathbf{q},\mathbf{p},t)$. For completeness, we will explicitly derive the correction term following the program in [DH06, Theorem 5]. The guided SDE (3.15) differs from the previous schemes in using the function $\phi_{t,T}:\mathbb{R}^{2dN}\rightarrow\mathbb{R}^{dN}$ to predict the endpoint, and, importantly, in that the diffusion field $\Sigma$ is not invertible resulting in a semi-elliptic diffusion. We handle the first issue by repeated application of the Girsanov theorem. This also accounts for the unboundedness of the drift term $b(\hat{\mathbf{q}},\hat{\mathbf{p}})$ coming from the momentum of the landmarks in the same way as [Mar11]. We do not here argue for the $t\rightarrow T$ limit of the expectation of the correction term that in the elliptic case follows from [DH06, Mar11].

Let $b(\mathbf{q},\mathbf{p})$ be the drift in (3.6) in Itô form. Because this drift is unbounded, we construct $\tilde{b}(\mathbf{q},\mathbf{p})$ in Theorem 3.4 to be an approximation so that the $\mathbf{q}$-part $\tilde{b}_{\mathbf{q}}$ is bounded on $\mathbb{R}^{dN}$. To construct a map $\phi_{t,T}:\mathbb{R}^{2dN}\rightarrow\mathbb{R}^{dN}$ satisfying the conditions of the theorem, let $\phi_{t,T}$ be the $\mathbf{q}$-part of the time $T$ solution to the ODE $\partial_{t}(\mathbf{q}_{t},\mathbf{p}_{t})=\tilde{b}(\mathbf{q}_{t},\mathbf{p}_{t})$ started at time $t$ with initial conditions $(\mathbf{q},\mathbf{p})$. This ODE corresponds to the deterministic ODE (2.3), however using the drift approximation to ensure $\partial_{t}\phi_{t,T}(\mathbf{q},\mathbf{p})$ is bounded. Then the process $\frac{\tilde{\phi}_{t,T}(\hat{\mathbf{q}},\hat{\mathbf{p}})-\hat{\mathbf{q}}}{T-t}$ is defined, bounded and continuous on $[0,T]$.

The SDE

[TABLE]

differs from the Itô form of the SDE (3.15) by

[TABLE]

As argued by [Mar11], A.1 has a unique solution satisfying $\lim_{t\rightarrow T}\hat{\mathbf{q}}=v$ a.s., and the processes A.1 and (3.15) are absolutely continuous with respect to each other. The correction term $\varphi(\mathbf{q},\mathbf{p},t)$ can be derived from [Mar11, Theorem 3] and the difference A.2. For completeness, we give the derivation in the landmark case that proves Theorem 3.4 below.

Proof of Theorem 3.4.

Let $f:W(\mathbb{R}^{2dN})\rightarrow\mathbb{R}$ be a non-negative measurable function on $[0,t]$, $t<T$. Following [DH06], we define

[TABLE]

noting that in the present case, we use the pseudo-inverse $\Sigma_{\mathbf{q}}(\mathbf{q})^{\dagger}$ in $h$ since $\Sigma_{\mathbf{q}}(\mathbf{q})$ is not invertible. Let now $(\tilde{\mathbf{q}},\tilde{\mathbf{p}})$ be a solution to the SDE

[TABLE]

From the Girsanov theorem with unbounded drift [DH06, Thm. 1], we have

[TABLE]

where

[TABLE]

We now define an intermediate function

[TABLE]

and compute

[TABLE]

Applying the product rule, we obtain for the second term

[TABLE]

Writing the process $\frac{1}{2}g(\tilde{\mathbf{q}},t)$ in integral form,

[TABLE]

Note that the first term $\int_{0}^{t}\frac{(\tilde{\mathbf{q}}-\mathbf{v})^{T}A(\tilde{\mathbf{q}})(\tilde{\mathbf{q}}-\mathbf{v})}{2(T-s)^{2}}ds$ of the right hand side is the negative of the term $-\frac{1}{2}\int_{0}^{t}\|h(\tilde{\mathbf{q}},\mathbf{p},s)\|^{2}ds$ in (A.5). The second term expands to

[TABLE]

where the second term of the right hand side is the negative of $\int_{0}^{t}h^{T}(\tilde{\mathbf{q}},\mathbf{p},s)dW$. Rearranging terms and inserting in (A.5),

[TABLE]

We can now use the Girsanov theorem again to change the drift from (A.3) to (A.1). For this, we define

[TABLE]

Then let $\varphi(\mathbf{q},\mathbf{p},t)$ be the function satisfying

[TABLE]

Now (A.4) and the definition of $\log\varphi$ gives

[TABLE]

with $(\hat{\mathbf{q}},\hat{\mathbf{p}})$ a solution to (3.6). Thus

[TABLE]

where convergence of the right-hand side limit to the conditioned process follows from Lemma A.10 below. The limit expression (3.18) for the density follows from Lemma A.11 and (A.8). ∎

The lemmas A.10 and A.11 below follow [DH06] and [Mar11] with minor modifications to clarify where the assumption 1 on the density of the process $(\mathbf{q}_{t},\mathbf{p}_{t})$ is needed in the semi-elliptic case.

Lemma A.1**.**

Let $(\mathbf{q},\mathbf{p})$ be a solution to (3.6) satisfying assumption 1 and the conditions of Theorem 3.4. Let $0<t_{1}<t_{2}<\cdots<t_{n}<T$ be a finite set of time point in $[0,T]$ and let $f\in C_{b}(\mathbb{R}^{n2dN})$. Let $\psi_{t}$ be the process

[TABLE]

with $g$ as defined above. Then

[TABLE]

Proof.

Following [DH06], we write

[TABLE]

with

[TABLE]

Note that $\Phi_{f}$ is continuous by assumption. We now apply a change of variable $\mathbf{q}=\mathbf{v}+(T-t)^{\frac{1}{2}}\mathbf{q}^{\prime}$ to get

[TABLE]

From assumption 1, $\Phi_{f}(t,\mathbf{q})$ is continuous and bounded. Because $\Sigma$ is bounded, $e^{-\frac{\|\Sigma_{\mathbf{q}}(\mathbf{v}+(T-t)^{\frac{1}{2}}\mathbf{q}^{\prime})^{\dagger}(\mathbf{q}^{\prime})\|^{2}}{2}}\leq e^{-\frac{c\|\mathbf{q}^{\prime}\|^{2}}{2}}$ for some constant $c$ and the dominated convergence theorem implies the limit. We conclude that

[TABLE]

The result now follows from the definition of $\Phi_{f}$, see [DH06]. ∎

If only the density $\mathbb{P}(\mathbf{q}_{0},\mathbf{p}_{0};\mathbf{q},t)$ is of interest, the following result holds assuming only continuity and boundedness of $\mathbb{P}(\mathbf{q}_{0},\mathbf{p}_{0};\mathbf{q},t)$ for fixed initial conditions $(\mathbf{q}_{0},\mathbf{p}_{0})$.

Lemma A.2**.**

Let $(\mathbf{q},\mathbf{p})$ be a solution to (3.6) with the conditions of Theorem 3.4 and assume the process has a density $\mathbb{P}(\mathbf{q}_{0},\mathbf{p}_{0};\mathbf{q},\mathbf{p},t)$ and that $\mathbb{P}(\mathbf{q}_{0},\mathbf{p}_{0};\mathbf{q},t)$ is continuous in $t$ and $\mathbf{q}$ and bounded on $\{(\mathbf{q},t)|T-\epsilon\leq t\leq T\}$. Let $\psi_{t}$ be the process defined above. Then

[TABLE]

Proof.

The result follows from the convergence of $\Phi_{1}(t,\mathbf{q})$ in Lemma A.10. ∎

Appendix B Moment equation for stochastic landmark

B.1. Cluster expansion method

We explain the basics of this method, which can be found in more details in, for example, [KK11] with application in the context of semiconductor physics. This method is used when one seeks the dynamics of the expected value of $N$ particles that we will write here $\Braket{N}$. One cannot solve the complete system, especially if the number of particles is large, thus we want to approximate the expected value of products in term of only a few independent variables. For this, we apply the cluster expansion, which begins by writing

[TABLE]

The next decomposition is

[TABLE]

and so on and so forth. We then only compute the dynamics for the singlets $\braket{1}$ and the correlations, up to some chosen order. In the sequel, we will only consider the doublet correlations $\Delta_{2}$, and in this case, we have the general decomposition

[TABLE]

In the context of quantum mechanics, where the particle operators do not commute, extra care is needed especially for the sign of the term. Here we will consider $q_{i}^{\alpha}$ and $p_{i}^{\alpha}$ as our particles, and as they commute, the expansions are simpler than in [KK11]. We directly compute two of them for illustration, up to quadratic order,

[TABLE]

This sort of expansion can fit a more geometrical framework, where the final equations for the first moment will preserve the original structure of the equations. This was developed first in [HLS90] and later in [HPT07, HPT10]. We will not use this method here for a good reason related to the form of the equations. A key step in these papers is to expand the expected value of the Hamiltonian in terms of a finite number of moments, to enable computation of the equation of motion. In our case, the Hamiltonian has a kernel function, which generally cannot be expanded in a finite sum of polynomial terms. By doing the computations directly, we will be able to do another approximation for the kernels, that is, we will assume that they commute with the operation of expectation. A more subtle approximation can be done using the Heaviside function but will give a much larger number of terms in the expansion, see appendix B.6 for no clear improvements of the solution.

To perform this expansion on the Fokker-Planck equation associated to the landmark dynamics we will use several simplifications:

• •

Gaussian noise fields $\sigma_{l}$ in (4.1),

• •

for a kernel $K(x)$, we will assume that $\braket{K(x)}\approx K(\braket{x})$ and

• •

only the second order correlations $\Delta_{2}$ will be considered in this expansion.

These assumptions can be relaxed but the resulting equation may be difficult to compute.

B.2. First moments

Recall the backward Kolmogorov operator on $q_{i}^{\alpha}$

[TABLE]

which is used to compute the time evolution of the singlet

[TABLE]

In this case, the equation only depends on the singlet of the momentum variable. We thus compute

[TABLE]

which similarly gives the time evolution of the momentum singlet in two terms as

[TABLE]

where

[TABLE]

We then expand $A_{p}$ further using the cluster-expansion method on the triplet to get

[TABLE]

Already this term depends on the mixed correlations which we will compute shortly, but we first compute the position correlation.

B.3. $\braket{qq}$ correlation

Recall the formula of the Kolmogorov operator applied to $q_{i}^{\alpha}q_{j}^{\beta}$,

[TABLE]

which together with (4.5) gives the time evolution of the position correlation in the form

[TABLE]

where

[TABLE]

We will denote by $A$ the terms corresponding to the drift, by $B$ the terms which are not present in the first moments equation, and by $C$ the other terms which only depend on the noise and the derivative of the noise fields. We proceed by first approximating the expectation of the kernels to get

[TABLE]

where we also used the explicit form of $\sigma_{l}$ as a Gaussian and its derivative. We will now approximate the $A_{qq}$ term to get

[TABLE]

It is now clear that the $B$ term will linearly increase the position correlation, which will then exponentially increase by the $C$ term and be affected by the momentum-position correlation by the $A$ term. We now proceed by computing the momentum correlation.

B.4. $\Braket{pp}$ correlation

We compute the Kolmogorov operator on $p_{i}^{\alpha}p_{j}^{\beta}$ to get

[TABLE]

and using (B.5), we obtain the time evolution of the correlation in three terms as

[TABLE]

where

[TABLE]

We first approximate

[TABLE]

The last two terms cancel as they are symmetric under the transpose operation because of the sum on the free indices, thus the $C$ term is

[TABLE]

We proceed with the $B_{pp}$ term, which is also symmetric under the transpose operation, thus giving the approximation

[TABLE]

We treat the two Hamiltonian terms separately by first writing them explicitly as

[TABLE]

We expand the first term to arrive at

[TABLE]

and the second term

[TABLE]

This term cancels the terms of the $A_{pp}^{1}$ proportional to $\braket{p_{i}^{\alpha}}$ to give the approximation

[TABLE]

We end this computation by approximating the dynamics of the mixed correlation.

B.5. $\Braket{pq}$ correlation

We compute

[TABLE]

Then, using (4.5) and (B.5) we obtain the time evolution of $\Delta\braket{p_{i}^{\alpha}q_{j}^{\beta}}$ as

[TABLE]

where

[TABLE]

We first approximate

[TABLE]

For the Hamiltonian term we obtain

[TABLE]

B.6. Kernel approximation

One of the approximations we did in the previous derivation of the moment equation is to replace the expectation of a kernel by the kernel of the expected values.

If we expand the kernel in powers of its argument, we could compute the errors up to any order. For example for a Gaussian kernel, the first few terms are

[TABLE]

The main problem with this approximation with polynomials is that the approximation to any order corresponds to having a kernel with unbounded values for large arguments. This results in non-physical and large interactions of particles far away, which should normally not interact. Obtaining a reliable expansion of a kernel function is thus a difficult task in the moment approximation.

Nevertheless, one could consider the following higher order approximation of the expected value of a Gaussian kernel

[TABLE]

where the function $\theta(\mathbf{x})$ is given by

[TABLE]

and the coefficient $f_{\alpha}$ is found such that this approximation is the best fit to the Gaussian. In practice, we have $f_{\alpha}\approx 0.6$, but this value depends on $\alpha$ in general. This cutoff function $\theta$ is necessary here, otherwise, this approximation will not be bounded, leading to large errors in the dynamics. The expected value of this approximation assumes that the $\theta$ function commutes with it, and only takes into account the approximation of the quadratic term. It turned out that for all our experiments, these correction terms did not substantially improve the result, thus we did not include them in the equations.