Simulation of Conditioned Diffusions on the Flat Torus

TL;DR

This paper introduces a new method for simulating conditioned diffusion processes on the flat torus, with proven convergence properties and a measure change technique that relates Euclidean diffusions to Brownian motion.

Contribution

It presents a novel approach for simulating conditioned diffusions on the flat torus using projection from Euclidean space, with convergence analysis and measure change insights.

Findings

Convergence of the projected diffusion process to the terminal point on the torus.

A change of measure makes the Euclidean diffusion locally resemble Brownian motion.

Method applicable to biological data involving toroidal structures.

Abstract

Diffusion processes are fundamental in modelling stochastic dynamics in natural sciences. Recently, simulating such processes on complicated geometries has found applications for example in biology, where toroidal data arises naturally when studying the backbone of protein sequences, creating a demand for efficient sampling methods. In this paper, we propose a method for simulating diffusions on the flat torus, conditioned on hitting a terminal point after a fixed time, by considering a diffusion process in R 2 which we project onto the torus. We contribute a convergence result for this diffusion process, translating into convergence of the projected process to the terminal point on the torus. We also show that under a suitable change of measure, the Euclidean diffusion is locally a Brownian motion.

Equations

$$d{X}_t=\frac{b-X_t}{T-t}dt+d{W}_t;0\le t<T\text{and}{X}_0=a,$$

$$d{X}_t=b(t,X_t)dt+\sigma (t,X_t)d{W}_t,X_0=a,$$

$$d{Y}_t$$

$$\tilde{b}(t,x)$$

$$\displaystyle p(s,x;t,y)=\frac{1}{\sqrt{2\pi(t-s)}}\exp\biggl{(}-\frac{||x-y||^{2}}{2(t-s)}\biggr{)},\qquad s<t,$$

$$d{X}_t=b(t,X_t)dt+\sigma d{W}_t;0\le t<T\text{ and }{X}_0=a\text{ a.s.},$$

$$h(W_t)\stackrel{m}{=}h(W_0)-\frac{1}{2}{\int }_0^t{\Delta }_{{\mathbb{R}}^2}h(W_s)ds,$$

$$\tilde{g}(W_t)$$

$$\{\omega \in \Omega :B_T(\omega )=a\}=\{\omega \in \Omega :W_T(\omega )\in {\pi }^{-1}(a)\},$$

$$h(t,z)=\sum _{y\in {\pi }^{-1}(a)}\frac{p(t,z;T,y)}{p(0,z_0;T,y)},$$

$$\begin{split}dZ_{t}&=\sigma^{2}\nabla_{z}\log\bigg{(}\sum_{y\in\pi^{-1}(a)}p(t,z;T,y)\bigg{)}\bigg{|}_{x=Z_{t}}dt+\sigma dW_{t}\\ &=\sum_{y\in\pi^{-1}(a)}g_{y}(t,Z_{t})\frac{y-Z_{t}}{T-t}dt+\sigma dW_{t},\qquad Z_{0}=z_{0},\end{split}$$

$$g_{y}(t,x)=\frac{\exp\big{(}-\frac{\lVert y-z\rVert^{2}}{2\sigma^{2}(T-t)}\big{)}}{\sum_{y\in\pi^{-1}(a)}\exp\big{(}-\frac{\lVert y-z\rVert^{2}}{2\sigma^{2}(T-t)}\big{)}}.$$

$$d{X}_t=1_{G^c}(X_t)\frac{\alpha (X_t)-X_t}{T-t}dt+\sigma d{W}_t,X_0=x_0$$

$$\alpha (X_t)$$

$$1_{G^c}(X_t)\frac{\alpha (X_t)-X_t}{T-t}=\{\begin{array}{ll}\frac{\alpha (X_t)-X_t}{T-t},& \text{if }{X}_t\in G^c\\ 0,& \text{otherwise, }\end{array}$$

$$G^{c}=\bigcup_{y\in\pi^{-1}(a)}\bigl{(}y_{1}-\tfrac{1}{2},y_{1}+\tfrac{1}{2}\bigr{)}\times\bigl{(}y_{2}-\tfrac{1}{2},y_{2}+\tfrac{1}{2}\bigr{)}=:\bigcup_{y\in\pi^{-1}(a)}V_{y}.$$

$$b^{-1}(U)=b^{-1}(U)\cap\big{(}[0,T)\times G^{c}\big{)}\cup b^{-1}(U)\cap\big{(}[0,T)\times G\big{)}.$$

$$b^{-1}(U)\cap ([0,T)\times G)=\{\begin{array}{ll}[0,T)\times G,& \text{if }(0,0)\in U\\ \varnothing ,& \text{otherwise, }\end{array}$$

$${\Vert 1_{G^c}(x)\frac{\alpha (x)-x}{T-t}\Vert }^2={\Vert \frac{y-x}{T-t}\Vert }^2\le \frac{C}{(T-S{)}^2}=C_S,$$

$$\mathcal{E}(L)_{t}=\exp\Bigl{(}-\int_{0}^{t}b(s,X_{s})dW_{s}-\frac{1}{2}\int_{0}^{t}\lVert b(s,X_{s})\rVert^{2}ds\Bigr{)},$$

$$\mathbb{E}\Bigl{[}\exp\Bigl{(}\int_{0}^{t}\lVert b(s,X_{s})\rVert^{2}ds\Bigr{)}\Bigr{]}\leq\exp(tC_{S})<\infty$$

$${\tilde{W}}_t=W_t+{\int }_0^{t}b(s,X_s)ds$$

$$x_{t_{i+1}}=x_{t_i}+\frac{{\mathrm{argmin}}_{y\in {\pi }^{-1}(a)}(\Vert y-x_{t_i}\Vert )-x_{t_i}}{T-t_i}\Delta t+\sigma \Delta W_{t_i},$$

Full text

1 Introduction

Stochastic differential equations are ubiquitous in models describing evolution of dynamical systems with, e.g. in modelling the evolution of DNA or protein structure, in pricing financial derivatives, or for modelling changes in landmark configurations which are essential in shape analysis and computational anatomy. In settings where the beginning and end values are known on some fixed time interval, the use of Brownian bridges becomes natural to evaluate the uncertainty on the intermediate time interval.

When the data elements are elements of non-linear spaces, here differentiable manifolds, methodology for simulating bridge processes is lacking. In particular, in cases where the transition probability densities are intractable, it is of interest to use simulation schemes that can numerically approximate the true densities. In this paper we propose a method for simulating diffusion bridges on the flat torus, $\mathbb{T}^{2}=\mathbb{R}^{2}/\mathbb{Z}^{2}$, i.e. we propose a process that can easily be simulated and satisfies that the distribution of the true bridge of interest is absolutely continuous with respect to the distribution of this proposal process. This specific case will serve as an example of the more general setting of simulating diffusion bridge processes on Riemannian manifolds. Because of the non-trivial topology of the torus $\mathbb{T}^{2}$, the conditioned process will be equivalent to a process in $\mathbb{R}^{2}$ that is conditioned on ending up in a set of points. Therefore, we will address the question of conditioning a process on infinitely many points. Secondly, we will handle the case when the process crosses the cut locus of the target point, i.e. the set of points with no unique distance minimizing geodesic.

It is a basic consequence of Doob’s h-transform that the distribution of a conditioned diffusion process is the same as another diffusion process with the drift depending on the transition density. However, as mentioned in [1], using this transform directly is undesirable for simulation purposes as the transition density is often intractable. Instead, the authors introduce a diffusion process which can easily be simulated and with the property that the distribution of the true conditioned diffusion is absolutely continuous wrt. the diffusion used for simulation. We here use this approach that in [1] covers the Euclidean case as the starting point for developing a simulation scheme on the torus.

@afterheading

Recent papers have considered diffusion processes on the torus, for example, Langevin diffusions on the torus were studied in [3] and [4], in the latter to describe protein evolution. In this paper, we introduce a diffusion process in $\mathbb{R}^{2}$ which can easily be simulated and projects onto a bridge process on the torus. More generally, Brownian bridges on manifolds have been studied for example in the context of landmark manifolds [9] and used for approximating the transition density of the Brownian motion. The present paper uses bridges on the flat torus to exemplify how some of the challenges of bridge simulation on Riemannian manifolds can be addressed, here in particular non-trivial topology of the manifold.

@afterheading

We begin in Section 2 with a short introduction to Brownian bridge processes in the standard Eucliden case and how it relates to the definition of a Brownian bridge process on the flat torus. At the end we introduce the stochastic differential equation (SDE) which will be used for simulating the bridge process. In Section 3 we argue that a strong solution of our proposed SDE exist. We show results about convergence and absolute continuity in Section 4. Numerical examples are presented in Section 5.

2 Theoretical Setup

This section will briefly review some Brownian bridges theory and discuss the torus case. A more general theory of diffusion bridges can be found in [1], constituting the main reference for this work. At the end, we introduce our proposal process.

@afterheading

Consider a Brownian motion $W=(W_{t})_{t\geq 0}$ in $\mathbb{R}^{n}$. By conditioning, it can be shown that $W$ will end up at a given point at a given time. For example, the process given by $B_{t}=W_{t}-\frac{t}{T}W_{T}$ defines a Brownian bridge conditioned to return to [math] at time $T$. It can be shown that the diffusion process given by

[TABLE]

for given $a,b\in\mathbb{R}^{d}$ and $W$ a $d$-dimensional standard Brownian motion, is a $d$-dimensional Brownian bridge from $a$ to $b$ on $[0,T]$ (see e.g. [6, sec. 5.6]). More generally, diffusion bridges can be defined through Doob’s $h$-transform, that is, the distribution of a diffusion

[TABLE]

conditioned on $X_{T}=b$ is the same as that of

[TABLE]

where $p(t,x;T,b)$ denotes the transition density of the process $X$. In the usual setting where $p$ is the transition density of a Brownian motion it has the form

[TABLE]

which yields (1).

We propose a method similar to the Euclidean scheme [1] for simulating Brownian bridges on the flat torus, which is of the form

[TABLE]

where $\sigma>0$, $a\in\mathbb{T}^{2}$ is given, and $W$ is a two-dimensional standard Brownian motion. The exact form of $b(t,x)$ will become apparent below. It is important here to note that in the particular case of the flat torus the transition density for the Brownian motion is known and therefore it is possible to simulate from the distribution of the true Brownian bridge on $\mathbb{T}^{2}$, however, it requires the calculation of the distance to infinitely many points which the proposed model does not. In Figure 4 is shown paths of the proposed model and the corresponding paths of the true bridge process.

@afterheading

Let $\pi\colon\mathbb{R}^{2}\rightarrow\mathbb{T}^{2}=\mathbb{R}^{2}/\mathbb{Z}^{2}$ denote the canonical projection onto the torus. The standard two-dimensional Brownian motion $W=(W^{1},W^{2})$, for two independent one-dimensional Brownian motions $W^{1}$ and $W^{2}$, is mapped to a Brownian motion $B=(B_{t})_{t\geq 0}$ on the flat torus $\mathbb{T}^{2}$ by the projection map $\pi$. Indeed, we can identify the torus $\mathbb{T}^{2}$ with the unit cube $Q=\{x\in\mathbb{R}^{2}:-\frac{1}{2}\leq x_{k}<\frac{1}{2},k=1,2\}$. Then for $g\in C^{\infty}(\mathbb{T}^{2})$ the Laplace-Beltrami operator, $\Delta_{\mathbb{T}^{2}}$, on $\mathbb{T}^{2}$ corresponds to the restriction to $Q$ of the usual Euclidean Laplacian, $\Delta_{\mathbb{R}^{2}}\tilde{g}$, where $\tilde{g}$ denotes the periodic extension of $g$, i.e. $\tilde{g}=g\circ\pi$ (see [8, Sec. 3.5]). Since $W$ is a Brownian motion in $\mathbb{R}^{2}$ if and only if it satisfies the diffusion equation

[TABLE]

for all smooth functions $h$, where $X\overset{m}{=}Y$ means that the difference $X-Y$ is a local martingale (see e.g. [2, Sec. 1.5]), it follows that, for $h=\tilde{g}$,

[TABLE]

As this holds for all smooth functions $g$ on $\mathbb{T}^{2}$, we get that $B$ is a Brownian motion on $\mathbb{T}^{2}$ in agreement with the definition of a manifold-valued Brownian motion given in [5, Sec. 3.2].

By conditioning $B$ on $\mathbb{T}^{2}$ to hit a given point $a\in\mathbb{T}^{2}$, at some fixed time $0\leq T<\infty$, it is seen that

[TABLE]

and so simulating a Brownian bridge on the flat torus $\mathbb{T}^{2}$ is equivalent to simulating a two-dimensional standard Brownian motion conditioned to end up in the set $\pi^{-1}(a)$ at time $T$. The diffusion given by (1) will not suffice as it is constructed to hit exactly one point. It will, however, provide one subset of sample paths of the Brownian bridge on $\mathbb{T}^{2}$, corresponding to subset of paths that will ”unwrap” the same number of times that it ”wraps” around the cut locus. This is illustrated in Figure 1. To give a precise meaning to this statement we consider the $h$-transform

[TABLE]

with $p$ denoting the transition density of the two-dimensional Brownian motion, which by Doob’s $h$-transform implies that the distribution of $W$ conditioned on $W_{T}\in\pi^{-1}(a)$ is the same as the distribution of the diffusion

[TABLE]

where

[TABLE]

Instead, we propose to consider the diffusion process on $[0,T)$, for some fixed positive $T$, defined by

[TABLE]

where $\sigma>0$ and $\alpha$ is defined by

[TABLE]

with $a\in\mathbb{T}^{2}$, and where $G$ is the set of ”straigt lines” of the form $\mathbb{R}\times\{x\}$ (resp. $\{x\}\times\mathbb{R}$) in $\mathbb{R}^{2}$ where $\alpha(X_{t})$ is not unique (see Figure 1). The indicator function removes the drift when the process does not have a natural attraction point.

3 Existence of Strong Solution

The drift term in equation (4) is discontinuous. However, we below show that it posses certain regularity conditions and use this to show that a strong solution to the SDE exist.

@afterheading

In order to ensure the existence of a solution to the diffusion in (4), we need some regularity of the drift term. The drift coefficient is given by

[TABLE]

for every $0\leq t<T$, where the superscript $c$ denotes the complement. It is a discontinuous process with the set of discontinuities being the set $G$ consisting of the set of straight lines in $\mathbb{R}^{2}$ where the argmin process is non-unique. It is not even clear that the drift term is suitably measurable as the argmin map in general is not.

Lemma 1.

Let $b\colon[0,T)\times\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ be the map given by (5). Then $b$ is \mathcal{B}([0,T))\otimes\mathcal{B}\big{(}\mathbb{R}^{2}\big{)}-\mathcal{B}\big{(}\mathbb{R}^{2}\big{)} measurable. Furthermore, the map $(s,\omega)\mapsto b(s,X_{s}(\omega))$ is $\mathcal{B}([0,t])\otimes\mathcal{F}^{0}_{t}$ measurable, for every $0\leq t<T$, where $(\mathcal{F}^{0}_{t})$ denotes the natural filtration generated by $X$. This is called progressive measurability.

Proof.

First note that $G^{c}$ is a Borel measurable set as we can write write it as a countable union of open sets, i.e., for $y=(y_{1},y_{2})$ we have

[TABLE]

Now, we need to show that for all A\in\mathcal{B}\big{(}\mathbb{R}^{2}\big{)}, the set $b^{-1}(A)$ is an element of \mathcal{B}([0,T))\otimes\mathcal{B}\big{(}\mathbb{R}^{2}\big{)}. It is enough to consider all open subsets $U\subseteq\mathbb{R}^{2}$ as these sets generate the Borel algebra on $\mathbb{R}^{2}$. So let $U$ be an arbitrary open subset, then we have that

[TABLE]

As $b$ is continuous on each of the sets $[0,T)\times V_{y}$ we have that b^{-1}(U)\cap\big{(}[0,T)\times G^{c}\big{)} is a countable union of open sets and therefore an element of \mathcal{B}([0,T))\otimes\mathcal{B}\big{(}\mathbb{R}^{2}\big{)}. For the second part we see that

[TABLE]

where both are elements of \mathcal{B}([0,T))\otimes\mathcal{B}\big{(}\mathbb{R}^{2}\big{)}. This shows that $b$ is Borel measurable.

@afterheading

Progressive measurability follows by a very similar argument. ∎

Usually, global or local Lipschitz conditions are imposed on the drift and diffusion coefficients in order to secure global (resp. local) strong solutions to an SDE. This is a too strong condition for the drift term in this case, however, it is bounded in the following sense.

Lemma 2.

The drift coefficient in (5) is uniformly bounded in $x$ and in $t$ on $[0,S]$, for any $0\leq S<T$.

Proof.

The first assertion is clear. Let $S\in[0,T)$ be arbitrary and $0\leq t\leq S$. For every $x\in G^{c}$ there exist a $y\in\pi^{-1}(a)$ such that we have

[TABLE]

for some positive constants $C>0$. ∎

We now come to the main result of this section.

Proposition 1.

There exist a strong solution of (4) on $[0,T)$, which is strongly unique.

Proof.

The drift term is Borel measurable and bounded on $[0,S]$ by Lemma 1 and Lemma 2. As indicated in [10, Thm. 2] and [11, Thm.1] (4) has a strong solution which is strongly unique. ∎

Remark 1.

The assumption in [10, Thm. 2] can be verified by using smooth bump functions.

4 Convergence and Absolute Continuity

The considerations above make the solution of (4) into a continuous semimartingale. If a semimartingale $X$ takes its values in an open set $U$ of $\mathbb{R}^{2}$ then Itô’s formula holds true for any $C^{1,2}([0,T)\times U)$ functions as well.

Proposition 2.

Let $X$ be a solution to (4) on the filtered probability space $(\Omega,\mathcal{F},(\mathcal{F}_{t}),P)$. For every $\omega\in\Omega$ for which there exist an $S<T$ such that $X_{t}(\omega)$ stays in $G^{c}$ on $[S,T)$, then $X$ converges pointwise almost surely to $\pi^{-1}(a)$.

Proof.

Assume that for some $\omega\in\Omega$ there exist some $S<T$ such that on $[S,T)$ the process $X_{t}(\omega)$ takes its values in $G^{c}$. By continuity of the process it will take it its values in some open neighborhood $V_{y}$ of the point $y\in\pi^{-1}$. The proof is then identical to the proof in [1, Lemma 4]. ∎

Remark 2.

It is of course of interest to show that for almost every path the process will converge. This can be obtained by showing that the process will not intersect $G$ infinitely many times close to $T$.

Consider the stochastic process $\mathcal{E}$ on $0\leq t\leq S$ defined by

[TABLE]

where $L$ is the local martingale in the exponential. This is known as the Doléans-Dade exponential. From Lemma 2 it follows that, for all $t\leq S$,

[TABLE]

The above is known as the Novikov condition (cf. [7]) which ensures that (6) is a martingale on $[0,T)$. Girsanov’s theorem ([6, Thm. 5.1 Chap. 3]) then provides that the process defined by

[TABLE]

is a Brownian motion under the new measure $Q$ introduced below.

Theorem 4.1.

Let $X$ defined on $(\Omega,\mathcal{F},(\mathcal{F}_{t}),P)$ be a solution of (4) on $[0,S]$ for $S<T$. The process in (6) defined on $0\leq t\leq S$ ($S<T$) is a true martingale and so there exists a measure $Q$ which is absolutely continuous wrt. $P$ such that $X$ is $Q$-Brownian motion.

Proof.

The martingale property of (6) on $[0,S]$ is a consequence of the Novikov condition. Then Girsanov’s theorem gives us that $X$ is a $Q$-Brownian motion on $[0,S]$. ∎

From the (perhaps obvious) fact that the distribution of the true Brownian bridge is locally equivalent to the distribution of the Brownian motion up to time $t<T$, it follows that the distribution of the Brownian bridge is absolutely continuous wrt. the proposed process up to time $t<T$.

Remark 3.

A bit of extra work is needed to obtain the correction term as in [1]. There are indications that it is possible to simulate from the true distribution of the Brownian bridge on the torus, however, Theorem 1 shows that (4) can approximate it.

5 Numerical Experiments.

For the numerical implementation of the proposed SDE in equation (4) we implemented the Euler-Maruyama scheme, i.e. taking $n$ equidistant discretization points of the time interval $t_{1},...,t_{n}$, with $t_{i+1}-t_{i}=\Delta t$, the numerical equation becomes

[TABLE]

where $\Delta W_{t_{i+1}}=W_{t_{i+1}}-W_{t_{i}}$ is equal in distribution to a normal random variable with mean zero and variance $\Delta t$.

Figure 2(a) shows the implementation of the numerical scheme on an embedded torus and Figure 2(b) its Euclidean counterpart. Figure 3(a) shows the behaviour of the drift term along a given path, illustrating that the attraction becomes stronger as time approaches the terminal time. The vector fields in Figure 3(b) shows the constant attraction to the center of the open subsets.

Acknowledgements. We acknowledge F. van der Meulen for discussions and insights on conditioned diffusions.