# Diffusion Schrödinger Bridges for Bayesian Computation

Jeremy<EMAIL_ADDRESS>[ Valentin De
<EMAIL_ADDRESS>[ Arnaud
<EMAIL_ADDRESS>[ Jeremy Heng is Assistant Professor,
ESSEC Business School, Singapore 139408, Singaporepresep= ]e1.  Valentin De
Bortoli is Research Scientist, Center for Sciences of Data, ENS Ulm, Paris,
Francepresep= ]e2.  Arnaud Doucet is Professor, Department of Statistics,
University of Oxford, Oxford OX1 3LB, U.K.presep= ]e3.

###### Abstract

Denoising diffusion models are a novel class of generative models that have
recently become extremely popular in machine learning. In this paper, we
describe how such ideas can also be used to sample from posterior
distributions and, more generally, any target distribution whose density is
known up to a normalizing constant. The key idea is to consider a forward
“noising” diffusion initialized at the target distribution which “transports”
this latter to a normal distribution for long diffusion times. The time-
reversal of this process, the “denoising” diffusion, thus “transports” the
normal distribution to the target distribution and can be approximated so as
to sample from the target. To accelerate simulation, we show how one can
introduce and approximate a Schrödinger bridge between these two
distributions, i.e. a diffusion which transports the normal to the target in
finite time.

Optimal transport,

Schrödinger bridge,

Score matching,

Stochastic differential equation,

Time-reversal,

###### keywords:

, and

##  1 Introduction

The use of diffusion processes as a mathematical model is ubiquitous in many
scientific disciplines. In differential form, such a process
$(X_{t})_{t\in[0,T]}$ in $\mathbb{R}^{d}$ is defined by a stochastic
differential equation (SDE) (see e.g. Øksendal (2003) and Klebaner (2012) for
textbooks on the subject)

$\displaystyle\mathrm{d}X_{t}=f(t,X_{t})\mathrm{d}t+\sigma(t,X_{t})\mathrm{d}B_{t}.$
(1)

The above describes infinitesimal changes in the process as a sum of
deterministic changes driven by a drift function
$f:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ for an infinitesimal
time step, and random fluctuations given by a diffusion matrix
$\sigma:[0,T]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$ and
infinitesimal changes of a standard Brownian motion in $\mathbb{R}^{d}$. The
interface between statistics and diffusion models have involved two main
threads: statistical inference for these models and their use for Monte Carlo
simulation. In the first, one considers inference for the parameters of $f$
and $\sigma$ when observations of the process are collected at various time
points. This has generated a large and comprehensive literature with
frequentist and Bayesian estimators, under various types of observation
regimes, and at various frequencies (see textbooks by Iacus (2008) and
Kessler, Lindner and Sørensen (2012) and references therein). In the second,
following ideas from statistical mechanics, the use of certain diffusion
processes and their long time behaviour have been adopted to simulate from
complex and high-dimensional distributions in computational statistics. For
example, by simulating Langevin dynamics whose equilibrium distribution is the
target distribution of interest, one progressively transforms an easy-to-
sample reference distribution to the desired target distribution. Indeed the
time-discretization of such Langevin dynamics have lead to efficient gradient-
based Markov chain Monte Carlo (MCMC) algorithms for sampling problems in
Bayesian statistics (Grenander and Miller, 1994; Roberts and Tweedie, 1996).

This article is based on denoising diffusions, which provide an alternative
approach to transform a simple reference distribution to a target
distribution. These diffusions were recently introduced in the machine
learning literature to obtain generative models (Sohl-Dickstein et al., 2015;
Ho, Jain and Abbeel, 2020; Song et al., 2021a). In contrast to running a
single realization of Langevin dynamics to obtain a large number of correlated
samples for MCMC, a single run of a denoising diffusion is used to simulate
one independent sample from the generative model. Given access to a large
dataset, say of images or proteins, generative models produce new synthetic
data with approximately the same distribution. Many generative modeling
techniques have been proposed in the literature, but the introduction of
denoising diffusions has recently revolutionized this field as these models
are easier to train than alternatives, and they provide state-of-the-art
performance in many domains. The main idea of these techniques is to
progressively transform the empirical data distribution into a normal
distribution through a forward “noising” diffusion process. By simulating an
approximation of the time-reversal of this process, a “denoising” diffusion,
we can generate new samples that closely follow the data distribution.

After introducing notation in Section 2, we review in Section 3 how one can
extend these ideas to perform approximate posterior simulation in scenarios
where one can only sample from the prior and simulate synthetic observations
of the likelihood (Song et al., 2021a). A limitation of this approach is that
it requires running the forward noising diffusion long enough to “transport”
the posterior distribution into an approximately normal distribution. Hence
the corresponding approximate time-reversal also needs to be run for a long
time to sample from the posterior. In Section 4, we show how a Schrödinger
bridge formulation of this problem can be exploited to transport in finite
time the posterior into a normal and vice-versa (Shi et al., 2022), hence
accelerating posterior simulation.

In scenarios where one only has access to an unnormalized density of a target
distribution, which might not correspond to the posterior distribution of a
Bayesian model, we show in Section 5 that the ideas behind denoising diffusion
models remain applicable. In this setting, one has to optimize an alternative
criterion to approximate the time-reversal of the forward noising diffusion
(Vargas, Grathwohl and Doucet, 2023). Finally we present a Schrödinger bridge
extension of this algorithm to speed up simulation in Section 6. To aid the
reading of this article, the features of algorithms covered are summarized in
Table 1. We conclude with some discussions in Section 7.

##  2 Notation

We will use throughout this paper the following notation. We will write
$x\in\mathcal{X}$ to denote latent states or parameters of a statistical
model, and $y\in\mathcal{Y}$ to denote observations. For example in Euclidean
spaces, we have $\mathcal{X}=\mathbb{R}^{d}$ and $\mathcal{Y}=\mathbb{R}^{p}$
for $p,d\in\mathbb{N}$. The Euclidean norm of $x\in\mathbb{R}^{d}$ is denoted
by $\|x\|$. For $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$, we denote its
gradient with respect to $x$ as $\nabla_{x}f(x)$. We write
$\mathcal{N}(\mu,\Sigma)$ to denote the normal distribution with mean $\mu$
and covariance $\Sigma$, and $x\mapsto\mathcal{N}(x;\mu,\Sigma)$ for its
probability density function. Given two probability measures $\mu$ and $\nu$
defined on a measurable space $\mathcal{Z}$, the Kullback–Leibler (KL)
divergence from $\mu$ to $\nu$ is defined as
$\mathrm{KL}(\mu|\nu)=\int_{\mathcal{Z}}\log(\mathrm{d}\mu/\mathrm{d}\nu)(z)\mathrm{d}\mu(z)$
if $\mu\ll\nu$ and $+\infty$ otherwise. We denote the set of integers
$[K]=\\{1,\dots,K\\}$. A path measure on $\mathcal{Z}$ is a probability
measure on $C([0,T],\mathcal{Z})$, where $C([0,T],\mathcal{Z})$ is the space
of continuous functions from $[0,T]$ to $\mathcal{Z}$. If $\mathbb{P}$ is a
path measure then we denote its _time-reversal_ as $\mathbb{P}^{R}$, defined
such that if $(X_{t})_{t\in[0,T]}\sim\mathbb{P}$ then
$(X_{T-t})_{t\in[0,T]}\sim\mathbb{P}^{R}$. We will also denote the time
$t\in[0,T]$ marginal distribution of $\mathbb{P}$ by $\mathbb{P}_{t}$, and
write its probability density function as $p_{t}$. Given a probability measure
$\mu$ on $\mathcal{Z}_{1}$, and a Markov kernel $\mathrm{K}$ from
$\mathcal{Z}_{1}$ to $\mathcal{Z}_{2}$, we define the probability measure
$\mu\otimes\mathrm{K}$ on $\mathcal{Z}_{1}\times\mathcal{Z}_{2}$ as
$(\mu\otimes\mathrm{K})(\mathsf{A}_{1}\times\mathsf{A}_{2})=\int_{\mathsf{A}_{1}}\mathrm{d}\mu(z_{1})\mathrm{K}(z_{1},\mathsf{A}_{2})$
for any Borel set $\mathsf{A}_{1}\times\mathsf{A}_{2}$. Finally, for
notational ease, we will use $(X_{t})_{t\in[0,T]}$ to denote _forward_
diffusion processes driven by Brownian motion $(B_{t})_{t\in[0,T]}$, and
$(Z_{t})_{t\in[0,T]}$ to denote _backward_ diffusion processes driven by
Brownian motion $(W_{t})_{t\in[0,T]}$.

Method name | Require samples | Require unnormalized density | Iterative | Section
---|---|---|---|---
Denoising Diffusion for Posterior Sampling (DDPS) | ✓ | ✗ | No | Section 3
Diffusion Schrödinger Bridge for Posterior Sampling (DSB-PS) | ✓ | ✗ | Yes | Section 4
Denoising Diffusion for General Sampling (DDGS) | ✗ | ✓ | No | Section 5
Diffusion Schrödinger Bridge for General Sampling (DSB-GS) | ✗ | ✓ | Yes | Section 6
Table 1: Summary of algorithms in this article.

##  3 Denoising diffusions for posterior sampling

We recall that $x\in\mathcal{X}$ denote latent states or parameters of a
statistical model and $y\in\mathcal{Y}$ denote observations. Given a prior
distribution $\mu(x)$ and likelihood function $x\mapsto g(y|x)$, Bayesian
inference is based on the posterior distribution

$\displaystyle p(x|y)=\frac{p(x,y)}{p(y)},\quad p(x,y)=\mu(x)g(y|x),$ (2)

where $p(y)=\int_{\mathcal{X}}\mu(x)g(y|x)\mathrm{d}x$ is the marginal
likelihood. We begin by describing an amortized variational inference
procedure that allows us to sample approximately from $p(x|y)$ for all
possible observations $y$. This method known as _denoising diffusions_ was
introduced by Song et al. (2021a) for $\mathcal{X}=\mathbb{R}^{d}$. We will
restrict ourselves here to this setup but it has been extended to general
spaces and discrete spaces (Austin et al., 2021; Hoogeboom et al., 2021;
Campbell et al., 2022; Benton et al., 2022; De Bortoli et al., 2022; Huang et
al., 2022). In this Section and Section 4, we assume that we can obtain
samples from the joint distribution $p(x,y)$. Practically we sample first
$X\sim\mu(x)$ followed by $Y\sim g(y|X)$.

We first define a forward “noising” process $(X_{t})_{t\in[0,T]}$ according to
the Ornstein–Uhlenbeck (OU) process

$\displaystyle\mathrm{d}X_{t}=-\tfrac{1}{2}X_{t}\mathrm{d}t+\mathrm{d}B_{t},\quad
X_{0}\sim p(x|y),$ (3)

where $(B_{t})_{t\in[0,T]}$ is a standard $d$-dimensional Brownian motion.
Writing the transition density of (3) as $p_{t|0}(x_{t}|x_{0})$, the
distribution of $X_{t}$ is given by
$p_{t}(x_{t}|y)=\int_{\mathcal{X}}p_{t|0}(x_{t}|x_{0})p(x_{0}|y)\mathrm{d}x_{0}$.
Since Equation (3) can be seen as the SDE of a Langevin diffusion with the
standard normal distribution $\mathcal{N}(0,I)$ as its stationary
distribution, by ergodicity properties of the process, $p_{T}(x_{T}|y)$ for
large time $T>0$ approaches $\mathcal{N}(0,I)$ for all $y$. The corresponding
time-reversed process $(Z_{t})_{t\in[0,T]}=(X_{T-t})_{t\in[0,T]}$ can be shown
to satisfy (weakly) the SDE (Anderson, 1982; Haussmann and Pardoux, 1986)

$\displaystyle\mathrm{d}Z_{t}=\tfrac{1}{2}Z_{t}\mathrm{d}t+\nabla_{z_{t}}\log
p_{T-t}(Z_{t}|y)\mathrm{d}t+\mathrm{d}W_{t},$ (4)

with another standard $d$-dimensional Brownian motion $(W_{t})_{t\in[0,T]}$
and $Z_{0}\sim p_{T}(x_{T}|y)$. If one could simulate the backward “denoising”
process $(Z_{t})_{t\in[0,T]}$, then $Z_{T}$ would be a sample from the
posterior distribution $p(x|y)$. Practically, we cannot sample from (4), so
consider instead a diffusion approximating it of the form

$\mathrm{d}Z_{t}=\tfrac{1}{2}Z_{t}\mathrm{d}t+s^{\theta}_{T-t}(Z_{t},y)\mathrm{d}t+\mathrm{d}W_{t},$
(5)

with $Z_{0}\sim\mathcal{N}(0,I)$. For the diffusion (5) to approximate (4), we
need i) the diffusion time $T$ to be large enough for
$p_{T}(x_{T}|y)\approx\mathcal{N}(x_{T};0,I)$; ii) an approximation of the
Stein score function $s^{\theta}_{t}(x_{t},y)\approx\nabla_{x_{t}}\log
p_{t}(x_{t}|y)$ for all
$(t,x_{t},y)\in[0,T]\times\mathcal{X}\times\mathcal{Y}$. We show next how such
an approximation can be obtained.

Let $\mathbb{P}_{y}$ denote the path measure on $C([0,T],\mathcal{X})$ induced
by (3) with observations $y\in\mathcal{Y}$ and let $\mathbb{Q}^{\theta}_{y}$
be the law on $C([0,T],\mathcal{X})$ induced by the time-reversal of (5). We
consider a parametric function class
$\\{s_{t}^{\theta}(x_{t},y):\theta\in\Theta\\}$ such as neural networks. We
obtain $\theta$ approximating the score by minimizing the expected KL
divergence between $\mathbb{P}_{y}$ and $\mathbb{Q}^{\theta}_{y}$ over $y\sim
p(y)$, which satisfies (Song et al., 2021b, Theorem 1)

$\displaystyle\mathcal{L}(\theta)$ $\displaystyle=2\mathbb{E}_{y\sim
p(y)}[\mathrm{KL}(\mathbb{P}_{y}|\mathbb{Q}^{\theta}_{y})]$
$\displaystyle\textstyle\equiv\int_{0}^{T}\mathbb{E}_{\overline{p}}\left[\|s_{t}^{\theta}(X_{t},Y)-\nabla_{x_{t}}\log
p_{t|0}(X_{t}|X_{0})\|^{2}\right]\mathrm{d}t,$ (6)

where ‘$\equiv$’ means equality up to an additive constant and the second
expectation in (6) is w.r.t.
$\overline{p}(x_{0},x_{t},y)=p(x_{0},y)p_{t|0}(x_{t}|x_{0})$. This loss
corresponds to a _denoising score matching_ (DSM) loss (Vincent, 2011). Under
the transition density of the OU process in (3), the gradient appearing in (6)
is $\nabla_{x_{t}}\log
p_{t|0}(x_{t}|x_{0})=\\{x_{0}\exp(-t/2)-x_{t}\\}/\\{1-\exp(-t)\\}$. The loss
$\mathcal{L}(\theta)$ can be minimized using stochastic gradient based
algorithms such as Adam (Kingma and Ba, 2014) if one can simulate
$(X_{0},Y)\sim p(x,y)$ from the model (2). Hence this approach is applicable
to problems where both the prior and the likelihood are intractable but one
can simulate parameters and synthetic data from them. From this perspective,
it is an alternative to Approximate Bayesian Computation (ABC) methods (Marin
et al., 2012; Sisson, Fan and Beaumont, 2018; Beaumont, 2019). Empirical
comparisons between denoising diffusions, MCMC, and ABC methods can be found
in (Benton et al., 2022; Sharrock et al., 2022; Shi et al., 2022).

Having found a minimizer $\theta$ of (6), the resulting denoising diffusion
posterior sampler (DDPS) described in Algorithm 1 involves simulating (5)
using a suitable numerical integrator (see e.g. Karras et al. (2022)) and
returning $Z_{T}$ as an approximate sample from the posterior distribution
$p(x|y)$ in (2). Various techniques that have been proposed to accelerate
denoising diffusion models are also applicable to DDPS (Dockhorn, Vahdat and
Kreis, 2022; Salimans and Ho, 2022). Finally, we note that one can tailor the
above procedure to specific observations $y$ in the case where the likelihood
function can be evaluated. This relies on the following alternative
decomposition to estimate the score

$\nabla_{x_{t}}\log
p_{t}(x_{t}|y)=\nabla_{x_{t}}\log\mu_{t}(x_{t})+\nabla_{x_{t}}\log
g_{t}(y|x_{t})$ (7)

In this setting, we have that
$\mu_{t}(x_{t})=\int_{\mathcal{X}}\mu(x_{0})p_{t|0}(x_{t}|x_{0})\mathrm{d}x_{0}$
is the diffused prior and
$g_{t}(y|x_{t})=\int_{\mathcal{X}}g(y|x_{0})p_{0|t}(x_{0}|x_{t})\mathrm{d}x_{0}$
is the modified likelihood. The first term on the r.h.s. of (7) can be
estimated by $s^{\theta}_{t}(x_{t})$ with $\theta$ obtained by minimizing
$\textup{KL}(\mathbb{P}|\mathbb{Q}^{\theta})$, where $\mathbb{P}$ is induced
by the noising diffusion (3) initialized using $X_{0}\sim\mu(x)$ and
$\mathbb{Q}^{\theta}$ is induced by the reversal of the form (5) with
$s^{\theta}_{t}(x_{t},y)$ replaced by $s^{\theta}_{t}(x_{t})$. Using
evaluations of $g(y|x)$, the term $g_{t}(y|x_{t})$ can be approximated using
regression, _conditional guidance_ techniques (Chung et al., 2023), or Monte
Carlo methods (Song et al., 2023).

The next section presents a principled framework based on Schrödinger bridges
to accelerate training and sampling in DDPS.

Algorithm 1 Denoising diffusion for posterior sampling

— Training procedure:

while not converged do

Sample $(X_{t}^{k})_{t\in[0,T]}$ using the SDE (3) where
$(X_{0}^{k},Y^{k})\sim p(x,y)$ for $k\in[K]$.

Approximate the loss (6) using $((X_{t}^{k})_{t\in[0,T]},Y^{k})_{k\in[K]}$.

Update $\theta$ in $s_{t}^{\theta}$ using Adam optimizer.

end while

— Sampling procedure:

Sample $(Z_{t})_{t\in[0,T]}$ using the SDE (5) for the observation $y$.

return: $Z_{T}$ approximately distributed according to $p(x|y)$.

##  4 Diffusion Schrödinger bridges for posterior sampling

For DDPS to perform well, the diffusion time $T$ has to be long enough so that
$p_{T}(x_{T}|y)\approx\mathcal{N}(x_{T};0,I)$. Shi et al. (2022) considered
the _Schrödinger bridge_ (SB) formulation described in (De Bortoli et al.,
2021; Vargas et al., 2021; Chen, Liu and Theodorou, 2022) to improve DDPS in
this setting. Let $\mathbb{D}$ denote the path measure on
$C([0,T],\mathcal{Y})$ induced by

$\displaystyle\mathrm{d}Y_{t}=0,\quad Y_{0}\sim p(y),$ (8)

i.e. $Y_{t}=Y_{0}$ for all $t\in[0,T]$. The SB is defined by the following
constrained KL minimization problem over probability measures on the path
space $C([0,T],\mathcal{X}\times\mathcal{Y})$:

$\displaystyle\textstyle\Pi^{\star}=\arg\min_{\Pi}\\{\mathrm{KL}(\Pi|\mathbb{P}):\Pi_{0}(x_{0},y_{0})=p(x_{0},y_{0}),$
$\displaystyle\textstyle\qquad\qquad\quad\quad\Pi_{T}(x_{T},y_{T})=\mathcal{N}(x_{T};0,I)p(y_{T})\\},$
(9)

where $\mathbb{P}=\mathbb{P}_{y_{0}}\otimes\mathbb{D}$ is the joint measure
associated to (3) and (8), and $\Pi_{t}$ denotes the time $t$ marginal
distribution under $\Pi$. As the minimizer has the form
$\Pi^{\star}=\mathbb{P}_{y_{0}}^{\star}\otimes\mathbb{D}$, by simulating a
backward process $(Z_{t})_{t\in[0,T]}$ with law $\mathbb{P}_{y}^{\star}$ and
$Z_{0}$ initialized from $\Pi_{T}^{\star}(x_{T}|y)=\mathcal{N}(x_{T};0,I)$,
the terminal state $Z_{T}$ returns a sample from the posterior distribution
$\Pi_{0}^{\star}(x_{0}|y)=p(x_{0}|y)$, $p(y)$-almost surely.

To solve for the SB, we apply the _iterative proportional fitting_ (IPF)
procedure (Fortet, 1940; Deming and Stephan, 1940; Kullback, 1968). This is
defined by the following recursions for $n\in\mathbb{N}_{0}$:

$\displaystyle\textstyle\Pi^{2n+1}=\arg\min_{\Pi}\\{\mathrm{KL}(\Pi|\Pi^{2n}):$
$\displaystyle\qquad\Pi_{T}(x_{T},y_{T})=\mathcal{N}(x_{T};0,I)p(y_{T})\\},$
$\displaystyle\textstyle\Pi^{2n+2}=\arg\min_{\Pi}\\{\mathrm{KL}(\Pi|\Pi^{2n+1}):$
$\displaystyle\qquad\qquad\Pi_{0}(x_{0},y_{0})=p(x_{0},y_{0})\\},$ (10)

initialized at $\Pi^{0}=\mathbb{P}$, which involves iterative KL-projections
to impose the marginal distribution $\Pi_{0}(x_{0},y_{0})=p(x_{0},y_{0})$ at
time $0$, and $\Pi_{T}(x_{T},y_{T})=\mathcal{N}(x_{T};0,I)p(y_{T})$ at time
$T$. Convergence of the iterates $(\Pi^{n})_{n\in\mathbb{N}}$ to the SB
$\Pi^{\star}$ has been studied under various sets of assumptions (Rüschendorf,
1995; Chen, Georgiou and Pavon, 2016; De Bortoli et al., 2021; Léger, 2021).
It can be shown that $\Pi^{2n+1}=\mathbb{P}_{y_{T}}^{2n+1}\otimes\mathbb{D}$
and $\Pi^{2n+2}=\mathbb{P}_{y_{0}}^{2n+2}\otimes\mathbb{D}$, where
$\mathbb{P}_{y_{T}}^{2n+1}$ and $\mathbb{P}_{y_{0}}^{2n+2}$ are path measures
on $C([0,T],\mathcal{X})$. The time-reversal path measure
$(\mathbb{P}_{y_{T}}^{2n+1})^{R}$ is induced by (11) while
$\mathbb{P}_{y_{0}}^{2n+2}$ is induced by (12)

$\displaystyle\mathrm{d}Z_{t}$
$\displaystyle=f_{T-t}^{2n+1}(Z_{t},y_{T})\mathrm{d}t+\mathrm{d}W_{t},\quad
Z_{0}\sim\mathcal{N}(0,I),$ (11) $\displaystyle\mathrm{d}X_{t}$
$\displaystyle=f_{t}^{2n+2}(X_{t},y_{0})\mathrm{d}t+\mathrm{d}B_{t},\quad
X_{0}\sim p(x|y_{0}),$ (12)

where $(B_{t})_{t\in[0,T]}$ and $(W_{t})_{t\in[0,T]}$ are $d$-dimensional
standard Brownian motions. Hence $(Z_{t})_{t\in[0,T]}$ given by (11)
corresponds to a backward process, while $(X_{t})_{t\in[0,T]}$ given by (12)
corresponds to a forward process. Sampling from the time-reversal of
$\Pi^{2n+1}=\mathbb{P}_{y_{T}}^{2n+1}\otimes\mathbb{D}$ requires additionally
sampling $Y_{T}\sim p(y)$ to ensure $(Z_{0},Y)\sim\Pi_{T}(z,y)$ while sampling
from $\Pi^{2n+2}=\mathbb{P}_{y_{0}}^{2n+2}\otimes\mathbb{D}$ requires sampling
$Y_{0}\sim p(y)$ to ensure $(X_{0},Y_{0})\sim\Pi_{0}(x,y)$. The above drift
functions satisfy the following recursions for $n\in\mathbb{N}_{0}$:

$\displaystyle f_{t}^{2n+1}(x_{t},y)$
$\displaystyle=-f_{t}^{2n}(x_{t},y)+\nabla_{x_{t}}\log\Pi_{t}^{2n}(x_{t}|y),$
(13) $\displaystyle f_{t}^{2n+2}(x_{t},y)$
$\displaystyle=-f_{t}^{2n+1}(x_{t},y)+\nabla_{x_{t}}\log\Pi_{t}^{2n+1}(x_{t}|y),$
(14)

initialized at $f_{t}^{0}(x_{t})=-x_{t}/2$. At iteration $n=0$, under the law
$\mathbb{P}_{y}^{1}$, the process (11) is exactly the denoising process (4)
with initialization $Z_{0}\sim\mathcal{N}(0,I)$. Hence we can proceed with DSM
to approximate $\nabla_{x_{t}}\log\Pi_{t}^{0}(x_{t}|y)$. This recovers DDPS of
Section 3; iterating further allows us to improve performance when $T$ is not
sufficiently large. The next iteration $n=1$ then defines a forward process
(12) under $\mathbb{P}_{y}^{2}$ by time-reversing the backward process
associated to $\mathbb{P}_{y}^{1}$ and initializing $(X_{0},Y_{0})\sim
p(x,y)$. Another application of DSM approximates the next score function
$\nabla_{x_{t}}\log\Pi_{t}^{1}(x_{t}|y)$ by simulating trajectories under
$\mathbb{P}_{y}^{1}$ using the previous score approximation. One can avoid
storing all score approximations by adopting the _mean-matching_ approach in
(De Bortoli et al., 2021; Shi et al., 2022). Iterating in this manner until
convergence yields a backward process (11) that we refer to as the Diffusion
Schrödinger Bridge Posterior Sampler (DSB-PS); see Algorithm 2 for an
algorithmic outline, and Shi et al. (2022, Section 4) for more implementation
details111Code for DDPS and DSB-PS is available here.. In practice, the drift
functions $f_{t}^{2n+1}$ and $f_{t}^{2n+2}$ are approximated by
$f_{t}^{\theta,2n+1}$ and $f_{t}^{\phi,2n+2}$.

Algorithm 2 Diffusion Schrödinger Bridge for posterior sampling

— Training procedure:

for $n=0,\dots,N$ do

– Approximation of $f_{t}^{2n+1}$ by $f_{t}^{\theta,2n+1}$

while not converged do

Sample $(X_{t}^{k})_{t\in[0,T]}$ using (12) where $(X_{0}^{k},Y_{0}^{k})\sim
p(x,y)$ with drift $f_{t}^{\phi,2n}$ for $k\in[K]$.

Set $Y_{t}^{k}=Y_{0}^{k}$ for $t\in(0,T]$ and $k\in[K]$.

Approximate the mean-matching loss (Shi et al., 2022, Section 3.1) using
$(X_{t}^{k},Y_{t}^{k})_{t\in[0,T],k\in[K]}$.

Update $\theta$ in $f_{t}^{\theta,2n+1}$ using Adam optimizer.

end while

– Approximation of $f_{t}^{2n+2}$ by $f_{t}^{\phi,2n+2}$

while not converged do

Sample $(Z_{t}^{k})_{t\in[0,T]}$ using (11) where
$(Z_{0}^{k},Y_{T}^{k})\sim\mathcal{N}(x;0,I)p(y)$ with drift
$f_{t}^{\theta,2n+1}$ for $k\in[K]$.

Set $Y_{t}^{k}=Y_{T}^{k}$ for $t\in[0,T)$ and $k\in[K]$.

Approximate the mean-matching loss (Shi et al., 2022, Section 3.1) using
$(Z_{t}^{k},Y_{t}^{k})_{t\in[0,T],k\in[K]}$.

Update $\phi$ in $f_{t}^{\phi,2n+2}$ using Adam optimizer.

end while

end for

— Sampling procedure:

Sample $(Z_{t})_{t\in[0,T]}$ using the SDE (11) for the observation $y$.

return: $Z_{T}$ approximately distributed according to $p(x|y)$.

##  5 Denoising diffusions for general sampling

As DDPS and DSB-PS are amortized variational inference procedures, we cannot
expect good performance for observations $y\in\mathcal{Y}$ that are very
unlikely under the model distribution $p(y)$, although this can be partially
mitigated by designing mechanisms to sample synthetic observations closer to
the available ones (Sharrock et al., 2022). These samplers are also not
applicable to general distributions that cannot be written as the posterior
distribution of a statistical model and for which we only have access to an
unnormalized density instead of samples from the joint model distribution. We
now present an algorithm proposed by Vargas, Grathwohl and Doucet (2023) that
can deal with these two limitations. Let $p(x)=\gamma(x)/\mathcal{Z}$ denote a
target distribution on $\mathcal{X}$, where
$\gamma:\mathcal{X}\rightarrow\mathbb{R}_{+}$ can be evaluated pointwise and
the normalizing constant $\mathcal{Z}=\int_{\mathcal{X}}\gamma(x)\mathrm{d}x$
is intractable.

Like before, we consider an OU process $(X_{t})_{t\in[0,T]}$ (3) with
initialization $X_{0}\sim p(x)$, and denote the induced path measure on
$C([0,T],\mathcal{X})$ as $\mathbb{P}$. The time-reversal
$(Z_{t})_{t\in[0,T]}=(X_{T-t})_{t\in[0,T]}$ satisfies (weakly) the SDE

$\displaystyle\mathrm{d}Z_{t}=\tfrac{1}{2}Z_{t}\mathrm{d}t+\nabla_{z_{t}}\log
p_{T-t}(Z_{t})\mathrm{d}t+\mathrm{d}W_{t},$ (15)

with $Z_{0}\sim p_{T}(x_{T})$ and
$p_{t}(x_{t})=\int_{\mathcal{X}}p_{t|0}(x_{t}|x_{0})p(x_{0})\mathrm{d}x_{0}$
denotes the marginal distribution of $X_{t}$. In contrast to Section 3, we
cannot approximate the score function $\nabla_{x_{t}}\log p_{t}(x_{t})$ by
minimizing a loss similar to (6) as it is infeasible to obtain or generate
samples from $p(x)$. Hence we shall seek an alternative representation of
$\mathbb{P}$. To do so, we introduce a reference path measure $\mathbb{M}$
induced by a stationary OU process, i.e. (3) with $X_{0}\sim\mathcal{N}(0,I)$.
This process is reversible, and its time-reversal satisfies

$\displaystyle\mathrm{d}Z_{t}=-\tfrac{1}{2}Z_{t}\mathrm{d}t+\mathrm{d}W_{t},\quad
Z_{0}\sim\mathcal{N}(0,I).$ (16)

For $s<t$, we denote the transition density of (16) as $m_{t|s}(z_{t}|z_{s})$.
Writing $\mathbb{P}=p\otimes\mathbb{M}_{|0}=\Phi\cdot\mathbb{M}$, where
$\mathbb{M}_{|0}$ denotes the law of $\mathbb{M}$ conditioned on $X_{0}$ and
$\Phi(x_{0})=p(x_{0})/\mathcal{N}(x_{0};0,I)$, we express $\mathbb{P}^{R}$ as
a _Doob’s $h$-transform_ (Rogers and Williams, 2000, p. 83) of (16) under
$\mathbb{M}$

$\displaystyle\mathrm{d}Z_{t}=-\tfrac{1}{2}Z_{t}\mathrm{d}t+\nabla_{z_{t}}\log
h_{T-t}(Z_{t})\mathrm{d}t+\mathrm{d}W_{t},$ (17)

with $Z_{0}\sim p_{T}(x_{T})$ and where we have

$\textstyle
h_{t}(x_{t})=\int_{\mathcal{X}}\Phi(x_{0})m_{T|T-t}(x_{0}|x_{t})\mathrm{d}x_{0},$

which can be characterized by a backward Kolmogorov equation and satisfies
$\nabla_{x_{t}}\log h_{t}(x_{t})=\nabla_{x_{t}}\log p_{t}(x_{t})+x_{t}$ from
(15). An implementation of (17) will involve setting $T$ sufficiently large
for $p_{T}\approx\mathcal{N}(0,I)$, and approximating $\nabla_{x_{t}}\log
h_{t}(x_{t})$ for all $(t,x_{t})\in[0,T]\times\mathcal{X}$ with a parametric
function class $\\{u_{t}^{\theta}(x_{t}):\theta\in\Theta\\}$. We write
$\mathbb{Q}^{\theta}$ for the law on $C([0,T],\mathcal{X})$ induced by the
time-reversal of the diffusion defined by $Z_{0}\sim\mathcal{N}(0,I)$ and

$\displaystyle\mathrm{d}Z_{t}=-\tfrac{1}{2}Z_{t}\mathrm{d}t+u_{T-t}^{\theta}(Z_{t})\mathrm{d}t+\mathrm{d}W_{t}.$
(18)

As we cannot sample from $\mathbb{P}$, we cannot easily obtain unbiased low-
variance estimates of the gradient of the forward KL loss
$\mathrm{KL}(\mathbb{P}|\mathbb{Q}^{\theta})$ w.r.t. $\theta$. So we instead
estimate $\theta$ by minimizing the reverse KL loss (Vargas, Grathwohl and
Doucet, 2023, Proposition 1)

$\displaystyle\textstyle\mathcal{L}(\theta)$
$\displaystyle=\mathrm{KL}(\mathbb{Q}^{\theta}|\mathbb{P})$
$\displaystyle=\textstyle\mathbb{E}_{\mathbb{Q}^{\theta}}[\frac{1}{2}\int_{0}^{T}\|u_{T-t}^{\theta}(Z_{t})\|^{2}\mathrm{d}t-\log\Phi(Z_{T})].$
(19)

This is a specific instance of a more general class of KL _optimal control_
problems with many connections to sampling (Kappen, Gómez and Opper, 2012;
Kappen and Ruiz, 2016; Heng et al., 2020; Zhang and Chen, 2022). We note that
intractability of the normalizing constant $\mathcal{Z}$ appearing in $\Phi$
is not an issue as it does not affect the minimizers of $\mathcal{L}(\theta)$.

With a minimizer $\theta$ of (5), we have a denoising diffusion general
sampler (DDGS, Algorithm 3) that gives approximate samples from $p(x)$ by
simulating (18) and returning $Z_{T}$, and an unbiased estimator of
$\mathcal{Z}$ using importance sampling with proposal law
$\mathbb{Q}^{\theta}$ and target law $\mathbb{P}$. An approximate sample from
$p(x)$ and an alternative estimator of $\mathcal{Z}$ can be obtained using a
_probability flow_ , i.e. an ordinary differential equation (ODE) that has the
same marginal distributions on $\mathcal{X}$ as $\mathbb{Q}^{\theta}$; see
Song et al. (2021a, Section 4.3), Vargas et al. (2023, Section 2.4), and
Appendix A. In practice, we have to consider time-discretizations of these
continuous-time processes and some care when choosing numerical integrators is
necessary (Vargas, Grathwohl and Doucet, 2023, Section 3). Empirical
comparisons between DDGS and Sequential Monte Carlo (SMC) methods can be found
in Vargas, Grathwohl and Doucet (2023).

Algorithm 3 Denoising diffusion for general sampling

— Training procedure:

while not converged do

Sample $(Z_{t}^{k})_{t\in[0,T]}$ using the SDE (18) for $k\in[K]$.

Approximate the loss (5) using $(Z_{t}^{k})_{t\in[0,T],k\in[K]}$.

Update $\theta$ in $u_{t}^{\theta}$ using Adam optimizer.

end while

— Sampling procedure:

Sample $(Z_{t})_{t\in[0,T]}$ using the SDE (18).

return: $Z_{T}$ approximately distributed according to $p(x)$.

##  6 Diffusion Schrödinger bridges for general sampling

Analogous to Section 4, we could also improve DDGS by adopting the following
SB formulation over probability measures on $C([0,T],\mathcal{X})$:

$\displaystyle\textstyle\Pi^{\star}=\arg\min_{\Pi}\\{\mathrm{KL}(\Pi|\mathbb{M}):$
$\displaystyle\qquad\Pi_{0}(x_{0})=p(x_{0}),\quad\Pi_{T}(x_{T})=\mathcal{N}(x_{T};0,I)\\}.$
(20)

We refer readers to Appendix B and De Bortoli et al. (2021, Section 3.1) for
connections to entropy-regularized _optimal transport_ problems and a
Monge–Kantorovich problem in the zero-noise limit (Mikami, 2004; Léonard,
2012, 2014). The IPF recursion solving (6) is

$\displaystyle\textstyle\Pi^{2n+1}=\arg\min_{\Pi}\\{\mathrm{KL}(\Pi|\Pi^{2n}):\Pi_{0}(x_{0})=p(x_{0})\\},$
$\displaystyle\textstyle\Pi^{2n+2}=\arg\min_{\Pi}\\{\mathrm{KL}(\Pi|\Pi^{2n+1}):\Pi_{T}(x_{T})=\pi_{T}(x_{T})\\},$

for $n\in\mathbb{N}_{0}$, with $\pi_{T}(x_{T})=\mathcal{N}(x_{T};0,I)$ and
$\Pi^{0}=\mathbb{M}$. The path measure $\Pi^{2n+1}$ is induced by the forward
process (21) while $(\Pi^{2n+2})^{R}$ is induced by the backward process (22)

$\displaystyle\mathrm{d}X_{t}$
$\displaystyle=f_{t}^{2n+1}(X_{t})\mathrm{d}t+\mathrm{d}B_{t},\quad X_{0}\sim
p(x),$ (21) $\displaystyle\mathrm{d}Z_{t}$
$\displaystyle=f_{T-t}^{2n+2}(Z_{t})\mathrm{d}t+\mathrm{d}W_{t},\quad
Z_{0}\sim\mathcal{N}(0,I),$ (22)

with drift functions

$\displaystyle f_{t}^{2n+1}(x_{t})$
$\displaystyle=-f_{t}^{2n}(x_{t})+\nabla_{x_{t}}\log\Pi_{t}^{2n}(x_{t}),$ (23)
$\displaystyle f_{t}^{2n+2}(x_{t})$
$\displaystyle=f_{t}^{2n}(x_{t})+\nabla_{x_{t}}\log h_{t}^{2n}(x_{t}),$ (24)

initialized at $f_{t}^{0}(x_{t})=-x_{t}/2$ and
$\Pi_{t}^{0}(x_{t})=\mathcal{N}(x_{t};0,I)$. Here

$\textstyle
h_{t}^{2n}(x_{t})=\int_{\mathcal{X}}\Phi^{2n}(x_{0})q_{T|T-t}^{2n}(x_{0}|x_{t})\mathrm{d}x_{0},$
(25)

where $\Phi^{2n}(x_{0})=p(x_{0})/\Pi_{0}^{2n}(x_{0})$ and
$q_{t|s}^{2n}(z_{t}|z_{s})$ for $s<t$ denotes the transition density of
$(Z_{t})_{t\in[0,T]}$ under $(\Pi^{2n})^{R}$.

Since the first two iterates are $\Pi^{1}=\mathbb{P}$ defined in Section 5 and
$\Pi^{2}=\mathcal{N}(0,I)\otimes\mathbb{P}_{|T}$, we recover DDGS of Section 5
by noticing that $\nabla_{x_{t}}\log h_{t}^{0}(x_{t})=\nabla_{x_{t}}\log
h_{t}(x_{t})$. As $p_{T}\to\mathcal{N}(0,I)$ when $T\to+\infty$, DDGS provides
an approximation of the SB when the time horizon $T$ is sufficiently large.
When the latter is not the case, further iterations can improve performance.
For $n\geq 1$, $\Pi^{2n+1}=p\otimes\Pi_{|0}^{2n}$ is the law of a forward
process (21) with a drift function (23) that involves
$\nabla_{x_{t}}\log\Pi_{t}^{2n}(x_{t})$. Simulating from our approximation of
$\Pi^{2n}$, DSM allows us to construct an approximation
$s_{t}^{\phi,2n}(x_{t})$ of this score function. Next by rewriting
$\Pi^{2n+2}=\mathcal{N}(0,I)\otimes\Pi_{|T}^{2n+1}=(h_{0}^{2n}/h_{T}^{2n})\cdot\Pi^{2n}$,
we see that its reversal is the law of a backward process (22) that is given
in terms of a Doob’s $h$-transform of the backward process associated to
$(\Pi^{2n})^{R}$; see Appendix C. To approximate $\nabla_{x_{t}}\log
h_{t}^{2n}(x_{t})$ with $u_{t}^{\theta,2n}(x_{t})$, we consider another
reverse KL minimization

$\displaystyle\textstyle\mathcal{L}^{2n+2}(\theta)$
$\displaystyle=\mathrm{KL}(\mathbb{Q}^{\theta}|\Pi^{2n+1})$
$\displaystyle\textstyle=\mathbb{E}_{\mathbb{Q}^{\theta}}[\frac{1}{2}\int_{0}^{T}\|u_{T-t}^{\theta,2n}(Z_{t})\|^{2}\mathrm{d}t-\log\Phi^{2n}(Z_{T})],$
(26)

where $(\mathbb{Q}^{\theta})^{R}$ is induced by the diffusion defined by
$Z_{0}\sim\mathcal{N}(0,I)$ and

$\displaystyle\mathrm{d}Z_{t}=f_{T-t}^{2n}(Z_{t})\mathrm{d}t+u_{T-t}^{\theta,2n}(Z_{t})\mathrm{d}t+\mathrm{d}W_{t}.$
(27)

In an implementation of (6) and (27), one has to replace $f_{t}^{2n}(x_{t})$
with the drift approximation from the previous KL minimizations, and
approximate $\Pi_{0}^{2n}(x_{0})$ in
$\Phi^{2n}(x_{0})=p(x_{0})/\Pi_{0}^{2n}(x_{0})$ by numerically integrating the
probability flow ODE using the score approximation $s_{t}^{\phi,2n}(x_{t})$.
One should iterate the above steps until convergence and avoid storing all
approximations in memory. The final backward process (22) yields the diffusion
Schrödinger bridge general sampler (DSB-GS), outlined in Algorithm 4. During
the training part of Algorithm 4, we can skip the approximation of
$\nabla\log\Pi_{t}^{2n}$ when $n=0$, since in that case
$\Pi_{t}^{0}=\mathcal{N}(0,I)$. Doing so, we see that the first iteration of
Algorithm 4 recovers Algorithm 3.

An alternative Schrödinger bridge formulation to sampling has been formulated
in (Föllmer, 1984; Dai Pra, 1991; Tzen and Raginsky, 2019). This corresponds
to (6) with the reference measure $\mathbb{M}$ defined using a pinned Brownian
motion running backwards in time, and the terminal distribution
$\Pi_{T}(x_{T})$ given by a Dirac measure at the origin. The appealing feature
of this simpler SB is that IPF converges in two iterations; various numerical
approximations have been developed for this case (Barr, Gispen and Lamacraft,
2020; Zhang, Sahai and Marzouk, 2021; Zhang and Chen, 2022; Vargas et al.,
2023). However, it was observed in Vargas, Grathwohl and Doucet (2023) that
the resulting scores one needs to estimate are very steep for times close to
$T$ due to the degenerate terminal distribution. As a result, this approach
can be numerically quite unstable.

Algorithm 4 Diffusion Schrödinger Bridge for general sampling

— Training procedure:

for $n=0,\dots,N$ do

– Approximation of $\nabla\log\Pi_{t}^{2n}$ by $s_{t}^{\phi,2n}$

while not converged do

Sample $(Z_{t}^{k})_{t\in[0,T]}$ using (22) with $f_{t}^{2n}$ for $k\in[K]$.

Approximate the DSM loss using $(Z_{t}^{k})_{t\in[0,T],k\in[K]}$

Update $\phi$ in $s_{t}^{\phi,2n}$ using Adam optimizer.

end while

– Approximation of $\nabla\log h_{t}^{2n}$ by $u_{t}^{\theta,2n}$

while not converged do

Sample $(Z_{t}^{k})_{t\in[0,T]}$ using the SDE (27) for $k\in[K]$.

Run probability flow ODE with $f_{t}^{2n}-\tfrac{1}{2}s_{t}^{\phi,2n}$ to
compute $\Pi_{0}^{2n}(Z_{T}^{k})$ for $k\in[K]$.

Approximate the KL loss (6) using $(Z_{t}^{k})_{t\in[0,T],k\in[K]}$.

Update $\theta$ in $u_{t}^{\theta,2n}$ using Adam optimizer.

end while

– Update drift $f_{t}^{2n+2}$

Set $f_{t}^{2n+2}=f_{t}^{2n}+u_{t}^{\theta,2n}$ as in (24).

end for

— Sampling procedure:

Sample $(Z_{t})_{t\in[0,T]}$ using the SDE (22).

return: $Z_{T}$ approximately distributed according to $p(x)$.

##  7 Discussion

We have presented a concise overview of how diffusion processes can be used to
sample approximately from posterior distributions and any general target
distributions. These methods have been successfully used to solve a wide class
of sampling problems, and are alternatives to standard MCMC, SMC, and ABC
techniques. However, compared to these well-established methods, while there
are convergence results available for such diffusion-based samplers (De
Bortoli, 2022; Chen et al., 2023), these results are based on assumptions
about the score function estimation error which remain difficult to verify.
Promising methodological alternatives to diffusion approaches have also been
recently proposed where one constructs a process between two distributions one
can sample from, but this process is built using an ordinary differential
equation whose drift function can be approximated by solving a simple
regression problem (Albergo and Vanden-Eijnden, 2023; Lipman et al., 2023;
Liu, Gong and Liu, 2023). The techniques developed in these works have been
further extended to compute the SB (Shi et al., 2023; Peluchetti, 2023) and
provide an alternative to DSB-PS.

##  Appendix A Probability Flow ODE

Consider the following SDE on $\mathcal{X}$

$\displaystyle\mathrm{d}X_{t}=f_{t}(X_{t})\mathrm{d}t+\mathrm{d}B_{t},\quad
X_{0}\sim p_{0}(x),$ (28)

where $(B_{t})_{t\in[0,T]}$ is a standard $d$-dimensional Brownian motion. The
marginal density $p_{t}$ of $X_{t}$ satisfies the Fokker–Planck–Kolmogorov
equation

$\displaystyle\partial_{t}p_{t}(x)$
$\displaystyle=-\mathrm{div}(f_{t}p_{t})(x)+\tfrac{1}{2}\Delta p_{t}(x)$ (29)
$\displaystyle=-\mathrm{div}(\bar{f}_{t}p_{t})(x)$ (30)

where $\mathrm{div}(u)(x):=\sum_{i=1}^{d}\tfrac{\partial u_{i}(x)}{\partial
x_{i}}$ for differentiable $u:\mathcal{X}\rightarrow\mathcal{X}$, $\Delta
a(x):=\sum_{i=1}^{d}\tfrac{\partial^{2}a(x)}{\partial x_{i}^{2}}$ for twice
differentiable $a:\mathcal{X}\rightarrow\mathbb{R}$, and

$\bar{f}_{t}(x)=f_{t}(x)-\tfrac{1}{2}\nabla\log p_{t}(x).$ (31)

Equation (30) shows that the ODE

$\displaystyle\mathrm{d}X_{t}=\bar{f}_{t}(X_{t})\mathrm{d}t,\quad X_{0}\sim
p_{0}(x),$ (32)

known as the _probability flow_ ODE, is such that $X_{t}\sim p_{t}(x)$ for all
$t\in[0,T]$, i.e. it admits the same marginal distributions as the SDE in
(28). By using the instantaneous change of variables formula (Chen et al.,
2018), we obtain

$\textstyle\log p_{0}(x_{0})=\log
p_{T}(x_{T})+\int_{0}^{T}\mathrm{div}(\bar{f}_{t})(x_{t})\mathrm{d}t,$ (33)

for $(x_{t})_{t\in[0,T]}$ obtained by solving (32) initialized at
$X_{0}=x_{0}$.

##  Appendix B Schrödinger Bridge as Entropy-Regularized Optimal Transport

Consider the following generic SB problem over probability measures on
$C([0,T],\mathcal{X})$:

$\displaystyle\textstyle\Pi^{\star}=\arg\min_{\Pi}\\{\mathrm{KL}(\Pi|\mathbb{S}):$
$\displaystyle\qquad\Pi_{0}(x_{0})=\nu_{0}(x_{0}),\quad\Pi_{T}(x_{T})=\nu_{T}(x_{T})\\},$
(34)

where $\mathbb{S}$ is a generic path measure on $C([0,T],\mathcal{X})$. We
note that the dynamic formulation (B) admits a static analogue. Let $s_{0,T}$
denote the marginal distribution of $(X_{0},X_{T})$ under $\mathbb{S}$, and
$\mathbb{S}_{|0,T}$ denote the law $\mathbb{S}$ conditioned on
$(X_{0},X_{T})$, i.e. we have $\mathbb{S}=s_{0,T}\otimes\mathbb{S}_{|0,T}$. By
decomposing the KL divergence in (B), we see that
$\Pi^{\star}=\pi_{0,T}^{\star}\otimes\mathbb{S}_{|0,T}$, where the static SB
$\pi_{0,T}^{\star}$ is defined by minimizing over probability measures on
$\mathcal{X}\times\mathcal{X}$ with the same marginal constraints

$\displaystyle\textstyle\pi_{0,T}^{\star}$
$\displaystyle=\textstyle\arg\min_{\pi_{0,T}}\\{\mathrm{KL}(\pi_{0,T}|s_{0,T}):$
$\displaystyle\pi_{0}(x_{0})=\nu_{0}(x_{0}),\quad\pi_{T}(x_{T})=\nu_{T}(x_{T})\\}.$
(35)

We can rewrite

$\mathrm{KL}(\pi_{0,T}|s_{0,T})=-\mathbb{E}_{\pi_{0,T}}[\log
s_{T|0}(X_{T}|X_{0})]-\text{H}(\pi_{0,T}),$

where $s_{T|0}(x_{T}|x_{0})$ denotes the transition density under $\mathbb{S}$
and $\text{H}(\pi_{0,T})=-\mathbb{E}_{\pi_{0,T}}[\log\pi_{0,T}(X_{0},X_{T})]$
is an entropy term. In particular, if $\mathbb{S}$ is the law of a standard
Brownian motion on $\mathcal{X}$, we have

$\displaystyle\textstyle\pi_{0,T}^{\star}$
$\displaystyle=\textstyle\arg\min_{\pi_{0,T}}\\{\mathbb{E}_{\pi_{0,T}}[||X_{T}-X_{0}||^{2}]-2T\text{H}(\pi_{0,T}):$
$\displaystyle\pi_{0}(x_{0})=\nu_{0}(x_{0}),\quad\pi_{T}(x_{T})=\nu_{T}(x_{T})\\}.$

The above can be seen as an entropy-regularized _optimal transport_ problem
that converges to a Monge–Kantorovich problem as $T\rightarrow 0$ (Mikami,
2004; Léonard, 2012, 2014).

##  Appendix C Iterative proportional fitting as Doob’s $h$-transforms

Recall from Section 6 that the IPF iterates satisfy

$\displaystyle\Pi^{2n+1}=p\otimes\Pi_{|0}^{2n},\quad\Pi^{2n+2}=\mathcal{N}(0,I)\otimes\Pi_{|T}^{2n+1}$
(36)

for $n\in\mathbb{N}_{0}$, which implies

$\displaystyle\Pi^{2n+2}=\frac{p(x_{0})\mathcal{N}(x_{T};0,I)}{\Pi_{0}^{2n}(x_{0})\Pi_{T}^{2n+1}(x_{T})}\Pi^{2n}.$
(37)

From the definition of $h_{t}^{2n}(x_{t})$ in (25), we note that

$\displaystyle
h_{0}^{2n}(x_{0})=\Phi^{2n}(x_{0})=p(x_{0})/\Pi_{0}^{2n}(x_{0}),$ (38)

and

$\displaystyle\textstyle
h_{T}^{2n}(x_{T})=\int_{\mathcal{X}}\Phi^{2n}(x_{0})q_{T|0}^{2n}(x_{0}|x_{T})\mathrm{d}x_{0}.$
(39)

Denoting the conditional distribution of $X_{T}$ given $X_{0}$ under
$\Pi^{2n}$ as $\pi^{2n}(x_{T}|x_{0})$, we have

$\displaystyle\textstyle\Pi_{T}^{2n+1}(x_{T})$
$\displaystyle=\textstyle\int_{\mathcal{X}}p(x_{0})\pi^{2n}(x_{T}|x_{0})\mathrm{d}x_{0}$
$\displaystyle\textstyle=\mathcal{N}(x_{T};0,I)\int_{\mathcal{X}}\Phi^{2n}(x_{0})q_{T|0}^{2n}(x_{0}|x_{T})\mathrm{d}x_{0},$
(40)

which gives

$\displaystyle\frac{\mathcal{N}(x_{T};0,I)}{\Pi_{T}^{2n+1}(x_{T})}=h_{T}^{2n}(x_{T})^{-1}.$
(41)

Combining (37), (38), and (41) gives

$\displaystyle\Pi^{2n+2}=\frac{h_{0}^{2n}(x_{0})}{h_{T}^{2n}(x_{T})}\Pi^{2n}.$
(42)

A.D. is partially supported by EPSRC grants EP/R034710/1 (CoSinES) and
EP/R018561/1 (Bayes4Health). J.H. was funded by CY Initiative of Excellence
(grant “Investissements d’Avenir” ANR-16-IDEX-0008).

## References

  * Albergo and Vanden-Eijnden (2023) [author] Albergo, Michael S.M. S. and Vanden-Eijnden, EricE. (2023). Building Normalizing Flows with Stochastic Interpolants. In The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.net.
  * Anderson (1982) [author] Anderson, Brian DOB. D. (1982). Reverse-time diffusion equation models. Stochastic Processes and Their Applications 12 313–326.
  * Austin et al. (2021) [author] Austin, JacobJ., Johnson, Daniel D.D. D., Ho, JonathanJ., Tarlow, DanielD. and van den Berg, RianneR. (2021). Structured Denoising Diffusion Models in Discrete State-Spaces. In Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual (Marc’AurelioM. Ranzato, AlinaA. Beygelzimer, Yann N.Y. N. Dauphin, PercyP. Liang and Jennifer WortmanJ. W. Vaughan, eds.) 17981–17993.
  * Barr, Gispen and Lamacraft (2020) [author] Barr, ArielA., Gispen, WillemW. and Lamacraft, AustenA. (2020). Quantum ground states from reinforcement learning. In Mathematical and Scientific Machine Learning.
  * Beaumont (2019) [author] Beaumont, Mark AM. A. (2019). Approximate Bayesian computation. Annual Review of Statistics and Its Application 6 379–403.
  * Benton et al. (2022) [author] Benton, JoeJ., Shi, YuyangY., De Bortoli, ValentinV., Deligiannidis, GeorgeG. and Doucet, ArnaudA. (2022). From denoising diffusions to denoising Markov models. arXiv preprint arXiv:2211.03595.
  * Campbell et al. (2022) [author] Campbell, AndrewA., Benton, JoeJ., De Bortoli, ValentinV., Rainforth, TomT., Deligiannidis, GeorgeG. and Doucet, ArnaudA. (2022). A Continuous Time Framework for Discrete Denoising Models. In Advances in Neural Information Processing Systems 35: Annual Conference on Neural Information Processing Systems 2022, NeurIPS 2022 (S.S. Koyejo, S.S. Mohamed, A.A. Agarwal, D.D. Belgrave, K.K. Cho and A.A. Oh, eds.).
  * Chen, Georgiou and Pavon (2016) [author] Chen, YongxinY., Georgiou, TryphonT. and Pavon, MicheleM. (2016). Entropic and displacement interpolation: a computational approach using the Hilbert metric. SIAM Journal on Applied Mathematics 76 2375–2396.
  * Chen, Liu and Theodorou (2022) [author] Chen, TianrongT., Liu, Guan-HorngG. and Theodorou, Evangelos A.E. A. (2022). Likelihood Training of Schrödinger Bridge using Forward-Backward SDEs Theory. In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net.
  * Chen et al. (2018) [author] Chen, Tian QiT. Q., Rubanova, YuliaY., Bettencourt, JesseJ. and Duvenaud, DavidD. (2018). Neural Ordinary Differential Equations. In Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada (SamyS. Bengio, Hanna M.H. M. Wallach, HugoH. Larochelle, KristenK. Grauman, NicolòN. Cesa-Bianchi and RomanR. Garnett, eds.) 6572–6583.
  * Chen et al. (2023) [author] Chen, SitanS., Chewi, SinhoS., Li, JerryJ., Li, YuanzhiY., Salim, AdilA. and Zhang, AnruA. (2023). Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. In The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.net.
  * Chung et al. (2023) [author] Chung, HyungjinH., Kim, JeongsolJ., McCann, Michael ThompsonM. T., Klasky, Marc LouisM. L. and Ye, Jong ChulJ. C. (2023). Diffusion Posterior Sampling for General Noisy Inverse Problems. In The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.net.
  * Dai Pra (1991) [author] Dai Pra, PaoloP. (1991). A stochastic control approach to reciprocal diffusion processes. Applied Mathematics and Optimization 23 313–329.
  * De Bortoli (2022) [author] De Bortoli, ValentinV. (2022). Convergence of denoising diffusion models under the manifold hypothesis. Transactions on Machine Learning Research 2022\.
  * De Bortoli et al. (2021) [author] De Bortoli, ValentinV., Thornton, JamesJ., Heng, JeremyJ. and Doucet, ArnaudA. (2021). Diffusion Schrödinger Bridge with Applications to Score-Based Generative Modeling. In Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual (Marc’AurelioM. Ranzato, AlinaA. Beygelzimer, Yann N.Y. N. Dauphin, PercyP. Liang and Jennifer WortmanJ. W. Vaughan, eds.) 17695–17709.
  * De Bortoli et al. (2022) [author] De Bortoli, ValentinV., Mathieu, EmileE., Hutchinson, MichaelM., Thornton, JamesJ., Teh, Yee WhyeY. W. and Doucet, ArnaudA. (2022). Riemannian score-based generative modeling. In Advances in Neural Information Processing Systems 35: Annual Conference on Neural Information Processing Systems 2022, NeurIPS 2022 (S.S. Koyejo, S.S. Mohamed, A.A. Agarwal, D.D. Belgrave, K.K. Cho and A.A. Oh, eds.).
  * Deming and Stephan (1940) [author] Deming, W EdwardsW. E. and Stephan, Frederick FF. F. (1940). On a least squares adjustment of a sampled frequency table when the expected marginal totals are known. The Annals of Mathematical Statistics 11 427–444.
  * Dockhorn, Vahdat and Kreis (2022) [author] Dockhorn, TimT., Vahdat, ArashA. and Kreis, KarstenK. (2022). Score-Based Generative Modeling with Critically-Damped Langevin Diffusion. In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net.
  * Föllmer (1984) [author] Föllmer, HansH. (1984). An entropy approach to the time reversal of diffusion processes. Lecture Notes in Control and Information Sciences 69 156–163.
  * Fortet (1940) [author] Fortet, RobertR. (1940). Résolution d’un système d’équations de M. Schrödinger. Journal de Mathématiques Pures et Appliquées 1 83–105.
  * Grenander and Miller (1994) [author] Grenander, UlfU. and Miller, Michael IM. I. (1994). Representations of knowledge in complex systems. Journal of the Royal Statistical Society: Series B (Methodological) 56 549–581.
  * Haussmann and Pardoux (1986) [author] Haussmann, Ulrich GU. G. and Pardoux, EtienneE. (1986). Time reversal of diffusions. The Annals of Probability 1188–1205.
  * Heng et al. (2020) [author] Heng, JeremyJ., Bishop, Adrian N.A. N., Deligiannidis, GeorgeG. and Doucet, ArnaudA. (2020). Controlled sequential Monte Carlo. The Annals of Statistics 48 2904 – 2929.
  * Ho, Jain and Abbeel (2020) [author] Ho, JonathanJ., Jain, AjayA. and Abbeel, PieterP. (2020). Denoising Diffusion Probabilistic Models. In Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual (HugoH. Larochelle, Marc’AurelioM. Ranzato, RaiaR. Hadsell, Maria-FlorinaM. Balcan and Hsuan-TienH. Lin, eds.).
  * Hoogeboom et al. (2021) [author] Hoogeboom, EmielE., Nielsen, DidrikD., Jaini, PriyankP., Forré, PatrickP. and Welling, MaxM. (2021). Argmax Flows and Multinomial Diffusion: Learning Categorical Distributions. In Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual (Marc’AurelioM. Ranzato, AlinaA. Beygelzimer, Yann N.Y. N. Dauphin, PercyP. Liang and Jennifer WortmanJ. W. Vaughan, eds.) 12454–12465.
  * Huang et al. (2022) [author] Huang, Chin-WeiC., Aghajohari, MiladM., Bose, JoeyJ., Panangaden, PrakashP. and Courville, Aaron C.A. C. (2022). Riemannian Diffusion Models. In Advances in Neural Information Processing Systems 35: Annual Conference on Neural Information Processing Systems 2022, NeurIPS 2022 (S.S. Koyejo, S.S. Mohamed, A.A. Agarwal, D.D. Belgrave, K.K. Cho and A.A. Oh, eds.).
  * Iacus (2008) [author] Iacus, Stefano MS. M. (2008). Simulation and Inference for Stochastic Differential Equations: with R Examples 486\. Springer.
  * Kappen, Gómez and Opper (2012) [author] Kappen, Hilbert JH. J., Gómez, VicençV. and Opper, ManfredM. (2012). Optimal control as a graphical model inference problem. Machine Learning 87 159–182.
  * Kappen and Ruiz (2016) [author] Kappen, Hilbert JohanH. J. and Ruiz, Hans ChristianH. C. (2016). Adaptive importance sampling for control and inference. Journal of Statistical Physics 162 1244–1266.
  * Karras et al. (2022) [author] Karras, TeroT., Aittala, MiikaM., Aila, TimoT. and Laine, SamuliS. (2022). Elucidating the Design Space of Diffusion-Based Generative Models. In Advances in Neural Information Processing Systems 35: Annual Conference on Neural Information Processing Systems 2022, NeurIPS 2022 (S.S. Koyejo, S.S. Mohamed, A.A. Agarwal, D.D. Belgrave, K.K. Cho and A.A. Oh, eds.).
  * Kessler, Lindner and Sørensen (2012) [author] Kessler, MathieuM., Lindner, AlexanderA. and Sørensen, MichaelM. (2012). Statistical Methods for Stochastic Differential Equations. Chapman & Hall.
  * Kingma and Ba (2014) [author] Kingma, Diederik PD. P. and Ba, JimmyJ. (2014). Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980.
  * Klebaner (2012) [author] Klebaner, Fima CF. C. (2012). Introduction to Stochastic Calculus with Applications. World Scientific Publishing Company.
  * Kullback (1968) [author] Kullback, SolomonS. (1968). Probability densities with given marginals. The Annals of Mathematical Statistics 39 1236–1243.
  * Léger (2021) [author] Léger, FlavienF. (2021). A gradient descent perspective on Sinkhorn. Applied Mathematics & Optimization 84 1843–1855.
  * Léonard (2012) [author] Léonard, ChristianC. (2012). From the Schrödinger problem to the Monge–Kantorovich problem. Journal of Functional Analysis 262 1879–1920.
  * Léonard (2014) [author] Léonard, ChristianC. (2014). A survey of the Schrödinger problem and some of its connections with optimal transport. Discrete and Continuous Dynamical Systems-Series A 34 1533–1574.
  * Lipman et al. (2023) [author] Lipman, YaronY., Chen, Ricky T. Q.R. T. Q., Ben-Hamu, HeliH., Nickel, MaximilianM. and Le, MatthewM. (2023). Flow Matching for Generative Modeling. In The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.net.
  * Liu, Gong and Liu (2023) [author] Liu, XingchaoX., Gong, ChengyueC. and Liu, QiangQ. (2023). Flow Straight and Fast: Learning to Generate and Transfer Data with Rectified Flow. In The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.net.
  * Marin et al. (2012) [author] Marin, Jean-MichelJ.-M., Pudlo, PierreP., Robert, Christian PC. P. and Ryder, Robin JR. J. (2012). Approximate Bayesian computational methods. Statistics and Computing 22 1167–1180.
  * Mikami (2004) [author] Mikami, ToshioT. (2004). Monge’s problem with a quadratic cost by the zero-noise limit of h-path processes. Probability Theory and Related Fields 129 245–260.
  * Øksendal (2003) [author] Øksendal, BerntB. (2003). Stochastic Differential Equations: An Introduction with Applications. Springer.
  * Peluchetti (2023) [author] Peluchetti, StefanoS. (2023). Diffusion Bridge Mixture Transports, Schrödinger Bridge Problems and Generative Modeling. arXiv preprint arXiv:2304.00917.
  * Roberts and Tweedie (1996) [author] Roberts, Gareth OG. O. and Tweedie, Richard LR. L. (1996). Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli 341–363.
  * Rogers and Williams (2000) [author] Rogers, L Chris GL. C. G. and Williams, DavidD. (2000). Diffusions, Markov Processes and Martingales: Volume 2, Itô calculus 2\. Cambridge university press.
  * Rüschendorf (1995) [author] Rüschendorf, LudgerL. (1995). Convergence of the iterative proportional fitting procedure. The Annals of Statistics 1160–1174.
  * Salimans and Ho (2022) [author] Salimans, TimT. and Ho, JonathanJ. (2022). Progressive Distillation for Fast Sampling of Diffusion Models. In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net.
  * Sharrock et al. (2022) [author] Sharrock, LouisL., Simons, JackJ., Liu, SongS. and Beaumont, MarkM. (2022). Sequential neural score estimation: Likelihood-free inference with conditional score based diffusion models. arXiv preprint arXiv:2210.04872.
  * Shi et al. (2022) [author] Shi, YuyangY., Bortoli, Valentin DeV. D., Deligiannidis, GeorgeG. and Doucet, ArnaudA. (2022). Conditional simulation using diffusion Schrödinger bridges. In Uncertainty in Artificial Intelligence, Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence, UAI 2022, 1-5 August 2022, Eindhoven, The Netherlands (JamesJ. Cussens and KunK. Zhang, eds.). Proceedings of Machine Learning Research 180 1792–1802. PMLR.
  * Shi et al. (2023) [author] Shi, YuyangY., De Bortoli, ValentinV., Campbell, AndrewA. and Doucet, ArnaudA. (2023). Diffusion Schrödinger Bridge Matching. arXiv preprint arXiv:2303.16852.
  * Sisson, Fan and Beaumont (2018) [author] Sisson, Scott AS. A., Fan, YananY. and Beaumont, MarkM. (2018). Handbook of Approximate Bayesian Computation. CRC Press.
  * Sohl-Dickstein et al. (2015) [author] Sohl-Dickstein, JaschaJ., Weiss, Eric A.E. A., Maheswaranathan, NiruN. and Ganguli, SuryaS. (2015). Deep Unsupervised Learning using Nonequilibrium Thermodynamics. In Proceedings of the 32nd International Conference on Machine Learning, ICML 2015, Lille, France, 6-11 July 2015 (Francis R.F. R. Bach and David M.D. M. Blei, eds.). JMLR Workshop and Conference Proceedings 37 2256–2265. JMLR.org.
  * Song et al. (2021a) [author] Song, YangY., Sohl-Dickstein, JaschaJ., Kingma, Diederik P.D. P., Kumar, AbhishekA., Ermon, StefanoS. and Poole, BenB. (2021a). Score-Based Generative Modeling through Stochastic Differential Equations. In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net.
  * Song et al. (2021b) [author] Song, YangY., Durkan, ConorC., Murray, IainI. and Ermon, StefanoS. (2021b). Maximum Likelihood Training of Score-Based Diffusion Models. In Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual (Marc’AurelioM. Ranzato, AlinaA. Beygelzimer, Yann N.Y. N. Dauphin, PercyP. Liang and Jennifer WortmanJ. W. Vaughan, eds.) 1415–1428.
  * Song et al. (2023) [author] Song, JiamingJ., Zhang, QinshengQ., Yin, HongxuH., Mardani, MortezaM., Liu, Ming-YuM.-Y., Kautz, JanJ., Chen, YongxinY. and Vahdat, ArashA. (2023). Loss-Guided Diffusion Models for Plug-and-Play Controllable Generation. In Proceedings of the 40th International Conference on Machine Learning (AndreasA. Krause, EmmaE. Brunskill, KyunghyunK. Cho, BarbaraB. Engelhardt, SivanS. Sabato and JonathanJ. Scarlett, eds.). Proceedings of Machine Learning Research 202 32483–32498. PMLR.
  * Tzen and Raginsky (2019) [author] Tzen, BelindaB. and Raginsky, MaximM. (2019). Theoretical guarantees for sampling and inference in generative models with latent diffusions. In Conference on Learning Theory, COLT 2019, 25-28 June 2019, Phoenix, AZ, USA (AlinaA. Beygelzimer and DanielD. Hsu, eds.). Proceedings of Machine Learning Research 99 3084–3114. PMLR.
  * Vargas, Grathwohl and Doucet (2023) [author] Vargas, FranciscoF., Grathwohl, Will SussmanW. S. and Doucet, ArnaudA. (2023). Denoising Diffusion Samplers. In The Eleventh International Conference on Learning Representations, ICLR 2023, Kigali, Rwanda, May 1-5, 2023. OpenReview.net.
  * Vargas et al. (2021) [author] Vargas, FranciscoF., Thodoroff, PierreP., Lamacraft, AustenA. and Lawrence, NeilN. (2021). Solving Schrödinger bridges via maximum likelihood. Entropy 23 1134\.
  * Vargas et al. (2023) [author] Vargas, FranciscoF., Ovsianas, AndriusA., Fernandes, DavidD., Girolami, MarkM., Lawrence, Neil DN. D. and Nüsken, NikolasN. (2023). Bayesian learning via neural Schrödinger–Föllmer flows. Statistics and Computing 33 3\.
  * Vincent (2011) [author] Vincent, PascalP. (2011). A connection between score matching and denoising autoencoders. Neural Computation 23 1661–1674.
  * Zhang and Chen (2022) [author] Zhang, QinshengQ. and Chen, YongxinY. (2022). Path Integral Sampler: A Stochastic Control Approach For Sampling. In The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net.
  * Zhang, Sahai and Marzouk (2021) [author] Zhang, BenjaminB., Sahai, TuhinT. and Marzouk, YoussefY. (2021). Sampling via controlled stochastic dynamical systems. In I (Still) Can’t Believe It’s Not Better! NeurIPS 2021 Workshop.