# Simulating Infinite-dimensional Nonlinear Diffusion Bridges

Gefan Yang1 Elizabeth Louise Baker1 Michael L. Severinsen2
Christy Anna Hipsley3 Stefan Sommer1
1Department of Computer Science, University of Copenhagen
2Globe Institute, University of Copenhagen
3Department of Biology, University of Copenhagen
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

The diffusion bridge is a type of diffusion process that conditions on hitting
a specific state within a finite time period. It has broad applications in
fields such as Bayesian inference, financial mathematics, control theory, and
shape analysis. However, simulating the diffusion bridge for natural data can
be challenging due to both the intractability of the drift term and continuous
representations of the data. Although several methods are available to
simulate finite-dimensional diffusion bridges, infinite-dimensional cases
remain unresolved. In the paper, we present a solution to this problem by
merging score-matching techniques with operator learning, enabling a direct
approach to score-matching for the infinite-dimensional bridge. We construct
the score to be discretization invariant, which is natural given the
underlying spatially continuous process. We conduct a series of experiments,
ranging from synthetic examples with closed-form solutions to the stochastic
nonlinear evolution of real-world biological shape data, and our method
demonstrates high efficacy, particularly due to its ability to adapt to any
resolution without extra training.

##  1 Introduction

Diffusion processes are commonly utilized in diverse scientific fields,
including mathematics, physics, evolutionary biology, and finance, to model
stochastic dynamics. Updating the posterior of the model through conditioning
on existing observed data is a crucial element of the process, and there are
several methods to achieve this. For example, the Gaussian process regression
can be used for finite-dimensional or infinite-dimensional linear models [Shi
and Choi, 2011], and finite-dimensional Doob’s $h$-transform is designed for
finite nonlinear processes [Rogers and Williams, 2000, Chapter 6]. On the
other hand, many data types, such as images and sound signals, are naturally
continuous and described by functions, and the stochasticity is also
nonlinear. One approach is to discretize and represent them by vectors using
finite-dimensional models. However, a more natural way is to work directly in
the infinite-dimensional function space, performing only finite-dimensional
projections during inference. This functional formulation offers resolution-
invariant features, which are more memory-efficient and demonstrate greater
generalization capacities. Simulating infinite-dimensional linear processes
has been studied in some research [Franzese et al., 2024, Lim et al., 2023,
Pidstrigach et al., 2023], while the extension to the conditioning of
nonlinear processes remains an open question.

Interest in the study of linear diffusion processes has surged recently,
largely driven by the advancements in diffusion generative models. In
diffusion generative models, the data undergoes perturbation through a forward
unconditional linear diffusion process, resulting in noise that follows a
simple distribution and is easy to sample. The sampled noise is then
transformed back to the clean data through a reverse process, which is also
linear and can be learned by the score-matching techniques [Hyvärinen and
Dayan, 2005, Vincent, 2011]. Recently, the study of diffusion models using
stochastic differential equations (SDEs) and their mathematical
interpretations has gained significant attention [Song et al., 2020, Huang et
al., 2021]. Based on the rich research on the time-reversed diffusion
processes and the design of neural network approximation, we show that similar
techniques can be applied to simulate the conditioning of nonlinear diffusion
processes.

Figure 1: The nonlinear stochastic bridge between two continuous butterfly
shapes (blue and red) is characterized by discrete landmarks. The nonlinearity
of the bridge prevents the shape from overlapping and the landmarks from
colliding, which could happen with a linear process.

In order to condition an infinite-dimensional nonlinear diffusion process on
function-valued observations, we use the theory of infinite-dimensional Doob’s
$h$-transform, as outlined in Baker et al. [2024]. This transformation
converts an unconditional SDE into a conditional one, which is commonly
referred to as the diffusion bridge [Delyon and Hu, 2006, Schauer et al.,
2017, Heng et al., 2021]. Theoretically, the existence of such conditional
processes has been established, and our research has focused on developing a
numerically efficient method for their simulation. This involved reversing the
time of the diffusion bridge, which adheres to an SDE characterized by a
learnable drift. We also propose a tractable loss function for learning such
drift, which is inspired by the denoising score-matching object [Vincent,
2011]. By directly tackling the learning of the infinite-dimensional process,
we employ a time-dependent neural operator as the approximator. Given the
temporal symmetry between the diffusion bridge and its reversal, we simulate
the reversed bridge backward in time, effectively recreating paths of the
original bridge and thereby resolving the issue.

We aim to incorporate the score-matching techniques into the diffusion bridge
simulation methodology, especially to handle infinite-dimensional nonlinear
bridges. Our contributions are outlined as follows:

  * •

Our study introduces a novel approach for simulating infinite-dimensional
nonlinear diffusion bridges utilizing a score-matching technique. Expanding
upon the research of [Heng et al., 2021] that concentrated on simulating
bridges in finite dimensions, we further develop this technique to encompass
infinite dimensions. By directly approximating the additional drift introduced
by the infinite-dimensional Doob’s $h$-transform, we are able to reverse the
diffusion bridge and simulate it backward using the learned drift operator.
This framework holds the potential for addressing simulating general infinite-
dimensional nonlinear diffusion bridges, especially those that involve
continuous data within infinite-dimensional function spaces.

  * •

In order to approximate the drift operator in the reversed bridge, we have
devised a continuous time-dependent neural operator architecture that draws
inspiration from the Fourier time modulation [Park et al., 2023] and the
U-shaped neural operator [Ashiqur Rahman et al., 2022]. This structure is able
to represent continuous time-dependent nonlinear operators that are shown to
be discretization invariant. Furthermore, by splitting the large nonlinear
transition step into small tractable Gaussian steps, we design the loss
function on the Gaussian steps and, therefore, holds a computable form.

  * •

We illustrate our proposed methodology using continuous function-data-valued
stochastic processes. This will involve a demonstration of the Brownian bridge
with a closed-form drift, as well as the simulation of stochastic nonlinear
biological shape evolution. To our best knowledge, this is the first time the
nonlinear diffusion shape bridges can be simulated using a discretization-
invariant framework.

##  2 Related work

Bridge simulations: The linear diffusion process always has a Gaussian
transition density, which gives the additional drift term in the corresponding
conditional process a closed form. The challenge lies in the nonlinear case,
where the transition is no longer Gaussian. Several methods have been
developed [Delyon and Hu, 2006, Van der Meulen and Schauer, 2020, Corstanje et
al., 2024] for simulating bridges in different scenarios. The essential idea
is to approximate the unknown drift with a proposed known one. The bridge,
including such an approximate drift, is called the proposed process. By
sampling from the proposed process and computing the likelihood ratio between
the proposed and the true ones, the MCMC is used to correct the sampling. Such
methods suffer from expensive matrix inversion in high-dimensional cases and
time-costly MCMC updating iterations. Another approach is to directly
approximate the drift, which is proposed by [Heng et al., 2021], where they
use denoising score matching to learn the score for the finite-dimensional
nonlinear bridge. Based on [Heng et al., 2021]’s work, we further generalize
it to the infinite-dimensional settings.

Infinite-dimensional diffusion models: (Score-based) diffusion generative
models (SGMs) are initially developed to generative samples in Euclidean
spaces [Sohl-Dickstein et al., 2015, Song and Ermon, 2019, Song et al., 2021].
Various studies have been done to generalize SGMs into infinite-dimensional
space with either Hilbert-space-defined score functions [Pidstrigach et al.,
2023, Baldassari et al., 2024, Lim et al., 2023] or finite-dimensional score
projections [Franzese et al., 2024, Hagemann et al., 2023]. However, due to
the generative purpose, their methods are all limited to the linear SDEs
(often referred to as variance preserving/exploring SDEs), i.e., with a
diffusion term independent of the value of the process, and such SDEs are
unconditional to the end states. While we consider a more general conditional
nonlinear case, where the drift and diffusion terms are more complex to derive
a closed-form solution.

Diffusion Schrödinger bridge: Another diffusion bridge related model is the
diffusion Schrödinger bridge [De Bortoli et al., 2021, Shi et al., 2024, Tang
et al., 2024, Thornton et al., 2022], where the transition between two
distributions are treated as the optimal transport (OT) problem. Among them,
the work in [Shi et al., 2024] is closely related to ours, as they also
minimize the KL divergence to match the unknown bridge. However, since their
focus is on generative purposes, a simpler finite-dimensional linear process
and the bridge points are considered. We are instead interested in the
imputation between the points.

##  3 Infinite-dimensional nonlinear bridges

###  3.1 Conditioning via infinite-dimensional Doob’s $h$-transform

The infinite-dimensional diffusion bridge can be achieved by applying the
infinite-dimensional Doob’s $h$-transform [Baker et al., 2024] on an
unconditional diffusion process. Let $(H,\langle\cdot,\cdot\rangle)$ denote a
separable Hilbert space with a chosen countable orthonormal base
$\\{e_{i}\\}_{i=1}^{\infty}$, and let $(\Omega,\mathcal{F},\mathbb{P})$ be a
probability space with natural filtration $\\{\mathcal{F}_{t}\\}$, let
$\mathcal{B}(H)$ the Borel algebra, $\mathbb{P}(H)$ be the set of probability
measures on $(H,\mathcal{B}(H))$ and $\mu_{0}$, $\mu_{T}$ be two measures in
$\mathbb{P}(H)$ that we want to bridge between. A $H$-valued diffusion process
$X=(X(t))_{t\in[0,T]}$ is defined as:

$\mathop{\mathrm{d}X(t)}=f(t,X(t))\mathop{\mathrm{d}t}+g(t,X(t))\mathop{\mathrm{d}W(t)},\quad
X_{0}\sim\mu_{0}$ (1)

where $W(t)$ is a Wiener process in another Hilbert space $U$ (where $U$ can
equal $H$) with a trace-class covariance operator $C$ (i.e.,
$\operatorname{Tr}(C)<\infty$), $f:[0,T]\times H\to H$, $g:[0,T]\times
H\to\operatorname{HS}(C^{1/2}(U),H)$, where $\operatorname{HS}(C^{1/2}(U),H)$
denotes the Hilbert-Schmidt operator from $C^{1/2}(U)$ to $H$. We can then
condition (1) on hitting $X(T)\sim\mu_{T}$ within a finite time $T$ by
applying infinite-dimensional Doob’s $h$-transform, such transformed diffusion
process $X^{c}=(X^{c}(t))_{t\in[0,T]}$ follows the SDE:

$\mathop{\mathrm{d}X^{c}(t)}=\\{f(t,X^{c}(t))+a(t,X^{c}(t))\nabla\log
h(t,X^{c}(t))\\}\mathop{\mathrm{d}t}+g(t,X^{c}(t))\mathop{\mathrm{d}W(t)}$ (2)

with $X^{c}(0)=x_{0}\sim\mu_{0},\;X^{c}(T)=x_{T}\sim\mu_{T}$, where
$a:=gg^{*}:[0,T]\times H\to\operatorname{HS}(C^{1/2}(U),H)$, $\nabla\log
h:[0,T]\times H\to C^{1/2}(U)$ is another Hilbert-Schmidt operator, we shall
call it “score operator” for short in the following text. $h:[0,T]\times H\to
H$ is the Doob’s $h$ function, defined as
$h(t,x):=\mathbb{P}(X(T)\sim\mu_{T}\mid X(t)=x)$. The composition
$a(\cdot,\cdot)\nabla\log h(\cdot,\cdot):[0,T]\times H\to H$ is then an
operator that maps within $H$ and, therefore, can be treated as the additional
drift term.

In finite dimension, the $h:[0,T]\times\mathbb{R}^{d}\to(0,\infty)$ can be
chosen as a function satisfying $h(t,x)=\int h(t+s,y)p_{(t+s)|t}(y\mid
x)\mathop{\mathrm{d}y}$ [Rogers and Williams, 2000, Chapter 6], where
$p_{(t+s)|t}(y\mid x)$ is the conditional probability density with respect to
the Lebesgue measure $\mathop{\mathrm{d}y}$ in $\mathbb{R}^{d}$. A common
choice of $h$ is $h(t,x):=p_{T|t}(y\mid x)/p_{T|0}(y\mid x_{0})$, and
therefore the score becomes $\nabla\log p_{T|t}(y\mid x_{t})$. However, since
the Lebesgue measure $dy$ and the probability density $p$ are no longer well-
defined in $H$, such a choice is unavailable for infinite dimension. The
definition of $\nabla\log h$ as an operator, as described above, can be used
to address the well-posedness of the score. However, simulating (2) is still
nontrivial even with the finite-dimensional projection as used in [Franzese et
al., 2024] because of the intractable conditional probability, which only has
a closed form when $g$ is $X$-independent. To simulate a more general
nonlinear bridge, inspired by [Heng et al., 2021], we apply the generalized
time reversal theorem on the bridge process, and it turns out that in the
time-reversed bridge process, the additional drift term can be approximated by
sampling from the unconditional process (1). In the following sections, we
shall detail the infinite-dimensional time reversal and a learnable object.

###  3.2 Time-reversed bridges

Under specific conditions, the finite-dimensional version of (1)
$X^{d}(t)\in\mathbb{R}^{d}$ has a time-reversed process
$Y^{d}(t)\in\mathbb{R}^{d}$ that satisfies [Haussmann and Pardoux, 1986]:

$\displaystyle\text{forward:}\quad\mathop{\mathrm{d}X^{d}(t)}=f^{d}(t,X^{d}(t))\mathop{\mathrm{d}t}+g^{d}(t,X^{d}(t))\mathop{\mathrm{d}W^{d}(t)},\quad
X^{d}(t)\in\mathbb{R}^{l},$ (3)
$\displaystyle\text{reversed:}\quad\mathop{\mathrm{d}Y^{d}(t)}=\bar{f}^{d}(t,Y^{d}(t))\mathop{\mathrm{d}t}+g^{d}(T-t,Y^{d}(t))\mathop{\mathrm{d}B^{d}(t)},\quad
Y^{d}(t)\in\mathbb{R}^{l},$ (4)

with
$f^{d}:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d},g^{d}:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{l}$
as finite-dimensional drift and diffusion terms, $W^{l}(t),B^{l}(t)$ are
Wiener processes in $\mathbb{R}^{l}$,
$\bar{f}^{d}(t,x)=-f^{d}(T-t,x)+\frac{1}{p_{T-t}(x)}\nabla[a^{d}(T-t,x)p_{T-t}(x)]$,
and $p_{t}(x)$ as the probability density in $\mathbb{R}^{d}$. Furthermore,
[Millet et al., 1989] shows that such a time reversal exists for infinite-
dimensional SDEs. Suppose the Hilbert space $U$, where the Wiener process
$W(t)$ lives, has a countable orthonormal basis $\\{k_{j}\\}_{j=1}^{\infty}$,
with the basis $\\{e_{i}\\}_{i=1}^{\infty}$ of $H$, we can write (1) in terms
of these bases:

$\mathop{\mathrm{d}X_{i}(t)}=f_{i}(t,X(t))\mathop{\mathrm{d}t}+\sum_{j}^{\infty}g_{ij}(t,X(t))\mathop{\mathrm{d}W_{j}(t)},\quad
i=1,\dots,\infty,$ (5)

where $f_{i}:[0,T]\times H\to\mathbb{R}$ is defined as $\langle
f(t,x),e_{i}\rangle$ and $g_{i,j}:[0,T]\times H\to\mathbb{R}$ is defined by
$\langle g(t,x)(e_{j}),e_{i}\rangle$.

This has a family of time-reversed processes:

$\mathop{\mathrm{d}Y_{i}(t)}=\bar{f}_{i}(T-t,Y(t))\mathop{\mathrm{d}t}+\sum^{\infty}_{j}g_{ij}(T-t,Y(t))\mathop{\mathrm{d}B_{j}(t)},\quad
i=1,\dots,\infty,$ (6)

where
$\bar{f}_{i}(t,x)=-f_{i}(t,x)+\frac{1}{p_{T-t}(x_{i}\mid\xi_{i})}\sum_{j\in
I(i)}\nabla_{j}[a_{ij}(T-t,x)p_{T-t}(x_{i}\mid\xi_{i})]$, where
$p_{t}(x_{i}\mid\xi_{i})$ is the conditional probability density of
$X^{i}(t)=(X_{j}(t),j\in I(i))$ given $(X_{j}(t),j\notin I(i))=\xi_{i}$, and
$I(i)$ is a finite set of indices characterized by $i$. Such a density is
assumed to exist with respect to the Lebesgue measure in $\mathbb{R}^{I}$ (See
[Millet et al., 1989, Section 5]), which is a finite subspace of $H$. We can
apply this time reversal upon the infinite-dimensional bridge (2) derived in
the previous section and obtain the infinite-dimensional time-reversed bridge.
We shall see that within this reversed bridge, the additional term does not
require the information of $\mu_{T}$, which makes it possible to learn from
the process that starts from $\mu_{0}$ without hitting $\mu_{T}$.

######  Theorem 3.1.

Let $f,g,$ and $a$ be as defined before, and let $B(t)$ be a Wiener process in
$U$. The conditional process $X^{c}(t)$ has a time-reversal
$Y^{c}=(Y^{c}(t))_{t\in[0,T]}:=X^{c}(T-t)$ following the SDE:

$\mathop{\mathrm{d}Y^{c}(t)}=\bar{f}(t,Y^{c}(t))\mathop{\mathrm{d}t}+g(T-t,Y^{c}(t))\mathop{\mathrm{d}B(t)}$
(7)

where $\bar{f}(t,x)=-f(T-t,x)+a(T-t,x)\nabla\log\bar{h}(T-t,x)+\nabla\cdot
a(T-t,x)$, and $\bar{h}(t,x)=\mathbb{P}(X(t)\in\mathop{\mathrm{d}x}\mid
X(0)=x_{0})=:P(0,x_{0},t,\mathop{\mathrm{d}x})$, which is the law of $X(t)$
given $X(0)=x_{0}$, and $\nabla\cdot a(t,x):=\sum_{j}\nabla_{j}a_{ij}(t,x)$.

The proof of Theorem 3.1 can be found in A.1. Note that (7) does not require
$\mu_{T}$. If the process starts from $\mu_{0}$, it is impossible to sample
from $\mu_{T}$ conditioned on the intermediate distribution $\mu_{t}$ without
constructing the bridge forehead. Therefore the $h$ in (2) is actually
inaccessible. In contrast, sampling from $\mu_{t}$ while conditioning on
$\mu_{0}$ is feasible, which can be achieved by sampling the unconditional
process that starts from $\mu_{0}$. This insight makes learning $\bar{h}$
possible instead. However, the general transition from $0$ to $t$ is still
intractable unless the original process has an $X$-independent diffusion term,
and then the transition shall be Gaussian, but we can address this issue by
using the Markovian property of the process, which enable us to split the
large transition into small steps with Gaussian-like transitions.

##  4 Learning the bridge

###  4.1 Learning the nonlinear drift

The idea of learning (7) is to minimize the distance between the estimated
bridge and the true bridge in the sense of KL divergence, which is sometimes
referred to as “bridge matching” [Shi et al., 2024, Liu et al., 2023a,
Delbracio and Milanfar, 2023]. We can write the object function as the KL
divergence between the estimated and the true bridge. However, the true bridge
still has no closed form, and the object function is thereby uncomputable.
Fortunately, the bridge is Markovian, which means any single transition step
is independent of the history, or we can say the whole information of the
process is “compressed” into a single small step. This property enables us to
learn with only small local steps and average it over the whole trajectory to
gain global information on the whole process.

Let $\mathbb{P}^{x_{0}}$ denote the path measure induced by (7) starting from
$Y^{c}(T)=x_{0}$, and $\mathbb{Q}_{\theta}^{x_{0}}$ be the path measure
induced by another time-reversed process $Y=(Y(t))_{t\in[0,T]}$ follows:

$\mathop{\mathrm{d}Y(t)}=f_{\theta}(t,Y(t))\mathop{\mathrm{d}t}+g(T-t,Y(t))\mathop{\mathrm{d}B_{t}}$
(8)

where $f_{\theta}(t,x)=-f(T-t,x)+\mathcal{G}_{\theta}(T-t,x)+\nabla\cdot
a(T-t,x)$, $\mathcal{G}_{\theta}:[0,T]\times H\to H$ is a parameterized
operator by $\theta$. Let
$\mathop{\mathrm{d}\mathbb{P}^{x_{0}}}/\mathop{\mathrm{d}\mathbb{Q}^{x_{0}}_{\theta}}$
denote the Radon-Nikodym derivative of $\mathbb{P}^{x_{0}}$ with respect to
$\mathbb{Q}^{x_{0}}_{\theta}$, then by minimizing the KL divergence
$D_{KL}(\mathbb{P}^{x_{0}}\mid\mathbb{Q}^{x_{0}}_{\theta})=\mathbb{E}_{X(t)}[\log\mathop{\mathrm{d}\mathbb{P}^{x_{0}}}/\mathop{\mathrm{d}\mathbb{Q}^{x_{0}}_{\theta}(X(t))}]$,
we shall approximate the bridge (7) using $\mathcal{G}_{\theta}$.

######  Theorem 4.1.

For any time partition $(t_{i})_{i=1}^{n}$ of the interval $[0,T]$, the loss
function which is defined as:

$\hat{L}(\theta):=\frac{1}{2}\sum^{n}_{i=1}\int^{t_{i}}_{t_{i-1}}\mathbb{E}_{X(t)}\left[\|A_{t}(x)\\{\mathcal{G}_{\theta}(t,x)-b(t,x,t_{i-1},x_{t_{i-1}})\\}\|_{U_{0}}\right]\mathop{\mathrm{d}t}$
(9)

is equivalent to $D_{KL}(\mathbb{P}^{x_{0}}\mid\mathbb{Q}^{x_{0}}_{\theta})$
up to a $\theta$-independent $C$, where
$b(t,x,t_{i-1},x_{t_{i-1}}):=a(t,x_{t})\nabla\log\mathbb{P}^{x_{0}}(X(t)\in\mathop{\mathrm{d}x}\mid
X(t_{i-1})=x_{t_{i-1}})$, $A_{t}(x)=C^{-1/2}g^{-1}(t,x):[0,T]\times H\to U$ is
a $\theta$-independent operator, $\|\cdot\|_{U}$ is the norm of $U$.

Furthermore, $A_{t}(x)$ is positive-definite because $g$ and $C$ are assumed
to be positive, so $\ker(A_{t}(x))=\\{0\\}$, and we can drop it when using
gradient descent as suggested in [Lim et al., 2023]. Finally, we train with
the following object function:

$L(\theta):=\frac{1}{2}\sum^{n}_{i=1}\int^{t_{i}}_{t_{i-1}}\mathbb{E}_{X(t)}\left\\{\|\mathcal{G}_{\theta}(t,x)-b(t,x,t_{i-1},x_{t_{i-1}})\|^{2}_{U}\right\\}\mathop{\mathrm{d}t}.$
(10)

###  4.2 Time-dependent neural operators

As described in the previous section, $\mathcal{G}_{\theta}:[0,T]\times H\to
H$ is a time-dependent operator. We use a time-dependent neural operator to
approximate it. Neural operators are initially developed to solve partial
differential equations where the mapping from initial conditions to solutions
is treated as a nonlinear smooth operator [Li et al., 2020a, b, Lu et al.,
2021, Kovachki et al., 2023]. Such a learned operator is shown to be
independent of the mesh or discretization used during the training [Kovachki
et al., 2023], thereby enabling them to operate on the continuous domain. We
choose the Fourier neural operator (FNO) proposed in [Li et al., 2020a] for
two main reasons: 1) FNO is efficient due to the Fast Fourier Transform (FFT)
and easy to implement; 2) Within the FNO, the functions are decomposed by the
Fourier basis, which can represent the continuous periodic signals. It is
especially suitable for our use since we model closed shapes that can be
treated as periodic.

Since our target operator $\mathcal{G}_{\theta}$ is time-dependent, we
incorporate the time information by using the continuous-in-time FNO (CTFNO)
proposed by [Park et al., 2023], which modulates the time embedding into the
Fourier modes. Besides the time step modulation, we also refer to the design
of UNO [Ashiqur Rahman et al., 2022], which is shown to outperform the
classical FNO structure used in [Li et al., 2020a] with less memory and higher
accuracy. By combining both, we propose the structure of U-shaped CTFNO, as
shown in Figure 2.

Figure 2: Continuous-time U-shaped FNO architecture, acting as a parameterized
operator $\mathcal{G}_{\theta}:(u,t)\mapsto v$, where both $u$ and $v$ are
evaluated on a discretized spatial-temporal grid $(\mathbf{x},\mathbf{t})$
with $\mathbf{x}\in\mathbb{R}^{m},\mathbf{t}\in\mathbb{R}^{n}$

In practice, we first specify the discrete spatial-temporal grid on which the
functions are evaluated. Normally, the grid resolution for training is lower
than that for inference. We assume for both input and output functions, the
dimension of the domains is 1, and the codomains are $d_{u}$ and $d_{v}$,
respectively. Then, the function evaluated at $m$-discretized domains can be
represented with the $m\times d_{u}$ matrix. In fact, in the context of bridge
approximation, $d_{u}=d_{v}$, we keep the notation for generality. After
discretizing the spatial-temporal grid as an
$\mathbb{R}^{m}\times\mathbb{R}^{n}$ space, we evaluate $u$ and obtain a
finite-dimensional vector $u(\mathbf{x},\mathbf{t})\in\mathbb{R}^{d_{v}\times
m\times n}$. Then $u(\mathbf{x},\mathbf{t})$ and $\mathbf{t}$ are fed into the
neural operator. Within the operator, $\mathcal{P}$ and $\mathcal{Q}$ are
lifting and projecting maps individually; the former increases the dimension
of the codomain, and the latter decreases. Both these two projections are
implemented as convolutional layers with 1-size kernels. At the same time, the
time t is transformed into two embeddings through $\varphi$ and $\psi$. We
implemented both as sinusoidal embeddings [Vaswani et al., 2017] as suggested
in [Park et al., 2023]. The lifted $v_{0}(\mathbf{x},\mathbf{t})$ and the
embeddings $\varphi(t),\psi(t)$ are fed into one of the continuous-time
Fourier layers $\mathcal{G}_{l}$. Within the block, the input $v_{0}$ is
transformed into the Fourier modes through the FFT; Meanwhile, the embedding
$\varphi(t)$ is transformed by a learnable linear transformation $A_{\theta}$
and will be used with a learnable linear transformation $R_{\theta}$ on the
Fourier modes through $A_{\theta}\varphi(t)R_{\theta}\mathcal{F}(v_{0})$.
Conversely, $\psi(t)$ and $v_{0}$ shall go through a similar process, except
all the actions are executed on the physical domain without Fourier
transformation. Note that the output spatial dimension does not have to be
equal to the input, i.e.,
$\mathbf{x_{out}}\in\mathbb{R}^{m_{out}},\mathbf{x_{in}}\in\mathbb{R}^{m_{in}},m_{out}\neq
m_{in}$; that is, a spatial subspace is focused, which is similar to the
classical U-net. Finally, the transformed Fourier component
$\mathcal{F}^{-1}[A_{\theta}\varphi(t)R_{\theta}\mathcal{F}[v_{0}]]$ and the
residual component $B_{\theta}\psi(t)W_{\theta}v_{0}$ are added and then
activated by a point-wise nonlinear activation function $\sigma$. Minic the
traditional U-net, our neural operator is equipped with a decreasing hierarchy
from top to bottom in two aspects: 1) the output spatial dimension is
decreasing; 2) the truncated Fourier modes in FFT are decreasing. Both the
designs are for extracting multi-level Fourier features [Ashiqur Rahman et
al., 2022]. The skip connections between the downsampling and upsampling paths
are kept as well. Compared to a normal FNO, our modified neural operator can
achieve better performance with fewer parameters and fast convergence speed,
thereby being suitable for the operator approximator in our case.

###  4.3 Numerical implementations

As we represent $\mathcal{G}_{\theta}$ as a neural operator, we implicitly
project the infinite-dimensional optimization object 10 into finite dimension.
Such projection $P^{d}:H\to\mathbb{R}^{d}$can be achieved by, for example,
taking the first $d$ components
$\\{e_{i}\\}^{d}_{i=1}\subset\\{e_{i}\\}^{\infty}_{i=1}$, and projecting the
element $h\in H$ by $h^{d}=\sum^{d}_{i=1}\langle h,e_{i}\rangle_{H}\cdot
e_{i}$. In the subspace $\mathbb{R}^{d}\subset H$, the conditional probability
$\mathbb{P}(X(t)\mid X(0)$) has a density function $p(X^{d}(t)\mid X^{d}(0))$
given the defined Lebesgue measure. $X^{d}(t)\in\mathbb{R}^{d}$ is a finite-
dimensional vector obtained by evaluating $X(t)$ on a finite-dimensional grid.
To sample $X^{d}(t)$, we need to solve the finite-dimensional version of (1):

$\mathop{\mathrm{d}X^{d}(t)}=f^{d}(t,X^{d}(t))\mathop{\mathrm{d}t}+g^{d}(t,X^{d}(t))\mathop{\mathrm{d}W^{k}(t)},\quad
X^{d}(0)\sim\mu^{d}_{0},$ (11)

where $f^{d}:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$,
$g^{d}:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d\times k}$, $W^{k}$ is a
standard Wiener process lives in $\mathbb{R}^{k}$, $\mu_{0}^{d}$ is the push-
forwards of $\mu_{0}$ by the projection map
$P^{d}:H\to\mathbb{R}^{d}$Pidstrigach et al. [2023], defined by:

$\mu^{d}:=(P^{d})_{*}\mu,\quad\text{where}\quad(P^{d})_{*}\mu(A)=\mu((P^{d})^{-1}(A)),A\subset
H.$ (12)

Solving (11) can be done using many finite-dimensional SDE numerical solving
approaches. We choose the Euler-Maruyama solver (See, for example, [Bayram et
al., 2018]). In the Euler-Maruyama scheme, the small-step transition is
modeled as a Gaussian density, and the time integral from $t_{i-1}$ to $t_{i}$
can be approximated by the discrete sum. The single-step update of Euler-
Maruyama is:

$X^{d}(t_{i})=X^{d}(t_{i-1})+f^{d}(t_{i-1},X^{d}(t_{i-1}))\delta
t+g^{d}(t_{i-1},X^{d}(t_{i-1}))(W^{k}(t_{i})-W^{k}(t_{i-1})),$ (13)

where $W^{k}(t_{i})-W^{k}(t_{i-1})\sim\mathcal{N}(0,\delta tI_{k})$ is the
independent Brownian increment. Then the approximate loss is:

$L^{d}(\theta):=\frac{1}{2}\delta
t\sum^{n}_{i=1}\mathbb{E}_{x_{0}}\\{\|\mathcal{G}_{\theta}(t_{i},X^{d}(t_{i}))-b^{d}(t_{i},X^{d}(t_{i}),t_{i-1},X^{d}(t_{i-1})\|^{2}\\}\\}.$
(14)

As $\delta t\to 0$, $L^{d}(\theta)\to L(\theta)$ by the property of Euler-
Maruyama solver. In practice, we choose the time steps as 0.01 within the
total time interval [0, 1], which we consider as a proper choice of balance
between accuracy and efficiency.
$b^{d}(t_{i},X^{d}(t_{i}),t_{i-1},X^{d}(t_{i-1}))=\Sigma^{d}(t_{i-1},X^{d}(t_{i-1}))\cdot\nabla\log
p(X^{d}(t_{i})\mid X^{d}(t_{i-1}))$.Since $p$ is Gaussian,
$b^{d}(t_{i},X^{d}(t_{i}),t_{i-1},X^{d}(t_{i-1}))$ can be easily computed as:

$b^{d}(t_{i},X^{d}(t_{i}),t_{i-1},X^{d}(t_{i-1}))=-\frac{1}{\delta
t}\left\\{X^{d}(t_{i})-X^{d}(t_{i-1})-\delta
tf^{d}(t_{i-1},X^{d}(t_{i-1}))\right\\}.$ (15)

Note that $b^{d}$ can be computed simultaneously as the SDE-solving step
iterates without additional computations afterward. So far, the parameterized
operator $\mathcal{G}_{\theta}$ and the computable loss function
$L^{d}(\theta)$ have been introduced. We summarize the algorithms for training
and inference as Algorithm 1 and 2 in the Appendix.

##  5 Experiments

###  5.1 Functional Brownian bridges

Quadratic functions: We first study the function-valued Brownian process,
whose score has a closed form when projected into the finite dimension. For
the object functions, we choose the quadratic function
$f(x)=ax^{2}+\varepsilon$, where $a=\\{1,-1\\}$ and
$\varepsilon\sim\mathcal{N}(0,10^{-4})$ is used to equip the set with non-zero
measure, as used in [Phillips et al., 2022]. We choose the function-valued
Brownian process as the unconditional process since its derived bridge has a
closed form. We choose the bounded grid of $[-1,1]$ to evaluate the function,
and the grid is uniformly discretized into a finite number of points, as shown
in Figure 3(a). Since the SDE is defined in $H$, such discretization stands
for choosing a finite base and projecting it into the finite dimension.
However, the projected process should be consistent under arbitrary
discretizations. We show it by training the neural operator under a low-
resolution scheme (8 points distributed within $[-1,1]$) and evaluate it in a
high-resolution case (128 points). Note that our learned bridge (Figure 3(b))
shows high consistency with the ground true Brownian bridge (Figure 3(c)),
even in a much finer grid. We also compared the MSE between the simulated
trajectories by the true score and the estimated one with the same random seed
initialization. The results indicate that the learned process is
discretization-invariant, i.e., it represents the infinite-dimensional
process. We further study the model’s performance under different
discretization levels, as shown in Figure 3(d); we observe that the model
shows great consistency under discretizations in terms of the MSE metrics.
More samples from the estimated bridge can be found in Figure 7

(a) Quadratic functions

(b) Learned bridge

(c) Ground truth

(d) MSE

Figure 3: Qualitative results for Brownian bridges between two quadratic
functions; (a) Two quadratic functions are chosen as the starting
($f(x)=x^{2}+\varepsilon$, blue) and the target ($f(x)=-x^{2}+\varepsilon$,
red). Both are evaluated at 8 points evenly distributed within $[-1,1]$ for
training (marked in cyan and orange respectively); (b) One sample from the
learned reversed Brownian bridge, evaluated at 128 evenly distributed points;
(c) One sample from the true reversed Brownian bridge, simulated with the same
random seed as (b) for comparison; (d) MSE between the estimated and the true
bridges under different resolutions.

Circles: We further test our method by simulating the Brownian bridge between
2D circles, which acts as the simplest 2D stochastic shape bridge and holds a
known form that we can compare. We shall see in the next experiment the true
shape bridge is inaccessible. The Circle can be treated as a function from
$\mathbb{R}$ to $\mathbb{R}^{2}$, and can be characterized by finite points
along the outlines. We are trying to bridge two circles with different radii.
Figure 4(a) gives a visualization of the forward unconditional process. Then
we train the model with only 16 points evenly sampled from the circle and
evaluate it in higher resolutions (See Figure 4(b) and 8). Also, Figure 4(c)
shows the ground true bridge simulated with the same random seed, and we could
expect consistency between these two simulations. Besides the visual
inspection, we show the score output at different stages of the process since
we have access to the ground true score that we can compare. The results are
shown in Figure 5. Since, in the Brownian case, there is no drift term in the
unconditional process, the additional drift introduced by the transformation
is the only force that drives the functions toward the target, which is
reflected in the visualization that the score arrows should always point
toward the center, despite of the discretization, which can be observed in the
Figure 5 as well.

(a) Unconditional process

(b) Learned bridge

(c) Ground truth

Figure 4: Qualitative results for Brownian bridges between two circles; (a)
The unconditional forward process starting from the red circle discretized
into 32 points, the training uses 16 points marked in orange crosses; (b) One
sample from the learned reversed Brownian bridge, evaluated at 128 evenly
distributed points; (c) One sample from the true reversed Brownian bridge,
simulated with the same random seed as (b) for comparison. Figure 5: The
inspections of the estimated score and the true score at different time steps
of the bridge process.

###  5.2 Biological shape evolution

We finally evaluate our method on real data. Specifically, we are interested
in modeling the stochastic change in butterfly morphometry over time. In the
phylogenetic analysis, such bridge processes model the transitions along the
edge between nodes in a phylogenetic tree, where the nodes store the observed
data that comes from specimens. Developing the fast bridge simulation approach
facilitates large-scale phylogenetic tree simulation and biological
inferences. However, when modeling the shape, the linear process can not be
directly applied since the topology of the shape must be preserved. We detail
our nonlinear setting for the shape process in B.4. By applying the proposed
method, we can now bridge between two butterfly shapes. It’s worth noticing
that the bridge is constructed in infinite rather than finite dimensions. To
prove this, we only train the score approximator with 32 points, corresponding
to a finite discretization of the shape function. During the inference, we
evaluate the bridge on different levels of discretization, and it shows that
the approximator can give reasonably enough interpolations between low-
discretized points, as shown in Figure 6. We point out that our method does
not constrain a specific shape since the deformation of the shape is modeled
instead of the shape itself. Also, since the bridge is constructed between two
sets with non-zero measures, in practice, it allows the observation to contain
noise. Both these features make our method suitable for modeling stochastic
biological shape evolution, where 1) Multiple similar shapes are considered,
and 2) The observations often contain noise. We test our model on more
butterfly species which are close on the phylogenetic tree, the results are
shown in Figure 9 and 10. This marks the first instance where the simulation
of nonlinear diffusion shape bridges has been achieved using a discretization-
invariant method, to the best of our knowledge.

(a) Evaluate on 64 landmarks

(b) Evaluate on 128 landmarks

(c) Evaluate on 256 landmarks

Figure 6: The stochastic shape bridge between two butterfly shapes (Papilio
polytes in red, Parnassius honrathi in blue). The shape is represented by
landmarks of different numbers. Our method only needs to train on a few sample
points (marked in cyan and orange crosses)and can generalize to arbitrary
resolutions.

##  6 Conclusion

We introduce a novel approach for simulating nonlinear diffusion bridges in
infinite-dimensional spaces. This method utilizes Doob’s $h$-transform applied
to an unconditional process, transforming it into a conditional diffusion
bridge. Subsequently, we reverse this bridge in time and simulate the reversed
sequence using a drift learned from the original, unconditional process.
Utilizing a neural operator, we establish that this learned drift retains
continuity across discretizations, affirming its role as a discretization-
invariant representation of the underlying infinite-dimensional dynamics. We
validate our approach through several functional diffusion bridge test cases,
especially the phylogenetic shape analysis problem, where our method showcases
its effectiveness of high-resolution zero-shot training and generalization to
various shapes, together with robustness under different levels of
discretizations.

However, simulating the time-reversed bridge poses computational challenges,
particularly the necessity of calculating the Jacobian of $a$. As outlined in
Heng et al. [2021], learning the forward bridge typically involves repeated
sampling from the reverse bridge, which significantly escalates computational
costs. To address these issues, we aim to develop a method that directly
learns the forward bridge, thereby bypassing the need for reversal. We also
anticipate applying this approach to additional applications, such as
stochastic image alignment and other contexts involving continuous data, to
further validate its versatility and effectiveness.

## References

  * Shi and Choi [2011] Jian Qing Shi and Taeryon Choi.  _Gaussian process regression analysis for functional data_.  CRC press, 2011.
  * Rogers and Williams [2000] Leonard CG Rogers and David Williams.  _Diffusions, markov processes, and martingales: Volume 1, foundations_ , volume 1.  Cambridge university press, 2000.
  * Franzese et al. [2024] Giulio Franzese, Giulio Corallo, Simone Rossi, Markus Heinonen, Maurizio Filippone, and Pietro Michiardi.  Continuous-time functional diffusion processes.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Lim et al. [2023] Jae Hyun Lim, Nikola B Kovachki, Ricardo Baptista, Christopher Beckham, Kamyar Azizzadenesheli, Jean Kossaifi, Vikram Voleti, Jiaming Song, Karsten Kreis, Jan Kautz, et al.  Score-based diffusion models in function space.  _arXiv preprint arXiv:2302.07400_ , 2023.
  * Pidstrigach et al. [2023] Jakiw Pidstrigach, Youssef Marzouk, Sebastian Reich, and Sven Wang.  Infinite-dimensional diffusion models for function spaces.  _arXiv e-prints_ , pages arXiv–2302, 2023.
  * Hyvärinen and Dayan [2005] Aapo Hyvärinen and Peter Dayan.  Estimation of non-normalized statistical models by score matching.  _Journal of Machine Learning Research_ , 6(4), 2005.
  * Vincent [2011] Pascal Vincent.  A connection between score matching and denoising autoencoders.  _Neural computation_ , 23(7):1661–1674, 2011.
  * Song et al. [2020] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  _International Conference on Learning Representations_ , 2020.
  * Huang et al. [2021] Chin-Wei Huang, Jae Hyun Lim, and Aaron C Courville.  A variational perspective on diffusion-based generative models and score matching.  _Advances in Neural Information Processing Systems_ , 34:22863–22876, 2021.
  * Baker et al. [2024] Elizabeth Louise Baker, Gefan Yang, Michael L Severinsen, Christy Anna Hipsley, and Stefan Sommer.  Conditioning non-linear and infinite-dimensional diffusion processes.  _arXiv preprint arXiv:2402.01434_ , 2024.
  * Delyon and Hu [2006] Bernard Delyon and Ying Hu.  Simulation of conditioned diffusion and application to parameter estimation.  _Stochastic Processes and their Applications_ , 116(11):1660–1675, 2006.
  * Schauer et al. [2017] Moritz Schauer, Frank van der Meulen, and Harry van Zanten.  Guided proposals for simulating multi-dimensional diffusion bridges.  _Bernoulli_ , 23(4A), November 2017.  ISSN 1350-7265.  doi: 10.3150/16-bej833.  URL http://dx.doi.org/10.3150/16-BEJ833.
  * Heng et al. [2021] Jeremy Heng, Valentin De Bortoli, Arnaud Doucet, and James Thornton.  Simulating diffusion bridges with score matching.  _arXiv preprint arXiv:2111.07243_ , 2021.
  * Park et al. [2023] Yesom Park, Jaemoo Choi, Changyeon Yoon, Myungjoo Kang, et al.  Learning pde solution operator for continuous modeling of time-series.  _arXiv preprint arXiv:2302.00854_ , 2023.
  * Ashiqur Rahman et al. [2022] Md Ashiqur Rahman, Zachary E Ross, and Kamyar Azizzadenesheli.  U-no: U-shaped neural operators.  _arXiv e-prints_ , pages arXiv–2204, 2022.
  * Van der Meulen and Schauer [2020] Frank Van der Meulen and Moritz Schauer.  Automatic backward filtering forward guiding for markov processes and graphical models.  _arXiv preprint arXiv:2010.03509_ , 2020.
  * Corstanje et al. [2024] Marc Corstanje, Frank van der Meulen, Moritz Schauer, and Stefan Sommer.  Simulating conditioned diffusions on manifolds.  _arXiv preprint arXiv:2403.05409_ , 2024.
  * Sohl-Dickstein et al. [2015] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli.  Deep unsupervised learning using nonequilibrium thermodynamics.  In _International conference on machine learning_ , pages 2256–2265. PMLR, 2015.
  * Song and Ermon [2019] Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  _Advances in neural information processing systems_ , 32, 2019.
  * Song et al. [2021] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=PxTIG12RRHS.
  * Baldassari et al. [2024] Lorenzo Baldassari, Ali Siahkoohi, Josselin Garnier, Knut Solna, and Maarten V de Hoop.  Conditional score-based diffusion models for bayesian inference in infinite dimensions.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Hagemann et al. [2023] Paul Hagemann, Sophie Mildenberger, Lars Ruthotto, Gabriele Steidl, and Nicole Tianjiao Yang.  Multilevel diffusion: Infinite dimensional score-based diffusion models for image generation.  _arXiv preprint arXiv:2303.04772_ , 2023.
  * De Bortoli et al. [2021] Valentin De Bortoli, James Thornton, Jeremy Heng, and Arnaud Doucet.  Diffusion schrödinger bridge with applications to score-based generative modeling.  _Advances in Neural Information Processing Systems_ , 34:17695–17709, 2021.
  * Shi et al. [2024] Yuyang Shi, Valentin De Bortoli, Andrew Campbell, and Arnaud Doucet.  Diffusion schrödinger bridge matching.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Tang et al. [2024] Zhicong Tang, Tiankai Hang, Shuyang Gu, Dong Chen, and Baining Guo.  Simplified diffusion schrödinger bridge.  _arXiv preprint arXiv:2403.14623_ , 2024.
  * Thornton et al. [2022] James Thornton, Michael Hutchinson, Emile Mathieu, Valentin De Bortoli, Yee Whye Teh, and Arnaud Doucet.  Riemannian diffusion schrödinger bridge.  _arXiv preprint arXiv:2207.03024_ , 2022.
  * Haussmann and Pardoux [1986] Ulrich G Haussmann and Etienne Pardoux.  Time reversal of diffusions.  _The Annals of Probability_ , pages 1188–1205, 1986.
  * Millet et al. [1989] Annie Millet, David Nualart, and Marta Sanz.  Time reversal for infinite-dimensional diffusions.  _Probability theory and related fields_ , 82(3):315–347, 1989.
  * Liu et al. [2023a] Guan-Horng Liu, Arash Vahdat, De-An Huang, Evangelos A Theodorou, Weili Nie, and Anima Anandkumar.  I2sb: Image-to-image schrödinger bridge.  _arXiv preprint arXiv:2302.05872_ , 2023a.
  * Delbracio and Milanfar [2023] Mauricio Delbracio and Peyman Milanfar.  Inversion by direct iteration: An alternative to denoising diffusion for image restoration.  _arXiv preprint arXiv:2303.11435_ , 2023.
  * Li et al. [2020a] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar, et al.  Fourier neural operator for parametric partial differential equations.  In _International Conference on Learning Representations_ , 2020a.
  * Li et al. [2020b] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Neural operator: Graph kernel network for partial differential equations.  _arXiv preprint arXiv:2003.03485_ , 2020b.
  * Lu et al. [2021] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  _Nature machine intelligence_ , 3(3):218–229, 2021.
  * Kovachki et al. [2023] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Neural operator: Learning maps between function spaces with applications to pdes.  _Journal of Machine Learning Research_ , 24(89):1–97, 2023.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Bayram et al. [2018] Mustafa Bayram, Tugcem Partal, and Gulsen Orucova Buyukoz.  Numerical methods for simulation of stochastic differential equations.  _Advances in Difference Equations_ , 2018:1–10, 2018.
  * Phillips et al. [2022] Angus Phillips, Thomas Seror, Michael Hutchinson, Valentin De Bortoli, Arnaud Doucet, and Emile Mathieu.  Spectral diffusion processes.  _arXiv preprint arXiv:2209.14125_ , 2022.
  * Da Prato and Zabczyk [2014] Giuseppe Da Prato and Jerzy Zabczyk.  _Stochastic equations in infinite dimensions_.  Cambridge university press, 2014.
  * Loges [1984] Wilfried Loges.  Girsanov’s theorem in hilbert space and an application to the statistics of hilbert space-valued stochastic differential equations.  _Stochastic processes and their applications_ , 17(2):243–263, 1984.
  * Younes [2010] Laurent Younes.  _Shapes and diffeomorphisms_ , volume 171.  Springer, 2010.
  * Bauer et al. [2014] Martin Bauer, Martins Bruveris, and Peter W Michor.  Overview of the geometries of shape spaces and diffeomorphism groups.  _Journal of Mathematical Imaging and Vision_ , 50:60–97, 2014.
  * Kunita and Kunita [1990] Hiroshi Kunita and Hiroshi Kunita.  _Stochastic flows and stochastic differential equations_ , volume 24.  Cambridge university press, 1990.
  * GBIF.Org User [2024] GBIF.Org User.  Occurrence download, 2024.  URL https://www.gbif.org/occurrence/download/0075323-231120084113126.
  * Kirillov et al. [2023] Alexander Kirillov, Eric Mintun, Nikhila Ravi, Hanzi Mao, Chloe Rolland, Laura Gustafson, Tete Xiao, Spencer Whitehead, Alexander C Berg, Wan-Yen Lo, et al.  Segment anything.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 4015–4026, 2023.
  * Liu et al. [2023b] Shilong Liu, Zhaoyang Zeng, Tianhe Ren, Feng Li, Hao Zhang, Jie Yang, Chunyuan Li, Jianwei Yang, Hang Su, Jun Zhu, et al.  Grounding dino: Marrying dino with grounded pre-training for open-set object detection.  _arXiv preprint arXiv:2303.05499_ , 2023b.
  * Adams and Otárola-Castillo [2013] Dean C Adams and Erik Otárola-Castillo.  geomorph: an r package for the collection and analysis of geometric morphometric shape data.  _Methods in ecology and evolution_ , 4(4):393–399, 2013.

##  Appendix A Proofs

In this section, we provide the proofs for two theorems in the paper.

###  A.1 Proof of Theorem 3.1

###### Proof.

We follow the setup in Millet et al. [1989]. Let
$I(i)\subset\mathbb{Z},\forall i\in\mathbb{Z}$ be a finite set of indices,
$\\{e_{i}\\}_{i=1}^{\infty}$ be an orthonormal basis for $H\ni X(t)$, and
$\\{k_{j}\\}_{j=1}^{\infty}$ be an orthonormal basis for $U\ni W(t)$. Then we
split $X(t)$ into an infinite number of finite vectors; specifically, we
denote $X^{i}(t)=(X_{j}(t),j\in I(i))\in\mathbb{R}^{I}$ to be the vector
consisting of all the components of $X$ indexed by the entries in $I(i)$, with
$X_{j}(t)=\langle X(t),e_{j}\rangle_{H}$. Similarly, let
$\hat{X}^{i}(t)=(X_{j}(t),j\notin I(i))$ be the complement of $X^{i}(t)$.
Then, by assumption, the conditional law of $X^{i}(t)$ given its complement
$\hat{X}^{i}(t)=\xi^{i}$ has a density $p_{t}(x_{i}\mid\xi^{i})$.

Consider the conditional bridge process (2), written in terms of the bases
$\\{e_{i}\\}$ and $\\{k_{j}\\}$:

$dX_{i}^{c}(t)=\left[f(t,X^{c}(t))+a(t,X^{c}(t))\nabla\log
h(t,X^{c}(t))\right]_{i}dt+\sum^{\infty}_{j=1}g_{ij}(t,X^{c}(t))dW_{j}(t).$
(16)

We apply the time reversal theorem [Millet et al., 1989], to obtain the
reversed bridge $Y^{c}(t)$:

$\displaystyle dY^{c}_{i}(t)$
$\displaystyle=-[f(T-t,Y^{c}(t))]_{i}dt-[a(T-t,Y^{c}(t))\nabla\log
h(T-t,Y^{c}(t))]_{i}dt$ (17a)
$\displaystyle+\left[\frac{1}{p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))}\sum_{j\in
I(i)}\nabla_{j}[a_{ij}(T-t,Y^{c}(t))p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))]\right]_{i}dt$
(17b) $\displaystyle+\sum^{\infty}_{j=1}g_{ij}(T-t,Y^{c}(t))dB_{j}(t).$ (17c)

In order to write this in the form appearing in the theorem, we focus on the
second term. By properties of the Fréchet derivative, it holds that

$\displaystyle\frac{1}{p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))}\sum_{j\in
I(i)}\nabla_{j}[a_{ij}(T-t,Y(t))p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))]$ (18a)
$\displaystyle=\sum_{j\in I(i)}\nabla_{j}[a_{ij}(T-t,Y(t))]$ (18b)
$\displaystyle+\frac{1}{p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))}\sum_{j\in
I(i)}a_{ij}(T-t,Y(t))\nabla_{j}[p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))]$ (18c)
$\displaystyle=\sum_{j\in I(i)}\nabla_{j}[a_{ij}(T-t,Y(t))]+\sum_{j\in
I(i)}a_{ij}(T-t,Y(t))\nabla_{j}\log[p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))].$
(18d)

We now compare this to the second term by expanding that. Given
$h(T-t,Y(t))=\mathbb{P}(Y(T)\mid Y(t))$, we expand it using conditional
densities:

$\displaystyle\mathbb{P}(Y^{c}(T)\mid
Y^{c}(t))=\sum_{i}^{\infty}p_{T|t}(Y^{c,i}(T)\mid
Y^{c,i}(t),\hat{Y}^{c,i}(t),\hat{Y}^{c,i}(T))$ (19)

Therefore,

$\displaystyle\begin{split}&[a(T-t,Y^{c}(t))\nabla\log h(T-t,Y^{c}(t))]_{i}\\\
=&\sum_{j\in I(i)}a_{ij}(T-t,Y^{c}(t))\nabla_{j}\log p_{T|t}(Y^{c,i}(T)\mid
Y^{c,i}(t),\hat{Y}^{c,i}(t),\hat{Y}^{c,i}(T))\end{split}$ (20)

Note that $p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))$ is a marginal density over
$0,T$, thereby can be written as:

$\displaystyle p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))$ (21a)
$\displaystyle=p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t),Y^{c,i}(0),\hat{Y}^{c,i}(0),Y^{c,i}(T),\hat{Y}^{c,i}(T))$
(21b)
$\displaystyle=\frac{p_{t|0}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t),Y^{c,i}(0)\mid\hat{Y}^{c,i}(0))p_{T|t}(Y^{c,i}(T)\mid\hat{Y}^{c,i}(T),Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))}{p_{T|0}(Y^{c,i}(T)\mid\hat{Y}^{c,i}(T),Y^{c,i}(0)\mid\hat{Y}^{c,i}(0))}$
(21c)

Taking the logarithm and the gradient with respect to $Y^{c,i}(t)$, we obtain:

$\displaystyle\nabla_{Y^{c,i}(t)}\log p_{t}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t))$
(22a) $\displaystyle=\nabla_{Y^{c,i}(t)}\log
p_{t|0}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t),Y^{c,i}(0),\hat{Y}^{c,i}(0))$ (22b)
$\displaystyle+\nabla_{Y^{c,i}(t)}\log
p_{T|t}(Y^{c,i}(T)\mid\hat{Y}^{c,i}(T),Y^{c,i}(t),\hat{Y}^{c,i}(t))$ (22c)

Therefore, the reversed bridge can be written with respect to the bases:

$\displaystyle dY^{c}_{i}(t)$ (23a)
$\displaystyle=\Bigg{[}-f(T-t,Y^{c}(t))+\sum_{j\in
I(i)}\nabla_{Y^{c,i}(t)}a_{ij}(T-t,Y^{c}(t))$ (23b) $\displaystyle-\sum_{j\in
I(i)}a_{ij}(T-t,Y^{c}(t))\nabla_{Y^{c,i}(t)}\log
p_{t|0}(Y^{c,i}(t)\mid\hat{Y}^{c,i}(t),Y^{c,i}(0),\hat{Y}^{c,i}(t))\Bigg{]}_{i}dt$
(23c) $\displaystyle+\sum^{\infty}_{j=1}g_{ij}(T-t,Y^{c}(t))dB_{j}(t)$ (23d)

By summing over indices $i$, we obtain the infinite-dimensional reversed
bridge as shown in Theorem 3.1. ∎

###  A.2 Proof of Theorem 4.1

###### Proof.

We first compute the Radon-Nikodym derivative of $\mathbb{P}^{x_{0}}$ with
respect to $\mathbb{Q}^{x_{0}}_{\theta}$ by Girsanov’s theorem (See Da Prato
and Zabczyk [2014], Theorem 10.14, Lemma 10.15, also see Loges [1984] Theorem
2), which states that there exists a $U$-valued stochastic process
$\psi(t,Y(t))$ that satisfies

$g(t,Y(t))C^{1/2}\psi(t,Y(t))=\tilde{f}(t,Y(t))-f_{\theta}(t,Y(t)),$ (24)

such that,

$\frac{\mathop{\mathrm{d}\mathbb{P}^{x_{0}}}}{\mathop{\mathrm{d}\mathbb{Q}^{x_{0}}_{\theta}}}(Y(t))=\exp\Big{\\{}\int^{T}_{0}\|\psi(t,Y(t))\|_{U}\mathop{\mathrm{d}B(t)}-\frac{1}{2}\int^{T}_{0}\|\psi(t,Y(t))\|_{U}^{2}\mathop{\mathrm{d}t}\Big{\\}}$

Thus, the KL divergence is:

$\displaystyle D_{KL}(\mathbb{P}^{x_{0}}\mid\mathbb{Q}^{x_{0}}_{\theta})$
(25a)
$\displaystyle=\mathbb{E}_{Y(t)}\left\\{\int^{T}_{0}\|\psi(t,Y(t)\|_{U}\mathop{\mathrm{d}B(t)}\right\\}+\mathbb{E}_{Y(t)}\left\\{-\frac{1}{2}\int^{T}_{0}\|\psi(t,Y(t))\|_{U}^{2}\mathop{\mathrm{d}t}\right\\}$
(25b)
$\displaystyle=\mathbb{E}_{X(t)}\left\\{\frac{1}{2}\int^{T}_{0}\|\psi(t,X(t))\|_{U}^{2}\mathop{\mathrm{d}t}\right\\}$
(25c)
$\displaystyle=\mathbb{E}_{X(t)}\left\\{\frac{1}{2}\int^{T}_{0}\left\|C^{-1/2}g^{-1}(t,X(t))\left\\{\tilde{f}(t,X(t))-f_{\theta}(t,X(t))\right\\}\right\|_{U}^{2}\mathop{\mathrm{d}t}\right\\}$
(25d)
$\displaystyle=\mathbb{E}_{X(t)}\left\\{\frac{1}{2}\int^{T}_{0}\left\|C^{-1/2}g^{-1}(t,X(t))\left\\{a(t,X(t))\nabla\log\tilde{h}(t,X(t))-\mathcal{G}_{\theta}(t,X(t))\right\\}\right\|_{U}^{2}\mathop{\mathrm{d}t}\right\\}$
(25e)

In the second equation, the first term vanishes because of the martingale
property of the Itô integral, while the second term flips the sign because of
the change of process from $Y(t)$ to $X(t)$ with the reversed time variable.
Denote $C^{-1/2}g^{-1}(t,x)$ by $A_{t}(x)$ for short, expand the squared norm
in (25):

$\displaystyle D_{KL}(\mathbb{P}^{x_{0}}\mid\mathbb{Q}^{x_{0}}_{\theta})$
$\displaystyle=\mathbb{E}_{X(t)}\left\\{\frac{1}{2}\int^{T}_{0}\left\|A_{t}(X(t))\left\\{a(t,X(t))\nabla\log\tilde{h}(t,X(t))-\mathcal{G}_{\theta}(t,X(t))\right\\}\right\|_{U}^{2}\mathop{\mathrm{d}t}\right\\}$
(26a)
$\displaystyle=\frac{1}{2}\int^{T}_{0}\int_{H}\left\|A_{t}(x)\left\\{a(t,x)\nabla\log\tilde{h}(t,x)-\mathcal{G}_{\theta}(t,x)\right\\}\right\|^{2}_{U}\mathop{\mathrm{d}\mu_{t}}(\mathop{\mathrm{d}x})\mathop{\mathrm{d}t}$
(26b)
$\displaystyle=\underbrace{\frac{1}{2}\int^{T}_{0}\int_{H}\left\|A_{t}(x)\mathcal{G}_{\theta}(t,x)\right\|^{2}_{U}\mathop{\mathrm{d}\mu_{t}}(\mathop{\mathrm{d}x})\mathop{\mathrm{d}t}}_{C_{1}}$
(26c)
$\displaystyle+\underbrace{\frac{1}{2}\int^{T}_{0}\int_{H}\left\|A_{t}(x)a(t,x)\nabla\log\tilde{h}(t,x)\right\|^{2}_{U}\mathop{\mathrm{d}\mu_{t}}(\mathop{\mathrm{d}x})\mathop{\mathrm{d}t}}_{C_{2}}$
(26d) $\displaystyle-\underbrace{\int^{T}_{0}\int_{H}\langle
A_{t}(x)\mathcal{G}_{\theta}(t,x),A_{t}(x)a(t,x)\nabla\log\tilde{h}(t,x)\rangle_{U}\mathop{\mathrm{d}\mu_{t}}(\mathop{\mathrm{d}x})\mathop{\mathrm{d}t}}_{C_{3}}$
(26e)

where $\mathop{\mathrm{d}\mu_{t}}(dx)$ denotes
$\mathbb{P}^{x_{0}}(X(t)\in\mathop{\mathrm{d}x}\mid
X(0))=P(0,x_{0},t,\mathop{\mathrm{d}x})$. Under the given partition
$(t_{i})_{i=1}^{I}$ of $[0,T]$, we shall split the time integral in $C_{3}$ by
finite parts:

$\displaystyle C_{3}$
$\displaystyle=\sum^{I}_{i=1}\int^{t_{i}}_{t_{i-1}}\int_{H}\langle
A_{t}(x)\mathcal{G}_{\theta}(t,x),A_{t}(x)a(t,x)\nabla\log\tilde{h}(t,x)\rangle_{U}\mathop{\mathrm{d}\mu_{t}}(\mathop{\mathrm{d}x})\mathop{\mathrm{d}t}$
(27a) $\displaystyle=\sum^{I}_{i=1}\int^{t_{i}}_{t_{i-1}}\int_{H}\langle
A_{t}(x)\mathcal{G}_{\theta}(t,x),A_{t}(x)a(t,x)\nabla\log
P(0,x_{0},t,\mathop{\mathrm{d}x})\rangle_{U}P(0,x_{0},t,\mathop{\mathrm{d}x})\mathop{\mathrm{d}t}$
(27b) $\displaystyle=\sum^{I}_{i=1}\int^{t_{i}}_{t_{i-1}}\int_{H}\langle
A_{t}(x)\mathcal{G}_{\theta}(t,x),A_{t}(x)a(t,x)\nabla
P(0,x_{0},t,\mathop{\mathrm{d}x})\rangle_{U}\mathop{\mathrm{d}t}$ (27c)

Using the Chapman-Kolmogorov equation [Da Prato and Zabczyk, 2014, Corollary
9.15], for any $0\leq t_{i-1}<t\leq T$, we show:

$\displaystyle\nabla P(0,x_{0},t,\mathop{\mathrm{d}x})$
$\displaystyle=\nabla_{X(t)}\int_{H}P(0,x_{0},t_{i-1},\mathop{\mathrm{d}x_{t_{i-1}}})P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})$
(28a)
$\displaystyle=\int_{H}\nabla_{X(t)}P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})P(0,x_{0},t_{i-1},\mathop{\mathrm{d}x_{t_{i-1}}})$
(28b) $\displaystyle=\int_{H}\nabla_{X(t)}\log
P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})P(0,x_{0},t_{i-1},\mathop{\mathrm{d}x_{t_{i-1}}})$
(28c)

Therefore, substitute (28) back into (27)

$\displaystyle C_{3}$
$\displaystyle=\sum^{I}_{i=1}\int^{t_{i}}_{t_{i-1}}\mathop{\mathrm{d}t}\int_{H}P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})P(0,x_{0},t_{i-1},\mathop{\mathrm{d}x_{t_{i-1}}})$
(29a) $\displaystyle\int_{H}\langle
A_{t}(x)\mathcal{G}_{\theta}(t,x),A_{t}(x)a(t,x)\nabla\log
P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})\rangle_{U}$ (29b)
$\displaystyle=\sum^{I}_{i=1}\int^{t_{i}}_{t_{i-1}}\mathop{\mathrm{d}t}P(0,x_{0},t,\mathop{\mathrm{d}x})\int_{H}\langle
A_{t}(x)\mathcal{G}_{\theta}(t,x),A_{t}(x)a(t,x)\nabla\log
P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})\rangle_{U}$ (29c)
$\displaystyle=\sum^{I}_{i=1}\int^{t_{i}}_{t_{i-1}}\int_{H}\langle
A_{t}(x)\mathcal{G}_{\theta}(t,x_{t}),A_{t}(x)a(t,x)\nabla\log
P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})\rangle_{U}\mathop{\mathrm{d}\mu_{t}}(\mathop{\mathrm{d}x})\mathop{\mathrm{d}t}$
(29d)

We then claim that:

$L(\theta):=\frac{1}{2}\sum^{I}_{i=1}\int^{t_{i}}_{t_{i-1}}\mathbb{E}_{X(t)}\left\\{\|A_{t}(x)\\{\mathcal{G}_{\theta}(t,x)-b(t,x,t_{i-1},x_{t_{i-1}})\\}\|_{U}^{2}\right\\}\mathop{\mathrm{d}t}=C_{1}+C_{4}-C_{3}$

where

$C_{4}=\frac{1}{2}\sum^{I}_{i=1}\int^{t_{i}}_{t_{i-1}}\int_{H}\left\|A_{t}(x)a(t,x)\nabla\log
P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})\right\|^{2}_{U}\mathop{\mathrm{d}\mu_{t}}(dx)\mathop{\mathrm{d}t}$

To show that, we expand the squared norm and get:

$\displaystyle L(\theta)$
$\displaystyle=\underbrace{\frac{1}{2}\sum_{i=1}^{I}\int^{t_{i}}_{t_{i-1}}\int_{H}\|A_{t}(x)\mathcal{G}_{\theta}(t,x)\|^{2}_{U}\mathop{\mathrm{d}\mu_{t}}(dx)\mathop{\mathrm{d}t}}_{C_{1}}$
(30a)
$\displaystyle+\underbrace{\frac{1}{2}\sum_{i=1}^{I}\int^{t_{i}}_{t_{i-1}}\int_{H}\|A_{t}(x)a(t,x_{t})\nabla\log
P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})\|^{2}_{U}\mathop{\mathrm{d}\mu_{t}}(\mathop{\mathrm{d}x})\mathop{\mathrm{d}t}}_{C_{4}}$
(30b)
$\displaystyle-\underbrace{\sum_{i=1}^{I}\int^{t_{i}}_{t_{i-1}}\int_{H}\langle
A_{t}(x)\mathcal{G}_{\theta}(t,x),A_{t}(x)a(t,x)\nabla\log
P(t_{i-1},x_{t_{i-1}},t,\mathop{\mathrm{d}x})\rangle_{H}\mathop{\mathrm{d}\mu_{t}}(\mathop{\mathrm{d}x})\mathop{\mathrm{d}t}}_{C_{3}}$
(30c)

Finally, by comparing with (26) and taking $C=C_{4}-C_{2}$, we finish the
proof. ∎

##  Appendix B Experiment Details

In this section, we provide the general algorithms for training and evaluating
the model, also more details on the three experiments. All the experiments are
done on the platform with one NVIDIA RTX 4090 GPU and Intel(R) Xeon(R) CPU
E5-2650 v4 @ 2.20GHz with 128GB memory and Ubuntu 22.04.4 LTS OS. The codes
are available on https://github.com/bookdiver/scoreoperator

###  B.1 Algorithms

Algorithm 1 Training

1:Bridging object measure $\mu^{d}_{0}$

2:Total time length $T$, time step $\delta t$

3:Batch size $B$

4:Initialized $\mathcal{G}_{\theta}$

5:Drift term $f^{d}$, diffusion term $g^{d}$

6:Form discrete time grid $\\{t_{0}=0,t_{1},\dots,t_{n-1},t_{n}=T\\},n=\lfloor
T/\delta t\rfloor$

7:while Not converged do

8: Sample $\\{X^{d}_{j}(0)\\}_{j=1}^{B}\sim\mu^{d}_{0}$, i.i.d.

9: for $\ \textbf{do}i\leftarrow 1,\dots,n$

10: Sample $\\{\varepsilon_{j}\\}_{j=1}^{B}\sim\mathcal{N}(0,I_{k})$, i.i.d

11: for $\ \textbf{do}j\leftarrow 1,\dots,B$

12: $X^{d}_{j}(t_{i})\leftarrow
X_{j}^{d}(t_{i-1})+f^{d}(t_{i-1},X_{j}^{d}(t_{i-1}))\delta
t+g^{d}(t_{i-1},X_{j}^{d}(t_{i-1}))\varepsilon_{j}$

13:
$b_{j}^{d}(t_{i},X_{j}^{d}(t_{i}),t_{i-1},X_{j}^{d}(t_{i-1}))\leftarrow-g^{d}(t_{i-1},X_{j}^{d}(t_{i-1}))\varepsilon_{j}/\delta
t$

14: end for

15: end for

16: $L^{d}(\theta)=\frac{1}{2B}\delta
t\sum^{B}_{j=1}\sum^{n}_{i=1}\|\mathcal{G}_{\theta}(t_{i},X_{j}^{d}(t_{i}))-b_{j}^{d}(t_{i},X_{j}^{d}(t_{i}),t_{i-1},X_{j}^{d}(t_{i-1})\|^{2}\\}$

17: Perform gradient descent step on $L^{d}(\theta)$

18:end while

Algorithm 2 Inference

1:Bridging object measures $\mu^{d}_{T}$

2:Total time length $T$, time step $\delta t$

3:Pretrained $\mathcal{G}_{\theta}$

4:Drift term $f^{d}$, diffusion term $g^{d}$

5:Form discrete time grid $\\{t_{0}=0,t_{1},\dots,t_{n-1},t_{n}=T\\},n=\lfloor
T/\delta t\rfloor$

6:Sample $X^{(c)d}(T)\sim\mu^{d}_{T}$, i.i.d.

7:for $\ \textbf{do}i\leftarrow n,\dots,1$

8: Sample $\varepsilon\sim\mathcal{N}(0,I_{k})$, i.i.d

9: $X^{(c)d}(t_{i-1})\leftarrow
X^{(c)d}(t_{i})+\\{f^{d}(t_{i},X^{(c)d}(t_{i}))+\mathcal{G}_{\theta}(t_{i},X^{(c)d}(t_{i}))+\nabla\cdot(g^{d}(g^{d})^{T})(t_{i},X^{(c)d}(t_{i}))\\}\delta
t+g^{d}(t_{i},X^{(c)d}(t_{i}))\varepsilon$

10:end for

###  B.2 Quadratic functions

We model the stochastic dynamics in the function space as the Brownian motion:

$\mathop{\mathrm{d}X(t)}=\sigma\mathop{\mathrm{d}W(t)}$ (31)

with $\sigma=0.1$, which has a closed-form time-reversed bridge:

$\mathop{\mathrm{d}Y(t)}=\sigma^{2}\frac{y_{T}-Y(t)}{T-t}\mathop{\mathrm{d}t}+\sigma\mathop{\mathrm{d}dB(t)}$
(32)

We designed our operator to have 6 Fourier layers. The details are shown in
Table 1. We train with 10,000 i.i.d. batched samples and with a batch size of
16. We chose the Adam optimizer with an initial learning rate of 0.001 and
cosine decay to 1e-5 as 80% of the training is finished.

Layer | Input | Output | Grid size | Fourier modes | Activation
---|---|---|---|---|---
Lifting | $u:\mathbb{R}\to\mathbb{R}$ | $v_{0}:\mathbb{R}\to\mathbb{R}^{16}$ | 8 | - | -
Down1 | $v_{0}:\mathbb{R}\to\mathbb{R}^{16}$ | $v_{1}:\mathbb{R}\to\mathbb{R}^{16}$ | 8 | 6 | GeLU
Down2 | $v_{1}:\mathbb{R}\to\mathbb{R}^{16}$ | $v_{2}:\mathbb{R}\to\mathbb{R}^{32}$ | 4 | 4 | GeLU
Down3 | $v_{2}:\mathbb{R}\to\mathbb{R}^{32}$ | $v_{3}:\mathbb{R}\to\mathbb{R}^{64}$ | 2 | 2 | GeLU
Up1 | $v_{3}:\mathbb{R}\to\mathbb{R}^{64}$ | $v_{4}:\mathbb{R}\to\mathbb{R}^{32}$ | 2 | 2 | GeLU
Up2 | $v_{4}:\mathbb{R}\to\mathbb{R}^{32}$ | $v_{5}:\mathbb{R}\to\mathbb{R}^{16}$ | 4 | 4 | GeLU
Up3 | $v_{5}:\mathbb{R}\to\mathbb{R}^{16}$ | $v_{6}:\mathbb{R}\to\mathbb{R}^{16}$ | 8 | 6 | GeLU
Projection | $v_{6}:\mathbb{R}\to\mathbb{R}^{16}$ | $v:\mathbb{R}\to\mathbb{R}$ | 8 | - | -
Table 1: Neural operator structure for quadratic function experiments Figure
7: More samples of estimated Brownian bridges between two quadratic functions

###  B.3 Circles

The SDE of circles follows the same as the one in the quadratic function
experiment. The bigger circle has a radius of 1.5, and the smaller circle has
a radius of 1.0. For the neural operator design, it is shown in Table 2. We
train with 10,000 i.i.d. batched samples and with a batch size of 16. We chose
the Adam optimizer with an initial learning rate of 0.001 and cosine decay to
1e-5 as 80% of the training is finished.

Layer | Input | Output | Grid size | Fourier modes | Activation
---|---|---|---|---|---
Lifting | $u:\mathbb{R}\to\mathbb{R}^{2}$ | $v_{0}:\mathbb{R}\to\mathbb{R}^{16}$ | 16 | - | -
Down1 | $v_{0}:\mathbb{R}\to\mathbb{R}^{16}$ | $v_{1}:\mathbb{R}\to\mathbb{R}^{16}$ | 16 | 8 | GeLU
Down2 | $v_{1}:\mathbb{R}\to\mathbb{R}^{16}$ | $v_{2}:\mathbb{R}\to\mathbb{R}^{32}$ | 8 | 6 | GeLU
Down3 | $v_{2}:\mathbb{R}\to\mathbb{R}^{32}$ | $v_{3}:\mathbb{R}\to\mathbb{R}^{64}$ | 4 | 4 | GeLU
Up1 | $v_{3}:\mathbb{R}\to\mathbb{R}^{64}$ | $v_{4}:\mathbb{R}\to\mathbb{R}^{32}$ | 4 | 4 | GeLU
Up2 | $v_{4}:\mathbb{R}\to\mathbb{R}^{32}$ | $v_{5}:\mathbb{R}\to\mathbb{R}^{16}$ | 8 | 6 | GeLU
Up3 | $v_{5}:\mathbb{R}\to\mathbb{R}^{16}$ | $v_{6}:\mathbb{R}\to\mathbb{R}^{16}$ | 16 | 8 | GeLU
Projection | $v_{6}:\mathbb{R}\to\mathbb{R}^{16}$ | $v:\mathbb{R}\to\mathbb{R}^{2}$ | 16 | - | -
Table 2: Neural operator structure for circle experiments Figure 8: More
samples of estimated Brownian bridges between two circles, evaluated under
different levels of discretization

###  B.4 Butterflies

We shall first briefly review the stochastic shape analysis. A shape
$s:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is usually modeled as an embedding in
$\mathbb{R}^{d}$ [Younes, 2010]. Matching two shapes $s_{0},s_{T}$ along time
is achieved by finding a diffeomorphism
$f:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$, such that
$f(0,s_{0}(\xi))=s_{0}(\xi)$ and $f(T,s_{0}(\xi))=s_{T}(\xi)$ for any
$\xi\in\mathbb{R}^{n}$. There are normally infinitely many numbers of such
$f$s that suffice, while only the one with the lowest energy is considered,
corresponding to the “simplest” deterministic transformation between the
shapes. Finding such a diffeomorphism is achieved by optimizing a given energy
functional [Bauer et al., 2014]. However, we are more interested in the
stochastic transformation between the shapes, that is, $f$ follows a function-
valued SDE that conditions on hitting close enough to $f(T):s_{0}(\xi)\mapsto
s_{T}(\xi)$ with the starting point as $f_{0}=\operatorname{id}$. One way of
defining $f$ is by the displacement of the shape, that is,
$f(t,s_{0}(\xi))=X(t,\xi)+s_{0}(\xi),x\in\mathbb{R}^{d}$, where $X\in
L^{2}([0,T]\times\mathbb{R}^{d},\mathbb{R}^{d})$ is a random variable in
Hilbert space. It is easy to show that $X(0,\cdot)=0,X(T,\cdot)=s_{T}-s_{0}$.
Furthermore, we define the unconditional SDE of $X$ as:

$dX(t)=\int^{t}_{0}\tilde{Q}(X(t))\mathop{\mathrm{d}W(t)},\quad X(0)=0$ (33)

Where for each shape $s\in
H=L^{2}(\mathbb{R}^{d},\mathbb{R}^{d})$,$\tilde{Q}(s):H\to H$ is a Hilbert-
Schmidt operator, defined by some smooth kernel $k\in
L^{2}(\mathbb{R}^{d}\times\mathbb{R}^{d},\mathbb{R}^{d}$) as:

$\tilde{Q}(s)(f(t,\xi))=\int_{\mathbb{R}^{d}}k(s(\xi)+\xi,\eta)f(t,\eta)\mathop{\mathrm{d}\eta},\quad\eta\in\mathbb{R}^{d}$
(34)

The SDE 33 is first defined in [Kunita and Kunita, 1990]. This formulation
ensures that the described process is diffeomorphic, which is required for
shape evolution since the shape change must be smooth and invertible. In our
experiment, we choose the butterfly shapes live in $\mathbb{R}^{2}$, and
choose the kernel $k$ to be 2D-Gaussian, i.e.,
$k(x,y)=\sigma\exp(-\|x-y\|^{2}/\kappa)$. Specifically, we set the domain in
$\mathbb{R}^{2}$ to be finite as $-0.5,1.5]^{2}$ and discretize it with the
resolution of $50\times 50$, we then choose the kernel with $\sigma=0.04$ and
$\kappa=0.02$, which allows only the correlations within small areas. We then
design the operator as Table 3. We train with 20,000 i.i.d. batched samples
and with a batch size of 16. We chose the Adam optimizer with an initial
learning rate of 0.001 and cosine decay to 1e-5 as 80% of the training is
finished.

The butterfly images are obtained from gbif.org [GBIF.Org User, 2024], and are
segmented with the Python package Segment Anything [Kirillov et al., 2023] and
Gounding Dino [Liu et al., 2023b]. The assigned landmarks are interpolated to
be evenly distributed along the outlines, then aligned by the R package
Germorph v.4.06 [Adams and Otárola-Castillo, 2013] and normalized into $[0,1]$
intervals.

Layer | Input | Output | Grid size | Fourier modes | Activation
---|---|---|---|---|---
Lifting | $u:\mathbb{R}\to\mathbb{R}^{2}$ | $v_{0}:\mathbb{R}\to\mathbb{R}^{16}$ | 32 | - | -
Down1 | $v_{0}:\mathbb{R}\to\mathbb{R}^{16}$ | $v_{1}:\mathbb{R}\to\mathbb{R}^{16}$ | 32 | 16 | GeLU
Down2 | $v_{1}:\mathbb{R}\to\mathbb{R}^{16}$ | $v_{2}:\mathbb{R}\to\mathbb{R}^{32}$ | 16 | 8 | GeLU
Down3 | $v_{2}:\mathbb{R}\to\mathbb{R}^{32}$ | $v_{3}:\mathbb{R}\to\mathbb{R}^{64}$ | 8 | 6 | GeLU
Down4 | $v_{3}:\mathbb{R}\to\mathbb{R}^{64}$ | $v_{4}:\mathbb{R}\to\mathbb{R}^{64}$ | 8 | 6 | GeLU
Up1 | $v_{4}:\mathbb{R}\to\mathbb{R}^{64}$ | $v_{5}:\mathbb{R}\to\mathbb{R}^{64}$ | 8 | 6 | GeLU
Up2 | $v_{5}:\mathbb{R}\to\mathbb{R}^{64}$ | $v_{6}:\mathbb{R}\to\mathbb{R}^{32}$ | 8 | 6 | GeLU
Up3 | $v_{6}:\mathbb{R}\to\mathbb{R}^{32}$ | $v_{7}:\mathbb{R}\to\mathbb{R}^{16}$ | 16 | 8 | GeLU
Up4 | $v_{7}:\mathbb{R}\to\mathbb{R}^{16}$ | $v_{8}:\mathbb{R}\to\mathbb{R}^{16}$ | 32 | 16 | GeLU
Projection | $v_{8}:\mathbb{R}\to\mathbb{R}^{16}$ | $v:\mathbb{R}\to\mathbb{R}^{2}$ | 32 | - | -
Table 3: Neural operator structure for butterfly experiments Figure 9: The
additional bridge simulations between different species, the bridges are
constructed between Archon apollinus (blue) and Allancastria cerisyi/Battus
polydamas (red). Figure 10: The additional bridge simulations between
different species, the bridges are constructed between Archon apollinus (blue)
and Cressida cressida/Graphium agamemon (red).