# Particle Denoising Diffusion Sampler

Angus Phillips  Hai-Dang Dau  Michael John Hutchinson  Valentin De Bortoli
George Deligiannidis  Arnaud Doucet

###### Abstract

Denoising diffusion models have become ubiquitous for generative modeling. The
core idea is to transport the data distribution to a Gaussian by using a
diffusion. Approximate samples from the data distribution are then obtained by
estimating the time-reversal of this diffusion using score matching ideas. We
follow here a similar strategy to sample from unnormalized probability
densities and compute their normalizing constants. However, the time-reversed
diffusion is here simulated by using an original iterative particle scheme
relying on a novel score matching loss. Contrary to standard denoising
diffusion models, the resulting Particle Denoising Diffusion Sampler (PDDS)
provides asymptotically consistent estimates under mild assumptions. We
demonstrate PDDS on multimodal and high dimensional sampling tasks.

##  1 Introduction

Consider a target probability density $\pi$ on $\mathbb{R}^{d}$ of the form

$\pi(x)=\frac{\gamma(x)}{\mathcal{Z}},\qquad\mathcal{Z}=\int_{\mathbb{R}^{d}}\gamma(x)\mathrm{d}x,\vspace{-0.2cm}$
(1)

where $\gamma:\mathbb{R}^{d}\rightarrow\mathbb{R}^{+}$ can be evaluated
pointwise but its normalizing constant $\mathcal{Z}$ is intractable. We
develop here a Monte Carlo scheme to sample approximately from $\pi$ and
estimate $\mathcal{Z}$.

We follow an approach inspired by denoising diffusion models (Ho et al., 2020;
Song et al., 2021; Song & Ermon, 2019) by considering a “noising” diffusion
progressively transporting the original target to a Gaussian. The time-
reversal of this diffusion, the “denoising” diffusion, allows us theoretically
to sample from the target starting from Gaussian noise. However, it is
impossible to simulate this process exactly as its drift depends on the
gradient of the logarithm of the intractable marginal densities of the noising
diffusion, i.e. the score. For generative modeling, where one has access to
samples from $\pi$, one can rely on neural networks and score matching
(Hyvärinen, 2005; Vincent, 2011). This strategy is not applicable in the Monte
Carlo sampling context as we cannot sample the “noising” diffusion since we do
not have access to any samples from $\pi$ to approximate the initial
distribution of the diffusion with.

The idea of using denoising diffusion models for Monte Carlo sampling has
already been explored by Berner et al. (2022); McDonald & Barron (2022);
Vargas et al. (2023); Huang et al. (2024); Richter et al. (2024). Berner et
al. (2022); Richter et al. (2024); Vargas et al. (2023) focus on the
minimization of a reverse Kullback–Leibler divergence or log-variance
criterion while McDonald & Barron (2022) rely on an importance sampling scheme
which scales poorly in high dimensions. Finally, Huang et al. (2024) relies on
a series of Markov chain Monte Carlo (MCMC) to estimate the score. This last
scheme does not provide estimates of normalizing constants.

We develop here an alternative approach inspired by denoising diffusion models
with guidance. Guided diffusions combine pre-trained diffusion models with a
guidance term derived from a likelihood to sample approximately from posterior
distributions; see e.g. Song et al. (2021); Chung et al. (2023); Song et al.
(2023); Corso et al. (2023). While they provide samples with appealing
perceptual properties, they rely on various approximations, in order of
importance: (1) approximation of the score and guidance terms, (2) time-
discretization of the diffusion and (3) approximate initialization of the
diffusion. Wu et al. (2023); Cardoso et al. (2024) have used particle methods
also known as Sequential Monte Carlo (SMC) (Doucet et al., 2001; Chopin &
Papaspiliopoulos, 2020) to obtain consistent estimates in the generative
modeling context.

Our contributions are as follows:
(1) we adapt guided diffusions to sampling problems,
(2) we provide theoretical results quantifying the error introduced by current
guided diffusions in a simple scenario,
(3) we develop an SMC scheme to provide consistent estimates in this setup and
establish limit theorems,
(4) we introduce an algorithm that reduces the variance of the SMC estimates
based on a novel score matching loss.
All proofs are postponed to the Appendix.

##  2 Denoising diffusions with guidance

###  2.1 Noising and denoising diffusions

Consider the following noising diffusion $(X_{t})_{t\in[0,T]}$,

$\differential X_{t}=-\beta_{t}X_{t}\differential
t+\sqrt{2\beta_{t}}\differential W_{t},\quad X_{0}\sim\pi,$ (2)

where $(W_{t})_{t\in[0,T]}$ is a $d$-dimensional Brownian motion and
$\beta_{t}>0$. The transition density of this diffusion is given by
$p(x_{t}|x_{0})=\mathcal{N}(x_{t};\sqrt{1-\lambda_{t}}x_{0},\lambda_{t}\mathrm{I})$
for $\textstyle\lambda_{t}=1-\exp[-2\int_{0}^{t}\beta_{s}\differential s]$. We
denote by $\pi_{t}$ the density of $X_{t}$ under (2). In practice, we consider
$\int_{0}^{T}\beta_{s}\differential s\gg 1$, and therefore
$\pi_{T}(x)\approx\mathcal{N}(x;0,\mathrm{I})$. The diffusion (2) thus
transforms $\pi_{0}=\pi$ into approximately $\mathcal{N}(0,\mathrm{I})$. If
instead we initialize (2) using $p_{0}(x)=\mathcal{N}(x;0,\mathrm{I})$, its
marginals $(p_{t})_{t\in[0,T]}$ satisfy $p_{t}=p_{0}$.

The time-reversal $(Y_{t})_{t\in[0,T]}=(X_{T-t})_{t\in[0,T]}$ of (2), the
denoising diffusion, satisfies $Y_{T}\sim\pi$ and

$\differential
Y_{t}=\left[\beta_{T-t}Y_{t}+2\beta_{T-t}\nabla\log\pi_{T-t}(Y_{t})\right]\differential
t+\sqrt{2\beta_{T-t}}\differential B_{t},$ (3)

where $(B_{t})_{t\in[0,T]}$ is another Brownian motion; see e.g. Haussmann &
Pardoux (1986); Cattiaux et al. (2023). The main idea of denoising diffusions
is to sample from $\pi$ by sampling (3) as $Y_{T}\sim\pi$ (Ho et al., 2020;
Song et al., 2021). However, we cannot simulate (3) exactly as, in order of
importance, (1) the score terms $(\nabla\log\pi_{t})_{t\in[0,T]}$ are
intractable, (2) it is necessary to time-discretize the diffusion and (3)
$\pi_{T}$ cannot be sampled. We can always use numerical integrators and
approximate $\pi_{T}$ with a unit Gaussian distribution to mitigate (2) and
(3). In generative modeling (1) is addressed by leveraging tools from the
score matching literature (Vincent, 2011; Hyvärinen, 2005) and using neural
network estimators. In our _sampling_ setting we do not have access to access
to samples from $\pi$ but only to its unnormalized density. Therefore,
alternative approximations must be developed.

###  2.2 Denoising diffusions with guidance

For generative modeling, the use of denoising diffusions with guidance terms
to sample approximately from posterior distributions has become prominent in
the inverse problem literature, see e.g. Song et al. (2021); Chung et al.
(2023); Song et al. (2023); Corso et al. (2023). We present here a simple
extension of this idea applicable to any target $\pi(x)$ by rewriting this
target as

$\pi(x)\propto
p_{0}(x)g_{0}(x),\quad\text{for~{}~{}}p_{0}(x)=\mathcal{N}(x;0,\mathrm{I}).$
(4)

It is always possible to perform such a decomposition, which provides an
alternative expression for the score.

######  Lemma 2.1.

The following identities hold

$\pi_{t}(x_{t})\propto
p_{0}(x_{t})g_{t}(x_{t}),\nabla\log\pi_{t}(x_{t})=-x_{t}+\nabla\log
g_{t}(x_{t}),$ (5)

where

$\textstyle g_{t}(x_{t})=\int g_{0}(x_{0})p(x_{0}|x_{t})\differential x_{0},$
(6)

and
$p(x_{0}|x_{t})=\mathcal{N}(x_{0};\sqrt{1-\lambda_{t}}x_{t},\lambda_{t}\mathrm{I})$
is the conditional density of $X_{0}$ given $X_{t}=x_{t}$ for the diffusion
(2) initialized using $X_{0}\sim p_{0}$.

From Lemma 2.1, it follows that (3) can be rewritten as

$\displaystyle\differential Y_{t}=$
$\displaystyle\left[-\beta_{T-t}Y_{t}+2\beta_{T-t}\nabla\log
g_{T-t}(Y_{t})\right]\differential t$ (7)
$\displaystyle+\sqrt{2\beta_{T-t}}\differential B_{t}.$

This is akin to having a diffusion model with tractable scores $\nabla\log
p_{t}(x_{t})=\nabla\log p_{0}(x_{t})=-x_{t}$ and with $g_{t}(x_{t})$ as a
guidance term.

###  2.3 Guidance approximation

The most important source of error when approximating (7) is the lack of a
closed form expression for $g_{t}$ as it involves an intractable integral. In
the context of inverse problems, a simple approximation used by (Chung et al.,
2023; Song et al., 2023) is given by

$\displaystyle\textstyle g_{t}(x_{t})\approx g_{0}(\int
x_{0}p(x_{0}|x_{t})\mathrm{d}x_{0})$ $\displaystyle=\textstyle
g_{0}(\sqrt{1-\lambda_{t}}x_{t})$ $\displaystyle:=\hat{g}_{t}(x_{t}).$ (8)

This approximation is good when $t$ is close to $0$ or $T$ but crude
otherwise, as established by the following result.

######  Proposition 2.2.

Let $\pi(x)=\mathcal{N}(x;\mu,\sigma^{2})$ and $(\beta_{t})_{t\in[0,T]}$ be
any schedule satisfying
$\lim_{T\rightarrow\infty}\int_{0}^{T}\beta_{s}\differential s=\infty$.
Consider the following approximation of (7)

$\displaystyle\differential Z_{t}^{(T)}$
$\displaystyle=[-\beta_{T-t}Z_{t}^{(T)}+2\beta_{T-t}\nabla\log\hat{g}_{T-t}(Z_{t}^{(T)})]\differential
t$ (9) $\displaystyle+\sqrt{2\beta_{T-t}}\differential B_{t},\quad
Z_{0}^{(T)}\sim\mathcal{N}(0,1).$

Then $\lim_{T\to\infty}\mathbb{E}[Z_{T}^{(T)}]=\mu$ and
$\lim_{T\to\infty}\operatorname{Var}(Z_{T}^{(T)})=1$ for $\sigma=1$ and
otherwise

$\displaystyle\lim_{T\to\infty}\mathbb{E}[Z_{T}^{(T)}]=\frac{\mu}{1-\sigma^{2}}(1-e^{-(1/\sigma^{2}-1)}),$
(10)
$\displaystyle\lim_{T\to\infty}\operatorname{Var}(Z_{T}^{(T)})=\frac{1-e^{-2(1/\sigma^{2}-1)}}{2(1/\sigma^{2}-1)}.$
(11)

Hence, even without considering any time discretization error,
$\lim_{T\to\infty}\mathbb{E}[Z_{T}^{(T)}]\neq\mu$ and
$\lim_{T\to\infty}\operatorname{Var}(Z_{T}^{(T)})\neq\sigma^{2}$; the error
getting larger as $\sigma^{2}$ moves away from $1$. If we consider the target
$\mathcal{N}(\mu,\sigma^{2})^{\otimes d}$, a similar result holds along each
dimension. Therefore, the Kullback–Leibler divergence between the target and
the result of running (9) grows linearly with $d$. As a result an
exponentially increasing number of samples is needed to obtain an importance
sampling approximation of the target of reasonable relative variance
(Chatterjee & Diaconis, 2018).

##  3 Particle denoising diffusion sampler

In this section, we propose a particle method to correct the discrepancy
between the distribution outputted by the guided diffusion and the target. Let
$[P]=\\{1,...,P\\}$ for $P\in\mathbb{N}$. We first present the exact joint
distribution of $(X_{t_{k}})_{k\in\\{0,...,K\\}}$ with $t_{k}=k\delta$ for a
fixed step size $\delta$ for the diffusion (2) where $K=T/\delta$ is an
integer. We then make explicit the corresponding reverse-time Markov
transitions for this joint distribution and show how they can be approximated.
Finally we show how these approximations can be corrected to obtain consistent
estimates using SMC. Instead of writing $X_{t_{k}}$, we will write $X_{k}$ to
simplify notation. Similarly we will write $g_{k}$ for $g_{t_{k}}$,
$\lambda_{k}$ for $\lambda_{t_{k}}$ etc.

###  3.1 From continuous time to discrete time

For $k\in[K]$, let
$\alpha_{k}=1-\exp[-2\int_{(k-1)\delta}^{k\delta}\beta_{s}\differential s]$.
The joint distribution of $X_{0:K}=(X_{0},X_{1},...,X_{K})$ under (2)
satisfies

$\textstyle\pi(x_{0:K})=\pi(x_{0})\prod\nolimits_{k\in[K]}p(x_{k}|x_{k-1}),$
(12)

for
$p(x_{k}|x_{k-1})=\mathcal{N}(x_{k};\sqrt{1-\alpha_{k}}x_{k-1},\alpha_{k}\mathrm{I})$.
This implies in particular

$\pi(x_{k},x_{k+1})=\pi_{k}(x_{k})p(x_{k+1}|x_{k}).$ (13)

We also denote by $p(x_{0:K})$ the joint distribution of (2) initialized at
$p_{0}(x_{0})$, in this case the Markov process is stationary, i.e.
$p_{k}(x_{k})=p_{0}(x_{k})$, and reversible w.r.t. $p_{0}$.

From (5), it can be easily checked that the backward transitions of (12)
satisfy

$\displaystyle\pi(x_{k}|x_{k+1})=\frac{p(x_{k}|x_{k+1})g_{k}(x_{k})}{g_{k+1}(x_{k+1})}$
$\displaystyle\approx p(x_{k}|x_{k+1})\exp[\langle\nabla\log
g_{k+1}(x_{k+1}),x_{k}-x_{k+1}\rangle]$
$\displaystyle=\textstyle\mathcal{N}(x_{k};\sqrt{1-\alpha_{k+1}}x_{k+1}$
$\displaystyle\qquad+\alpha_{k+1}\nabla\log
g_{k+1}(x_{k+1}),\alpha_{k+1}\mathrm{I}),$ (14)

where we used $g_{k}\approx g_{k+1}$ and a Taylor expansion of $\log
g_{k+1}(x_{k})$ around $x_{k+1}$; both of which are reasonable when $\delta\ll
1$. Since $\alpha_{k+1}\approx 2\beta_{k+1}\delta$ and
$\sqrt{1-\alpha_{k+1}}\approx 1-\beta_{k+1}\delta$, this approximation of
$\pi(x_{k}|x_{k+1})$ corresponds to a form of exponential integrator for the
time-reversal (7).

While this discrete-time representation is interesting, it is clearly
typically impossible to exploit it to obtain exact samples from $\pi$ since,
like its continuous time counterpart (7), it relies on the intractable
potentials $(g_{k})_{k=1}^{K}$. In the following Section 3.2, given an
approximation $\hat{g}_{k}$ of $g_{k}$, we propose a particle mechanism to
sample exactly from $\pi$ as the number of particles goes to infinity.

###  3.2 From discrete time to particles

We use a particle method to sample from $\pi$. The key idea is to break the
difficult problem of sampling from $\pi$ into a sequence of simpler
intermediate sampling problems. Ideally we would sample from $\pi_{K}$ first
then $\pi_{K-1},\pi_{K-2},...$ until $\pi_{0}=\pi$. Unfortunately this is not
possible as this requires knowing $g_{k}$. Suppose that we have an
approximation $\hat{g}_{k}$ of $g_{k}$ (for instance, the guidance
approximation defined in (2.3)). Inspired by (5), we sample instead for
$k\in[K]$ from the sequence of densities

$\textstyle\hat{\pi}_{k}(x_{k})\propto
p_{0}(x_{k})\hat{g}_{k}(x_{k}),\quad\hat{\mathcal{Z}}_{k}=\int
p_{0}(x_{k})\hat{g}_{k}(x_{k})\mathrm{d}x_{k},$ (15)

backward in time. At $k=K$, the function $g_{K}$ is almost constant so we
choose $\hat{g}_{K}\equiv 1$ and sampling from
$\hat{\pi}_{K}=\mathcal{N}(0,\mathrm{I})$ is easy.

Given particles approximately distributed according to $\hat{\pi}_{k+1}$, we
aim to obtain sample from $\hat{\pi}_{k}$. Drawing inspiration from (14), we
sample according to the proposal

$\hat{\pi}(x_{k}|x_{k+1}):=\mathcal{N}(x_{k};\sqrt{1-\alpha_{k+1}}x_{k+1}\\\
+\alpha_{k+1}\nabla\log\hat{g}_{k+1}(x_{k+1}),\alpha_{k+1}\mathrm{I}).$ (16)

We then reweight these pairs using the weights

$\displaystyle w_{k}(x_{k},x_{k+1})$
$\displaystyle:=\frac{\hat{\pi}_{k}(x_{k})p(x_{k+1}|x_{k})}{\hat{\pi}_{k+1}(x_{k+1})\hat{\pi}(x_{k}|x_{k+1})}$
$\displaystyle\propto\frac{\hat{g}_{k}(x_{k})p(x_{k}|x_{k+1})}{\hat{g}_{k+1}(x_{k+1})\hat{\pi}(x_{k}|x_{k+1})},$
(17)

where the last line follows from a simple application of (15). The first line
of (17) means that the $x_{k}$ marginal of the weighted system consistently
approximates $\hat{\pi}_{k}$. It should be noted that $w_{k}$ quantifies the
error in (16).

Finally, we resample these particles with weights proportional to (17). This
resampling operation allows us to focus the computational efforts on promising
regions of the space but some particles are replicated multiple times,
reducing the population diversity. Therefore, we then optionally perturb the
resampled particles using a MCMC kernel of invariant distribution
$\hat{\pi}_{k}$. The resulting Particle denoising diffusion sampler (PDDS) is
summarized in Algorithm 1.

Algorithm 1 Particle Denoising Diffusion Sampler

0: Schedule $(\beta_{t})_{t\in[0,T]}$ as in (2); Approximations
$(\hat{g}_{k})_{k=0}^{K}$ s.t. $\hat{g}_{0}=g_{0},\hat{g}_{K}=1$; Number of
particles $N$

Sample $X_{K}^{i}\overset{\mathrm{iid}}{\sim}\mathcal{N}(0,\mathrm{I})$ for
$i\in[N]$

Set $\hat{\mathcal{Z}}_{K}\leftarrow 1$ and $\omega_{K}^{i}\leftarrow 1/N$ for
$i\in[N]$

for $k=K-1,\ldots,0$ do

Move. Sample $\tilde{X}_{k}^{i}\sim\hat{\pi}(\cdot|X_{k+1}^{i})$ for $i\in[N]$
(see (16))

Weight. $\omega_{k}^{i}\leftarrow\omega_{k}(\tilde{X}_{k}^{i},X_{k+1}^{i})$
for $i\in[N]$ (see (17))

Set
$\hat{\mathcal{Z}}_{k}\leftarrow\hat{\mathcal{Z}}_{k+1}\times\frac{1}{N}\sum_{i\in[N]}\omega_{k}^{i}$

Normalize
$\omega_{k}^{i}\leftarrow\omega_{k}^{i}/\sum_{j\in[N]}\omega_{k}^{j}$

Resample.
$X_{k}^{1:N}\leftarrow\operatorname{resample}(\tilde{X}_{k}^{1:N},\omega_{k}^{1:N})$
(see Section 3.3)

MCMC (Optional). Sample $X_{k}^{i}\leftarrow\mathfrak{M}_{k}(X_{k}^{i},\cdot)$
for $i\in[N]$ using a $\hat{\pi}_{k}$-invariant MCMC kernel
$\mathfrak{M}_{k}$.

end for

Estimates $\hat{\pi}^{N}=\frac{1}{N}\sum_{i\in[N]}\delta_{X_{0}^{i}}$ of
$\pi$, $\hat{\mathcal{Z}}^{N}_{0}$ of $\mathcal{Z}$

###  3.3 Algorithm settings

Reparameterization. In practice, we would like to have $p_{0}$ to be such that
$g_{0}$ is bounded or has bounded moments. To achieve this, it can be
desirable to obtain a variational approximation $\mathcal{N}(x;\mu,\Sigma)$ of
$\pi$ then do the change of variables $x^{\prime}=\Sigma^{-1/2}(x-\mu)$,
sample in this space using PDDS before mapping the samples back using
$x=\mu+\Sigma^{1/2}x^{\prime}$.

Resampling. The idea of resampling is to only propagate particles in promising
regions of the state space. Given $N$ particles $\tilde{X}_{k}^{1:N}$ and $N$
weights $\omega_{k}^{1:N}$ summing to $1$, resampling selects $N$ output
particles $X_{k}^{1:N}$ such that, for any function
$\varphi:\mathbb{R}^{d}\to\mathbb{R}$,

$\mathbb{E}\Bigl{[}\tfrac{1}{N}{\textstyle\sum}_{i\in[N]}\varphi(X_{k}^{i})|\tilde{X}_{k}^{1:N},\omega_{k}^{1:N}\Bigr{]}={\textstyle\sum}_{i\in[N]}\omega_{k}^{i}\varphi(\tilde{X}_{k}^{i}).$

Popular schemes satisfying this identity are multinomial, stratified,
residual, and systematic resampling (Douc & Cappé, 2005).

We employ systematic resampling in all our simulations as it provides the
lowest variance estimates.

Resampling can however reduce the diversity of the particle population by
introducing identical particles in the output. As such, a popular recipe is to
trigger resampling at time $k$ only when the Effective Sample Size (ESS), a
measure of particle diversity defined by
$(\sum_{i\in[N]}(\omega^{i}_{k})^{2})^{-1}$, is below a certain threshold (Del
Moral et al., 2012; Dai et al., 2022). This is implemented using Algorithm 3
presented in Appendix A.

MCMC kernel. We want to design an MCMC kernel of invariant distribution
$\hat{\pi}_{k}(x_{k})$ defined in (15). As we have access to
$\nabla\log\hat{\pi}_{k}(x_{k})=-x_{k}+\nabla\log\hat{g}_{k}(x_{k})$, we can
use a Metropolis-adjusted Langevin algorithm (MALA) or Hamiltonian Monte
Carlo; e.g. MALA considers a proposal

$\textstyle
x^{\star}_{k}=x_{k}+\gamma\nabla\log\hat{\pi}_{k}(x_{k})+\sqrt{2\gamma}\epsilon,\quad\epsilon\sim\mathcal{N}(0,\mathrm{I}),$

for a step size $\gamma$. This proposal is accepted with probability

$\min\biggl{\\{}1,\frac{\hat{\pi}_{k}(x_{k}^{\star})\mathcal{N}(x_{k};x^{\star}_{k}+\gamma\nabla\log\hat{\pi}_{k}(x^{\star}_{k}),2\gamma\mathrm{I})}{\hat{\pi}_{k}(x_{k})\mathcal{N}(x^{\star}_{k};x_{k}+\gamma\nabla\log\hat{\pi}_{k}(x_{k}),2\gamma\mathrm{I})}\biggr{\\}}.$

###  3.4 Theoretical results

##### Fixed number of discretization steps.

We show below that the estimates $\hat{\mathcal{Z}}^{N}_{0}$ and
$\pi^{N}f=\tfrac{1}{N}\sum_{i=1}^{N}f(X^{i}_{0})$ returned by Algorithm 1 both
satisfy a central limit theorem. This follows from standard SMC theory (Del
Moral, 2004; Webber, 2019).

######  Proposition 3.1.

Assume that $\mathbb{E}[w_{k}(X_{k},X_{k+1})^{2}]<\infty$ where the
expectation is w.r.t. $\hat{\pi}(x_{k+1})\hat{\pi}(x_{k}|x_{k+1})$ and that
$\int\pi_{k}(x)(g_{k}/\hat{g}_{k})(x)\mathrm{d}x<\infty$. Then
$\hat{\mathcal{Z}}^{N}_{0}$ is an unbiased estimate of $\mathcal{Z}$ and has
finite variance. If multinomial resampling is used at every step,
$\sqrt{N}(\hat{\mathcal{Z}}^{N}_{0}/\mathcal{Z}-1)$ is asymptotically normal
with asymptotic variance

$\displaystyle\sigma_{K}^{2}=\chi^{2}(\pi_{K}||\hat{\pi}_{K})+$
$\displaystyle\sum_{k=0}^{K-1}\chi^{2}(\pi_{k}(x_{k})\pi(x_{k}|x_{k+1})||\hat{\pi}_{k}(x_{k})\hat{\pi}(x_{k}|x_{k+1})),$

with $\chi^{2}(\cdot||\cdot)$ the chi-squared divergence between two
distributions. Moreover, for any bounded function $f$, we also have asymptotic
normality of $\sqrt{N}(\hat{\pi}^{N}f-\pi f)$.

The finiteness assumptions on $\mathbb{E}[w_{k}(X_{k},X_{k+1})^{2}]$ and
$\int\pi_{k}(x)(g_{k}/\hat{g}_{k})(x)\mathrm{d}x$ require that $\hat{g}_{k}$
is not too far from $g_{k}$ but still allow enough freedom in the choice of
$\hat{g}_{k}$. The expression for $\sigma_{K}^{2}$ suggests that choosing
$\hat{g}_{k}$ close to $g_{k}$ will reduce the asymptotic variance as (16)
will better approximate (14).

##### Infinitely fine discretization limit.

We now investigate the performance of PDDS as the number of discretization
time steps goes to infinity. Even when all the weights are equal, multinomial
resampling still kills over a third of particles on average (Chopin &
Papaspiliopoulos, 2020). Hence a fine discretization with repeated
applications of multinomial resampling leads to the total collapse of the
particle approximation. This has been formalized in the continuous-time
setting (Chopin et al., 2022). In our case, it can be readily checked that
indeed $\sigma_{K}^{2}\to\infty$ as $K\to\infty$ in general. This justifies
using a more sophisticated resampling strategy. Indeed, when using the sorted
stratified resampling strategy of Gerber et al. (2019), the following results
show that our particle approximations remain well behaved as
$K\rightarrow\infty$.

######  Proposition 3.2.

Consider the setting of Proposition 3.1 with sorted stratified resampling.
Then there exists a sequence of sets $(B_{N})$ such that $\mathbb{P}(B_{N})\to
1$ and

$\limsup_{N\to\infty}\mathbb{E}[N(\hat{\mathcal{Z}}^{N}_{0}/\mathcal{Z}-1)^{2}\mathbbm{1}_{B_{N}}]\leq\zeta^{2}_{K}$

with

$\displaystyle\zeta^{2}_{K}:=\chi^{2}(\pi_{K}||\hat{\pi}_{K})+$
$\displaystyle\sum_{k=0}^{K-1}\int\frac{\pi_{k+1}(x_{k+1})^{2}}{\hat{\pi}_{k+1}(x_{k+1})}\chi^{2}(\pi(x_{k}|x_{k+1})||\hat{\pi}(x_{k}|x_{k+1}))\differential
x_{k+1}.$

The next result bounds the limit when $K\to\infty$, i.e.
$\delta=T/K\rightarrow 0$

######  Proposition 3.3.

Under the setting of Proposition 3.2, assume that $\beta_{t}\equiv 1$ and that
the target distribution satisfies the regularity conditions in Appendix C.4.1.
Let $\bar{g}_{t}:=g_{t}/\mathcal{Z}$ and $\tilde{g}_{t}:=\hat{g}_{t}/\int
p_{0}(x)\hat{g}_{t}(x)\differential x$. Assume further that the approximations
$\hat{g}_{t},\bar{g}_{t},\tilde{g}_{t}$ satisfy
$C_{1}^{-1}\leq\bar{g}_{t}/\tilde{g}_{t}\leq C_{1}$ and
$\norm{\nabla\log\hat{g}_{t}(x_{t})}\leq C_{2}(1+\norm{x_{t}})$ for some
$C_{1},C_{2}\geq 0$. Then

$\displaystyle\textstyle\limsup_{K\to\infty}\zeta^{2}_{K}\leq\chi^{2}(\pi_{T}||\hat{\pi}_{T})+$
$\displaystyle\textstyle
2\int_{0}^{T}\mathbb{E}_{X_{t}\sim\pi_{t}}\left[\frac{\bar{g}_{t}}{\tilde{g}_{t}}(X_{t})\norm{\nabla\log
g_{t}(X_{t})-\nabla\log\hat{g}_{t}(X_{t})}^{2}\right]\differential t.$

This result highlights precisely how the asymptotic error depends on the
quality of the estimates of $g_{t}$ and its logarithmic gradient. Let
$\pi^{\hat{g}}(\differential x_{[0,T]})$ be the path measure induced by
running (7) with some approximation $\hat{g}_{t}$, i.e. the distribution of
$(Y_{t})_{t\in[0,T]}$ given by (7). If importance sampling were directly used
to correct between $\pi^{\hat{g}}(\differential x_{[0,T]})$ and
$\pi(\differential x_{[0,T]})$, the asymptotic error for
$\hat{\mathcal{Z}}^{N}_{0}/\mathcal{Z}$ would be equal to
$\chi^{2}(\pi||\pi^{\hat{g}})$ which is greater than
$\exp{\mathrm{KL}(\pi|\pi^{\hat{g}})}-1$. In contrast, ignoring the negligible
first term $\chi^{2}(\pi_{T}||\hat{\pi}_{T})$, the error in this proposition
is upper bounded by $2C_{1}\mathrm{KL}(\pi|\pi^{\hat{g}})$. This shows how
PDDS helps reduce the error of naive importance sampling.

##  4 Learning potentials via score matching

The previous theoretical results establish the consistency of PDDS estimates
for any reasonable approximation $\hat{g}_{k}$ of $g_{k}$ but also show that
better approximations lead to lower Monte Carlo errors. We show here how to
use the approximation of $\pi$ outputted by PDDS to learn a better neural
network (NN) approximation $\hat{g}_{\theta}(k,x_{k})$ of the potential
functions $g_{k}$ by leveraging score matching ideas. Once we have learned
those approximations, we can then run again PDDS (Algorithm 1) with the new
learned potentials. In our experiments (Section 6), we demonstrate that this
can drastically improve results in challenging scenarios.

###  4.1 Different score identities

From (5), we have the identity

$\nabla\log g_{k}(x_{k})=x_{k}+\nabla\log\pi_{k}(x_{k}).$ (18)

Hence, if we obtain an approximation $\nabla\log\pi_{\theta}(k,x_{k})$ of
$\nabla\log\pi_{k}(x_{k})$, then we get an approximation
$\log\hat{g}_{\theta}(k,x_{k})=\tfrac{1}{2}||x_{k}||^{2}+\log\hat{\pi}_{\theta}(k,x_{k})$
of $\log g_{k}(x_{k})$.

To learn the score, we rely on the following result.

######  Proposition 4.1.

The score satisfies the standard Denoising Score Matching (DSM) identity

$\textstyle\nabla\log\pi_{k}(x_{k})=\int\nabla\log
p(x_{k}|x_{0})~{}\pi(x_{0}|x_{k})\differential x_{0}.$ (19)

Moreover, if
$\int\norm{\nabla\pi(x_{0})}e^{-\eta\norm{x_{0}}^{2}}\differential
x_{0}<\infty,\forall\eta>0$, the Novel Score Matching (NSM) identity

$\textstyle\nabla\log\pi_{k}(x_{k})=\kappa_{k}\int\nabla\log
g_{0}(x_{0})~{}\pi(x_{0}|x_{k})\differential x_{0}-x_{k}$ (20)

holds for $\kappa_{k}=\sqrt{1-\lambda_{k}}$ and $\nabla\log
g_{0}(x_{0})=\nabla\log\pi_{0}(x_{0})+x_{0}$. Hence, we can approximate the
score by minimizing one of the two following loss functions

$\displaystyle\ell_{\textup{DSM}}(\theta)=\sum_{k\in[K]}\mathbb{E}\norm{\nabla\log\hat{\pi}_{\theta}(k,X_{k})-\nabla\log
p(X_{k}|X_{0})}^{2},$
$\displaystyle\ell_{\textup{NSM}}(\theta)=\sum_{k\in[K]}\mathbb{E}||\nabla\log\hat{\pi}_{\theta}(k,X_{k})+X_{k}$
$\displaystyle\hskip 71.13188pt-\kappa_{k}\nabla\log g_{0}(X_{0})||^{2},$

where the first loss is applicable if $\pi$ has finite second moment and the
second loss is applicable if additionally
$\mathbb{E}_{\pi}[\norm{\nabla\log\pi(X)}^{2}]<\infty$. All the expectations
are taken w.r.t. $\pi(x_{0})p(x_{k}|x_{0})$. For expressive neural networks,
both $\ell_{\textup{DSM}}(\theta)$ and $\ell_{\textup{NSM}}(\theta)$ are such
that $\nabla\log\hat{\pi}_{\theta}(k,x_{k})=\nabla\log\pi_{k}(x_{k})$ at the
minimizer.

The benefit of this novel loss is that it is much better behaved compared to
the standard denoising score matching loss when $\delta\ll 1$ as established
below, this is due to the fact that the variance of the terms appearing in
$\ell_{\textup{DSM}}$ for $k$ close to zero become very large. This loss is
not applicable to generative modeling as $\pi$ is only available through
samples.

###  4.2 Benefits of alternative score matching identity

We establish here formally the benefits of NSM over DSM. NSM score matching
prevents variance blow up as $k$ is close to time $0$ and $\delta\rightarrow
0$. This is more elegantly formalized in continuous-time. In this case, the
renormalized loss functions $\ell_{\textup{DSM}}$ and $\ell_{\textup{NSM}}$
introduced in Proposition 4.1 become for
$s_{\theta}(t,x_{t})=\nabla\log\hat{\pi}_{\theta}(t,x_{t})$

$\displaystyle\textstyle\ell_{\textup{DSM}}(\theta)=\textstyle\int^{T}_{0}\mathbb{E}\norm{s_{\theta}(t,X_{t})-\nabla\log
p(X_{t}|X_{0})}^{2}\differential t,$
$\displaystyle\textstyle\ell_{\textup{NSM}}(\theta)=\textstyle\int^{T}_{0}\mathbb{E}||s_{\theta}(t,X_{t})+X_{t}-\kappa_{t}\nabla\log
g_{0}(X_{0})||^{2}\differential t,$

where the expectations are w.r.t. $\pi(x_{0})p(x_{t}|x_{0})$. Unbiased
estimates of the gradient of these losses
$\hat{\nabla\ell}_{\textup{DSM}}(\theta)$ and
$\hat{\nabla\ell}_{\textup{NSM}}(\theta)$ are obtained by computing the
gradient w.r.t. $\theta$ of the argument within the expectation at a sample
$\tau\sim\operatorname{Unif}[0,T]$, $X_{0}\sim\pi$ and $X_{\tau}\sim
p(x_{\tau}|X_{0})$.

The following result clarifies the advantage of NSM score matching.

######  Proposition 4.2.

Let $d=1$, $\beta_{t}\equiv 1$, and
$\pi(x_{0})=\mathcal{N}(x_{0};\mu,\sigma^{2})$. Suppose that the score network
$s_{\theta}(t,x_{t})$ is a continuously differentiable function jointly on its
three variables. Moreover, assume that
$\absolutevalue{s}+\norm{\nabla_{\theta}s}\leq
C(1+|\theta|+|t|+|x_{t}|)^{\alpha}$ for some $C,\alpha>0$; and that
$\mathbb{E}_{\pi}[\norm{\nabla_{\theta}s_{\theta}(0,X_{0})}^{2}]>0$ for all
$\theta$. Then, for all $\theta$, $\ell_{\textup{DSM}}(\theta)$ and
$\mathbb{E}[\hat{||\nabla\ell}_{\textup{DSM}}(\theta)||^{2}]$ are infinite
whereas $\ell_{\textup{NSM}}(\theta)$ and
$\mathbb{E}[||\hat{\nabla\ell}_{\textup{NSM}}(\theta)||^{2}]$ are finite.

###  4.3 Neural network parametrization

Contrary to standard practice for generative modeling (Ho et al., 2020; Song
et al., 2021), we parameterize a function and not a vector field. This is
necessary to run PDDS (Algorithm 1) as it relies on potentials to weight
particles.

We will use a parameterization of the form
$\log\hat{\pi}_{\theta}(k,x_{k})=\log\hat{g}_{\theta}(k,x_{k})-\tfrac{1}{2}||x||^{2}$
where we parametrize $\hat{g}_{\theta}$ using two neural networks $r_{\eta}$
and $N_{\gamma}$ such that $\theta=(\eta,\gamma)$. The network $r_{\eta}$
returns a scalar whereas $N_{\gamma}$ returns a vector in $\mathbb{R}^{d}$.
The precise expression is given by

$\log\hat{g}_{\theta}(k,x_{k})=\left[r_{\eta}(k)-r_{\eta}(0)\right]\langle
N_{\gamma}(k,x_{k}),x_{k}\rangle+\\\
+\left[1-r_{\eta}(k)+r_{\eta}(0)\right]\log g_{0}(\sqrt{1-\lambda_{k}}x_{k}).$

This parametrization takes advantage of the simple approximation (2.3) with
which it coincides at time $0$. Zhang & Chen (2022) used a similar
parametrization to incorporate gradient information in their control policy.
Although $\log\hat{g}_{\theta}(k,x_{k})$ is a scalar, we deliberately let
$N_{\gamma}(k,x_{k})$ return a vector which is then scalar-multiplied with
$x_{k}$. This is usually done in the literature when the scalar potential
instead of just its gradient is learned (see e.g. Salimans & Ho, 2021) and
helps improve model expressiveness.

###  4.4 Training the neural network

Both $\ell_{\textup{DSM}}(\theta)$ and $\ell_{\textup{NSM}}(\theta)$ can be
written as $\sum_{k\in[K]}\mathbb{E}\left[\ell(\theta,k,X_{0},X_{k})\right]$
for appropriately defined local loss functions $\ell$. Algorithm 2 describes
the gradient updates according to the batch size $B$.

Algorithm 2 Potential neural network training

0: Particle approx. $\hat{\pi}^{N}$ outputted by Algorithm 1; Potential NN
$\hat{g}_{k}(\theta,x_{k})$; Initialization $\theta_{0}$; Local loss functions
$\ell(\theta,k,x_{0},x_{k})$; Batch size $B$; Number of gradient updates
$N_{\mathrm{up}}$; Learning rate $\eta>0$.

for $i=1,2,\ldots,N_{\mathrm{up}}$ do

for $b=1,2,\ldots,B$ do

Sample $X^{b}_{0}\sim\hat{\pi}^{N}(\cdot)$; $k_{b}\sim\operatorname{Unif}[K]$;

Sample $X_{k_{b}}\sim p_{k_{b},0}(\cdot|X_{0}^{b})$

end for

$\theta_{i}:=\theta_{i-1}-\tfrac{\eta}{B}\sum_{b\in[B]}\nabla_{\theta}\ell(\theta_{i-1},k_{b},X_{0}^{b},X_{k_{b}})$

end for

Potential $\hat{g}_{\theta_{N_{\mathrm{up}}}}(k,x_{k})$

In practice, we first run Algorithm 1 with a simple approximation such as
(2.3). We then use Algorithm 2 to learn $\hat{g}_{\theta}(k,x_{k})$ and can
then execute Algorithm 1 again with
$\hat{g}_{k}(x_{k})=\hat{g}_{\theta}(k,x_{k})$ to obtain lower variance
estimates of $\pi$ and $\mathcal{Z}$. We can further refine the approximation
of $g_{k}$ using Algorithm 2. In practice, we found that one or two iterations
are sufficient.

##  5 Related work

SMC samplers (Del Moral et al., 2006; Dai et al., 2022) are a general
methodology to sample sequentially from a sequence of distributions. It relies
on the notion of forward and backward kernels in order to move from one
distribution to another. PDDS can be cast in this framework where the forward
kernel is chosen to be the forward noising diffusion and the backward kernel
the approximate time-reversal. The novelty of our work is the exploitation of
the special structure induced by such a choice to come up with efficient
backward kernel estimates. Other standard methods include Annealed Importance
Sampling (Neal, 2001) and Parallel Tempering (Geyer, 1991). Unlike our work,
all these methods approximate a sequence of tempered versions of the target.
While they are standard, it is also well-known that tempering strategies can
exhibit poor performance for multimodal targets (Woodard et al., 2009; Tawn et
al., 2020; Syed et al., 2022) as tempering can change dramatically the masses
of distribution modes depending on their widths. Adding noise does not suffer
from this issue (Máté & Fleuret, 2023). It only perturbs the distribution
locally, preserving the weights of the modes; see Figure 1 for an
illustration.

Figure 1: Tempered (top) and noised (bottom) sequences of distributions for
the target $\pi(x)=0.8\mathcal{N}(x;-4,0.5^{2})+0.2\mathcal{N}(x;4,1)$. The
tempered sequence follows
$\pi_{t}(x)\propto\pi(x)^{(1-\eta_{t})}\phi(x)^{\eta_{t}}$ where $\phi$ is the
standard normal and $\eta_{t}$ increases from $0$ to $1$. The noising sequence
follows the forward diffusion in Equation 2. Red dots indicate the position of
modes in the original target. The tempered sequence suffers from mode
switching, i.e. the low mass large width mode becomes dominant across the
tempered path. The noised sequence does not suffer from this problem.

Our sampler can also be interpreted as sampling from a sequence of
distributions using so-called twisted proposals. The general framework for
approximating these twisted proposals using SMC has been considered in
Guarniero et al. (2017); Heng et al. (2020); Lawson et al. (2022). In
particular, Wu et al. (2023) and Cardoso et al. (2024) recently apply such
ideas to conditional simulation in generative modeling given access to a
pretrained score network. Wu et al. (2023) rely on the simple approximation
(2.3) and do not quantify its error, while Cardoso et al. (2024) use a
different proposal kernel which leverages the structure of a linear inverse
problem. Our setup is here different as we consider general Monte Carlo
sampling problems. We do not rely on any pretrained network and refine our
potential approximations using an original loss function. We additionally
provide theoretical results in particular when the discretization time step
goes to zero.

##  6 Experimental results

###  6.1 Normalizing constant estimation

Figure 2: Log normalizing constant estimation results for our method with one
and two iterations of potential approximation (PDDS1 and PDDS2 respectively),
compared with SMC, CRAFT, DDS and PIS. Dotted black represents analytic ground
truth where available, otherwise long-run SMC. Variation is displayed over
both training and sampling seeds (2000 total). The y-axes on Sonar and LGCP
have been cropped and outliers (present in all methods) removed for clarity,
uncurated samples are presented in Section D.4.

We evaluate here the quality of normalizing constant estimates produced by
PDDS on a variety of sampling tasks. We consider two synthetic target
densities and two posterior distributions, with results on additional targets
included in Section D.4. The synthetic targets are `Mixture`, a 2-dimensional
mixture of Gaussian distributions with separated modes (Arbel et al., 2021)
and `Funnel`, a 10-dimensional target displaying challenging variance
structure (Neal, 2003). The posterior distributions are `Sonar`, a logistic
regression posterior fitted to the Sonar (61-dimensional) dataset and `LGCP`,
a log Gaussian Cox process (Møller et al., 1998) modelling the rate parameter
of a Poisson point process on a $40\times 40=1600$-point grid, fitted to the
Pines dataset. Precise specification of these targets can be found in Section
D.1.

We compare our performance to a selection of strong baselines. We consider two
annealing algorithms, firstly an SMC sampler (Del Moral et al., 2006) with HMC
kernels and secondly CRAFT (Matthews et al., 2022), which uses normalizing
flows to transport particles at each step of an SMC algorithm. We also
consider two diffusion-based sampling methods, Path Integral Sampler (PIS)
(Zhang & Chen, 2022) and Denoising Diffusion Sampler (DDS) (Vargas et al.,
2023). As mentioned in Section 3.3, we reparameterize the target using a
variational approximation for all methods. We include full details of
hyperparameter and optimizer settings, run times and experimental procedures
in Section D.2.

We present the normalizing constant estimation results in Figure 2. We report
results from our method with both one and two iterations of potential
approximation (PDDS1 and PDDS2 respectively). PDDS1 is given the same training
budget as PIS and DDS. Considering the posterior sampling tasks, we uniformly
outperform both diffusion-based approaches (PIS and DDS) in terms of
estimation bias and variance. We find that in the `Sonar` task CRAFT
outperforms our approach in few step regimes but we close this gap and
sometimes outperform CRAFT given a larger number of steps, while the reverse
is true on the `LGCP` task. Considering the synthetic target densities, our
approach similarly outperforms the diffusion-based approaches. Our approach
produces higher variance estimates than CRAFT in low step regimes, but
produces lower variance estimates in larger step regimes. CRAFT clearly
outperforms all competing methods on the synthetic `Funnel` task.

###  6.2 Sample quality

We also visually assess the quality of samples from each method in Figure 3.
We choose the multi-modal `Mixture` task in 2-dimensions for ease of
visualization. Unsurprisingly we find that both PIS and DDS do not capture all
6 modes of the distribution due to the mode-seeking behaviour of the reverse
KL objective. Samples from our approach appear of similar quality to those of
CRAFT and SMC.

Figure 3: Samples from each method using 4 steps on the Gaussian Mixture task.

###  6.3 Iterations of potential approximation

Figure 4: ESS curves for our method after 0, 1, and 2 iterations of Algorithm
2 on the Sonar task with 200 steps.

We demonstrate here the improvement in normalizing constant estimates due to
our iterative potential approximation scheme. Figure 4 shows the improvement
in ESS and reduction in number of resampling steps. In the top pane of Figure
5, we see that the simple potential approximation (2.3) provides very poor
normalizing constant estimates. Algorithm 2 considerably improves upon this
performance.

In the second pane of Figure 5, we show the evolution of
$\hat{\mathcal{Z}}^{N}_{0}$ during training for each of the diffusion-based
approaches. While PDDS initially falls below PIS and DDS due to the simple
initial approximation, we quickly recover and exceed each of these methods
after around $7\times 10^{7}$ density evaluations. We also observe that PDDS
improves incrementally due to the iterative potential approximation.

Figure 5: Top: normalizing constant estimation after 0, 1 and 2 iterations of
Algorithm 2 on the Sonar task. Bottom: evolution of normalizing constant
estimates during training, measured in terms of number of target density
evaluations, showing one realization on the Sonar task with 16 steps.

##  7 Discussion

This paper contributes to the growing literature on the use of denoising
diffusion ideas for Monte Carlo sampling (Berner et al., 2022; McDonald &
Barron, 2022; Vargas et al., 2023; Huang et al., 2024; Richter et al., 2024).
It proposes an original iterative SMC algorithm which provides an unbiased
estimate of the normalizing constant for any finite number of particles by
leveraging an original score matching technique. This algorithm also provides
asymptotically consistent estimates of the normalizing constant and of
expectations with respect to the target. One limitation of PDDS is that it
practically relies on $g_{0}(x)$ being a well-behaved potential function.
While our approach using a variational approximation to guide a
reparameterization of the target has been effective in our examples, more
sophisticated techniques might have to be implemented (Hoffman et al., 2019).

##### Impact Statement.

Sampling is a ubiquitous problem. Therefore, PDDS can be applied to a wide
range of applications. While we have obtained limit theorems for this scheme,
it is important to exercise caution if the output were to be used for decision
making as we practically only use a finite number of particles.

## References

  * Arbel et al. (2021) Arbel, M., Matthews, A., and Doucet, A.  Annealed flow transport Monte Carlo.  In _International Conference on Machine Learning_ , 2021.
  * Berner et al. (2022) Berner, J., Richter, L., and Ullrich, K.  An optimal control perspective on diffusion-based generative modeling.  In _NeurIPS 2022 Workshop on Score-Based Methods_ , 2022.
  * Bradbury et al. (2018) Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., and Zhang, Q.  JAX: composable transformations of Python+NumPy programs, 2018\.
  * Cardoso et al. (2024) Cardoso, G., Idrissi, Y. J. E., Corff, S. L., and Moulines, E.  Monte Carlo guided diffusion for Bayesian linear inverse problems.  In _International Conference on Learning Representations_ , 2024.
  * Cattiaux et al. (2023) Cattiaux, P., Conforti, G., Gentil, I., and Léonard, C.  Time reversal of diffusion processes under a finite entropy condition.  _Annales de l’Institut Henri Poincare (B) Probabilites et statistiques_ , 59(4):1844–1881, 2023.
  * Chatterjee & Diaconis (2018) Chatterjee, S. and Diaconis, P.  The sample size required in importance sampling.  _Ann. Appl. Probab._ , 28(2):1099–1135, 2018\.
  * Chopin & Papaspiliopoulos (2020) Chopin, N. and Papaspiliopoulos, O.  _An Introduction to Sequential Monte Carlo_.  Springer Series in Statistics. 2020.
  * Chopin et al. (2022) Chopin, N., Singh, S. S., Soto, T., and Vihola, M.  On resampling schemes for particle filters with weakly informative observations.  _Ann. Stat._ , 50(6):3197–3222, 2022.  ISSN 0090-5364.
  * Chung et al. (2023) Chung, H., Kim, J., Mccann, M. T., Klasky, M. L., and Ye, J. C.  Diffusion posterior sampling for general noisy inverse problems.  In _The Eleventh International Conference on Learning Representations_ , 2023.
  * Corso et al. (2023) Corso, G., Xu, Y., De Bortoli, V., Barzilay, R., and Jaakkola, T.  Particle guidance: non-iid diverse sampling with diffusion models.  _arXiv preprint arXiv:2310.13102_ , 2023.
  * Dai et al. (2022) Dai, C., Heng, J., Jacob, P. E., and Whiteley, N.  An invitation to sequential monte carlo samplers.  _Journal of the American Statistical Association_ , 117(539):1587–1600, 2022.
  * De Bortoli et al. (2021) De Bortoli, V., Thornton, J., Heng, J., and Doucet, A.  Diffusion Schrödinger bridge with applications to score-based generative modeling.  In _Advances in Neural Information Processing Systems_ , 2021.
  * Del Moral (2004) Del Moral, P.  _Feynman-Kac Formulae: Genealogical and Interacting Particle Approximations_.  Springer, 2004.
  * Del Moral et al. (2006) Del Moral, P., Doucet, A., and Jasra, A.  Sequential Monte Carlo samplers.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 68(3):411–436, 2006.
  * Del Moral et al. (2012) Del Moral, P., Doucet, A., and Jasra, A.  On adaptive resampling strategies for sequential Monte Carlo methods.  _Bernoulli_ , 18(1):252–278, 2012.
  * Douc & Cappé (2005) Douc, R. and Cappé, O.  Comparison of resampling schemes for particle filtering.  In _Proceedings of the 4th International Symposium on Image and Signal Processing and Analysis_ , pp. 64–69. IEEE, 2005.
  * Doucet et al. (2001) Doucet, A., De Freitas, N., and Gordon, N. J.  _Sequential Monte Carlo Methods in Practice_.  Information Science and Statistics. New York, NY: Springer, New York, 2001.
  * Gerber et al. (2019) Gerber, M., Chopin, N., and Whiteley, N.  Negative association, ordering and convergence of resampling methods.  _Ann. Stat._ , 47(4):2236–2260, 2019.  ISSN 0090-5364.
  * Geyer (1991) Geyer, C.  Markov chain Monte Carlo maximum likelihood.  In _Computing science and statistics: Proceedings of 23rd Symposium on the Interface Interface Foundation, Fairfax Station, 1991_ , pp. 156–163, 1991.
  * Guarniero et al. (2017) Guarniero, P., Johansen, A. M., and Lee, A.  The iterated auxiliary particle filter.  _Journal of the American Statistical Association_ , 112(520):1636–1647, 2017.
  * Haussmann & Pardoux (1986) Haussmann, U. G. and Pardoux, E.  Time reversal of diffusions.  _The Annals of Probability_ , 14(3):1188–1205, 1986.
  * Heng et al. (2020) Heng, J., Bishop, A. N., Deligiannidis, G., and Doucet, A.  Controlled sequential Monte Carlo.  _The Annals of Statistics_ , 48(5):2904–2929, 2020.
  * Ho et al. (2020) Ho, J., Jain, A., and Abbeel, P.  Denoising diffusion probabilistic models.  In _Advances in Neural Information Processing Systems_ , 2020.
  * Hoffman et al. (2019) Hoffman, M., Sountsov, P., Dillon, J. V., Langmore, I., Tran, D., and Vasudevan, S.  Neutra-lizing bad geometry in Hamiltonian Monte Carlo using neural transport.  _arXiv preprint arXiv:1903.03704_ , 2019.
  * Huang et al. (2024) Huang, X., Dong, H., Hao, Y., Ma, Y., and Zhang, T.  Reverse diffusion Monte Carlo.  In _International Conference on Learning Representations_ , 2024.
  * Hyvärinen (2005) Hyvärinen, A.  Estimation of non-normalized statistical models by score matching.  _The Journal of Machine Learning Research_ , 6:695–709, 2005.
  * Kingma & Ba (2015) Kingma, D. and Ba, J.  Adam: A method for stochastic optimization.  In _ICLR_ , 2015.
  * Kloeden & Platen (1992) Kloeden, P. E. and Platen, E.  _Numerical solution of stochastic differential equations_ , volume 23 of _Appl. Math. (N. Y.)_.  Berlin: Springer-Verlag, 1992.  ISBN 3-540-54062-8.
  * Lai et al. (2022) Lai, C.-H., Takida, Y., Murata, N., Uesaka, T., Mitsufuji, Y., and Ermon, S.  Regularizing score-based models with score fokker-planck equations.  _arXiv preprint arXiv:2210.04296_ , 2022.
  * Lawson et al. (2022) Lawson, D., Raventós, A., Linderman, S., et al.  Sixo: Smoothing inference with twisted objectives.  _Advances in Neural Information Processing Systems_ , 35:38844–38858, 2022.
  * Liptser & Shiryayev (1977) Liptser, R. S. and Shiryayev, A. N.  _Statistics of random processes. I. General theory. Translated by A. B. Aries_ , volume 5 of _Appl. Math. (N. Y.)_.  Springer, New York, 1977.
  * Máté & Fleuret (2023) Máté, B. and Fleuret, F.  Learning deformation trajectories of Boltzmann densities.  _arXiv preprint arXiv:2301.07388_ , 2023.
  * Matthews et al. (2022) Matthews, A. G. D. G., Arbel, M., Rezende, D. J., and Doucet, A.  Continual repeated annealed flow transport Monte Carlo.  In _International Conference on Machine Learning_ , 2022.
  * McDonald & Barron (2022) McDonald, C. J. and Barron, A. R.  Proposal of a score based approach to sampling using Monte Carlo estimation of score and oracle access to target density.  In _NeurIPS 2022 Workshop on Score-Based Methods_ , 2022.
  * Møller et al. (1998) Møller, J., Syversveen, A. R., and Waagepetersen, R. P.  Log Gaussian Cox processes.  _Scandinavian Journal of Statistics_ , 25(3):451–482, 1998.
  * Neal (2001) Neal, R. M.  Annealed importance sampling.  _Statistics and Computing_ , 11(2):125–139, 2001\.
  * Neal (2003) Neal, R. M.  Slice sampling.  _The Annals of Statistics_ , 31:705–767, 06 2003.
  * Nichol & Dhariwal (2021) Nichol, A. Q. and Dhariwal, P.  Improved denoising diffusion probabilistic models.  In _International Conference on Machine Learning_ , 2021.
  * Richter et al. (2024) Richter, L., Berner, J., and Liu, G.-H.  Improved sampling via learned diffusions.  In _International Conference on Learning Representations_ , 2024.
  * Salimans & Ho (2021) Salimans, T. and Ho, J.  Should EBMs model the energy or the score?  In _Energy Based Models Workshop-ICLR 2021_ , 2021.
  * Song et al. (2023) Song, J., Vahdat, A., Mardani, M., and Kautz, J.  Pseudoinverse-guided diffusion models for inverse problems.  In _International Conference on Learning Representations_ , 2023.
  * Song & Ermon (2019) Song, Y. and Ermon, S.  Generative modeling by estimating gradients of the data distribution.  _Advances in neural information processing systems_ , 32, 2019.
  * Song et al. (2021) Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B.  Score-based generative modeling through stochastic differential equations.  In _International Conference on Learning Representations_ , 2021.
  * Sountsov et al. (2020) Sountsov, P., Radul, A., and contributors.  Inference gym, 2020.  URL https://pypi.org/project/inference_gym.
  * Syed et al. (2022) Syed, S., Bouchard-Côté, A., Deligiannidis, G., and Doucet, A.  Non-reversible parallel tempering: A scalable highly parallel MCMC scheme.  _Journal of the Royal Statistical Society Series B_ , 84(2):321–350, 2022.
  * Tawn et al. (2020) Tawn, N. G., Roberts, G. O., and Rosenthal, J. S.  Weight-preserving simulated tempering.  _Statistics and Computing_ , 30(1):27–41, 2020\.
  * Vargas et al. (2023) Vargas, F., Grathwohl, W., and Doucet, A.  Denoising diffusion samplers.  In _International Conference on Learning Representations_ , 2023.
  * Vincent (2011) Vincent, P.  A connection between score matching and denoising autoencoders.  _Neural Computation_ , 23(7):1661–1674, 2011\.
  * Webber (2019) Webber, R. J.  Unifying sequential Monte Carlo with resampling matrices.  _arXiv preprint arXiv:1903.12583_ , 2019.
  * Woodard et al. (2009) Woodard, D. B., Schmidler, S. C., and Huber, M.  Sufficient conditions for torpid mixing of parallel and simulated tempering.  _Electronic Journal of Probability_ , 14:780–804, 2009\.
  * Wu et al. (2023) Wu, L., Trippe, B. L., Naesseth, C. A., Blei, D., and Cunningham, J. P.  Practical and asymptotically exact conditional sampling in diffusion models.  In _Thirty-seventh Conference on Neural Information Processing Systems_ , 2023.
  * Zhang & Chen (2022) Zhang, Q. and Chen, Y.  Path integral sampler: a stochastic control approach for sampling.  In _International Conference on Learning Representations_ , 2022.

## Appendix

The Appendix is organized as follows. In Appendix A we detail the PDDS
algorithm when adaptive resampling is used. In Appendix B we reformulate the
results of Webber (2019) in terms of chi-squared divergences. In Appendix C,
we expose our proofs. Finally, in Appendix D we present details and additional
work relating to our experiments.

##  Appendix A Particle Denoising Diffusion Sampler with Adaptive Resampling

Algorithm 3 Particle Denoising Diffusion Sampler with Adaptive Resampling

0: Schedule $(\beta_{t})_{t\in[0,T]}$ as in (2); Approximations
$(\hat{g}_{k})_{k=0}^{K}$ s.t. $\hat{g}_{0}=g_{0},\hat{g}_{K}=1$; Number of
particles $N$; ESS resampling threshold $\alpha$

Sample $X_{K}^{i}\overset{\mathrm{iid}}{\sim}\mathcal{N}(0,\mathrm{I})$ for
$i\in[N]$

Set $\hat{\mathcal{Z}}_{K}\leftarrow 1$ and $\omega_{K}^{i}\leftarrow 1/N$ for
$i\in[N]$

for $k=K-1,\ldots,0$ do

Move. Sample $\tilde{X}_{k}^{i}\sim\hat{\pi}(\cdot|X_{k+1}^{i})$ for $i\in[N]$
(see (16))

Weight.
$\omega_{k}^{i}\leftarrow\omega_{k}^{i}~{}\omega_{k}(\tilde{X}_{k}^{i},X_{k+1}^{i})$
for $i\in[N]$ (see (17))

Set
$\hat{\mathcal{Z}}_{k}\leftarrow\hat{\mathcal{Z}}_{k+1}\times\frac{1}{N}\sum_{i\in[N]}\omega_{k}^{i}$

Normalize
$\omega_{k}^{i}\leftarrow\omega_{k}^{i}/\sum_{j\in[N]}\omega_{k}^{j}$

Resample and MCMC.

if $(\sum_{i\in[N]}(\omega_{k}^{i})^{2})^{-1}<\alpha N$ then

$X_{k}^{1:N}\leftarrow\operatorname{resample}(\tilde{X}_{k}^{1:N},\omega_{k}^{1:N})$
(see Section 3.3)

(Optional). Sample $X_{k}^{i}\leftarrow\mathfrak{M}_{k}(X_{k}^{i},\cdot)$ for
$i\in[N]$ using a $\hat{\pi}_{k}$-invariant MCMC kernel.

Reset $\omega_{k}^{1:N}\leftarrow 1/N$

else

$X_{k}^{1:N}\leftarrow\tilde{X}_{k}^{1:N}$

end if

end for

Particle estimates $\hat{\pi}^{N}=\frac{1}{N}\sum_{i\in[N]}\delta_{X_{0}^{i}}$
of $\pi$ and $\hat{\mathcal{Z}}^{N}_{0}$ of $\mathcal{Z}$

##  Appendix B Asymptotic error formulae for SMC

The results of Webber (2019) play a fundamental role in our work. Compared to
traditional SMC literature, they put more stress on the normalizing constant
estimate, have results for different resampling schemes, and do not require
weight boundedness.

In this section, we present these results using the language of chi-squared
divergence. This gives alternative expressions which are easier to manipulate
in our context.

###  B.1 Chi-squared divergence

For two probability distributions $p$ and $q$ such that $q\ll p$, the chi-
squared divergence of $q$ with respect to $p$ is defined as

$\chi^{2}(q||p)=\int_{\mathcal{X}}p(\differential
x)\left[\left(\frac{\differential q}{\differential p}\right)^{2}(x)-1\right].$

We state without proof two simple properties of this distance.

######  Lemma B.1.

Let $g$ be a nonnegative function from $\mathcal{X}$ to $\mathbb{R}$ such that
$p(g):=\mathbb{E}_{p}[g(X)]<\infty$. Define the probability distribution $\pi$
by $\pi(\differential x)\propto p(x)g(x)$. Then

$\operatorname{Var}_{p}(g)=\left(p(g)\right)^{2}\chi^{2}(\pi||p).$

######  Lemma B.2.

The chi-squared divergence between two probability distributions $p$ and $q$
on $\mathcal{X}\times\mathcal{Y}$ satisfies the decomposition

$\chi^{2}(q(\differential x,\differential y)||p(\differential x,\differential
y))=\chi^{2}(q(\differential x)||p(\differential x))+\int p(\differential
x)\left(\frac{\differential q}{\differential
q}\right)^{2}(x)\chi^{2}(q(\differential y|x)||p(\differential y|x)).$

###  B.2 Generic Feynman-Kac formula and the associated SMC algorithm

Let

$\mathbb{M}(\differential x_{0:T})=\mathbb{M}(\differential
x_{0})\mathbb{M}(\differential x_{1}|x_{0})\ldots\mathbb{M}(\differential
x_{T}|x_{T-1})$

be a Markov measure defined on
$\mathcal{X}_{0}\times\ldots\times\mathcal{X}_{T}$. Let $G_{0}(x_{0})$,
$G_{1}(x_{0},x_{1})$, $\ldots$, $G_{T}(x_{T-1},x_{T})$ be strictly positive
functions. For any $t<T$, assume that there exists $Z_{t}>0$ such that

$\mathbb{Q}_{t}(\differential
x_{0:T})=\frac{1}{Z_{t}}\mathbb{M}_{0}(\differential
x_{0})G_{0}(x_{0})\mathbb{M}_{1}(\differential
x_{1}|x_{0})G_{1}(x_{0},x_{1})\ldots\mathbb{M}_{t}(\differential
x_{t}|x_{t-1})G_{t}(x_{t})\mathbb{M}(\differential x_{t+1:T}|x_{t})$

is a probability measure. Then $\mathbb{Q}_{t}(\differential x_{0:t})$ is
called the Feynman-Kac model associated with the Markov measure
$\mathbb{M}(\differential x_{0:t})$ and the weight functions $(G_{s})_{s\leq
t}$. Given a number of particles $N$, Algorithm 4 approximates
$\mathbb{Q}_{T}(\differential x_{T})$ and the normalizing constant $Z_{T}$. It
is called the SMC algorithm associated to the Feynman-Kac model
$\mathbb{Q}_{T}$.

0: Markov kernels $\mathbb{M}(\differential x_{t}|x_{t-1})$; Functions
$G_{t}(x_{t-1},x_{t})$; Number of particles $N$

Sample $X_{0}^{1:N}\overset{\mathrm{iid}}{\sim}\mathbb{M}_{0}(\differential
x_{0})$

Set $\omega_{0}^{n}=G_{0}(X_{0}^{n})$ and
$Z_{0}^{N}=\frac{1}{N}\sum\omega_{0}^{n}$

Normalize $\omega_{0}^{n}\leftarrow\omega_{0}^{n}/\sum_{m}\omega_{0}^{m}$

for $t=1,\ldots,T$ do

Resample particles $\tilde{X}_{t-1}^{1:N}$ among particles $X_{t-1}^{1:N}$
with weights $\omega_{t-1}^{1:N}$

Sample $X_{t}^{n}\sim\mathbb{M}_{t}(\differential x_{t}|\tilde{X}_{t-1}^{n})$

Set $\omega_{t}^{n}=G_{t}(\tilde{X}_{t-1}^{n},X_{t}^{n})$ and
$Z_{t}^{n}=Z_{t-1}^{n}\frac{1}{N}\sum\omega_{t}^{n}$

Normalize $\omega_{t}^{n}\leftarrow\omega_{t}^{n}/\sum\omega_{t}^{m}$

end for

Empirical measure $\sum\omega_{T}^{n}\delta_{X_{T}^{n}}$ approximating
$\mathbb{Q}_{T}(\differential x_{T})$ and estimate $Z_{T}^{N}$ approximating
$Z_{T}$

Algorithm 4 Generic SMC algorithm

###  B.3 Asymptotic error

We recall the following definition from Webber (2019).

######  Definition B.3.

The notation $|Y_{n}-c|\lesssim U_{n}$ means that there exists a sequence of
sets $(B_{n})$ such that $\mathbb{P}(B_{n})\to 1$ and

$\limsup_{n\to\infty}\mathbb{E}\left[\mathbbm{1}_{B_{n}}\absolutevalue{\frac{Y_{n}-c}{U_{n}}}^{2}\right]\leq
1.$

We are now ready to restate parts of Theorem 3.2 and Example 3.4 of Webber
(2019) in terms of chi-squared divergences.

######  Theorem B.4.

Given the Feynman-Kac model defined in Section B.2, assume that
$\chi^{2}(\mathbb{Q}_{T}(\differential x_{0})||\mathbb{M}_{0}(\differential
x_{0}))<\infty$ and $\chi^{2}(\mathbb{Q}_{T}(\differential
x_{t})||\mathbb{Q}_{t}(\differential x_{t}))<\infty$, $\forall t$. Then
$\sqrt{N}(Z_{T}^{N}/Z_{T}-1)$ is asymptotically normal with variance
$\sigma^{2}_{\mathrm{mult}}$ if multinomial resampling is used;
$\absolutevalue{\sqrt{N}(Z_{T}^{N}/Z_{T}-1)}^{2}\lesssim\sigma^{2}_{\mathrm{sort}}$
if sorted resampling (Gerber et al., 2019) is used; with

$\begin{split}\sigma^{2}_{\mathrm{mult}}&=\chi^{2}(\mathbb{Q}_{T}(\differential
x_{0})||\mathbb{M}_{0}(\differential
x_{0}))+\sum_{t=1}^{T}\chi^{2}(\mathbb{Q}_{T}(\differential
x_{t-1},\differential x_{t})||\mathbb{Q}_{t-1}(\differential
x_{t-1})M_{t}(x_{t-1},\differential x_{t}))\\\
&=\underbrace{\chi^{2}(\mathbb{Q}_{T}(\differential
x_{0})||\mathbb{M}_{0}(\differential
x_{0}))+\sum_{t=1}^{T}\int\mathbb{Q}_{t-1}(\differential
x_{t-1})\left(\frac{\mathbb{Q}_{T}(\differential
x_{t-1})}{\mathbb{Q}_{t-1}(\differential
x_{t-1})}\right)^{2}\chi^{2}(\mathbb{Q}_{T}(\differential
x_{t}|x_{t-1})||\mathbb{M}_{t}(\differential
x_{t}|x_{t-1}))}_{\sigma^{2}_{\mathrm{sort}}}+\\\
&+\sum_{t=1}^{T}\chi^{2}(\mathbb{Q}_{T}(\differential
x_{t-1})||\mathbb{Q}_{t-1}(\differential x_{t-1}))\end{split}$

where we have decomposed $\sigma^{2}_{\mathrm{mult}}$ using the chain rule
(Lemma B.2).

###### Proof.

The original formulation of Theorem 3.2 (Webber, 2019) is written in terms of
the following quantities

$\displaystyle\tilde{G}_{t}$
$\displaystyle:=\mathbb{E}_{\mathbb{M}}\left[\prod_{s=0}^{t-1}G_{s}\right]G_{t}/\mathbb{E}\left[\prod_{s=0}^{t}G_{s}\right]$
$\displaystyle h_{t}(x_{t})$
$\displaystyle:=\mathbb{E}\left[\left.{\prod_{s=t+1}^{T}G_{s}}\right|{X_{t}=x_{t}}\right].$

To translate these notations into our case, note that
$\tilde{G}_{t}=\mathbb{Q}_{t}(\differential
x_{0:T})/\mathbb{Q}_{t-1}(\differential x_{0:T})$ and thus

$\begin{split}\min_{c\in\mathbb{R}}\mathbb{E}_{\mathbb{M}}\left[\prod_{s=0}^{t}\tilde{G}_{s}|h_{t}-c|^{2}\right]&=\min_{c\in\mathbb{R}}\mathbb{E}_{\mathbb{Q}_{t}}\left[|h_{t}-c|^{2}\right]=\operatorname{Var}_{\mathbb{Q}_{t}}(h_{t})=\left[\mathbb{Q}_{t}(h_{t})\right]^{2}\chi^{2}(\mathbb{Q}_{T}(\differential
x_{t})||\mathbb{Q}_{t}(\differential x_{t}))\end{split}$

using Lemma B.1. Moreover,

$\operatorname{Var}_{\mathbb{M}}(G_{t+1}(X_{t},X_{t+1})h_{t+1}(X_{t+1})|X_{t})=h_{t}^{2}(X_{t})\chi^{2}(\mathbb{Q}_{T}(\differential
x_{t+1}|x_{t})||\mathbb{M}_{T}(\differential x_{t+1}|x_{t}))$

and thus, using $h_{t}(x_{t})=\frac{\mathbb{Q}_{T}(\differential
x_{t})}{\mathbb{Q}_{t}(\differential x_{t})}\mathbb{Q}_{t}(h_{t})$ we get

$\mathbb{E}\left[\prod_{s=0}^{t}G_{s}\operatorname{Var}(G_{t+1}h_{t+1}|X_{t})\right]=Z_{t}\mathbb{Q}_{t}(h_{t})^{2}\mathbb{E}_{\mathbb{Q}_{t}}\left[\left\\{\frac{\mathbb{Q}_{T}(\differential
x_{t})}{\mathbb{Q}_{t}(\differential
x_{t})}(X_{t})\right\\}^{2}\chi^{2}(\mathbb{Q}_{T}(\differential
x_{t+1}|X_{t})||\mathbb{M}_{T}(\differential x_{t+1}|X_{t})\right].$

The identity $Z_{t}\mathbb{Q}_{t}(h_{t})=Z_{T}$ helps conclude the proof. ∎

##  Appendix C Proofs

###  C.1 Proof of Proposition 2.1

###### Proof.

Write

$\begin{split}\pi_{t}(x_{t})&=\int\pi_{0}(x_{0})p(x_{t}|x_{0})\differential
x_{0}=\frac{1}{\mathcal{Z}}\int
g_{0}(x_{0})p_{0}(x_{0})p(x_{t}|x_{0})\differential x_{0}\\\
&=\frac{1}{\mathcal{Z}}\int
g_{0}(x_{0})p_{0}(x_{t})p(x_{0}|x_{t})\differential x_{0}\\\
&=\frac{1}{\mathcal{Z}}p_{0}(x_{t})g_{t}(x_{t}).\end{split}$ (21)

The proof is concluded by taking the log gradient of the obtained identity
with respect to $x_{t}$. ∎

###  C.2 Proof of Lemma 2.2

We first state the following elementary lemma on the solution of a linear SDE
(Kloeden & Platen, 1992, Chapter 4.2).

######  Lemma C.1.

Let $a:[0,T]\to\mathbb{R}$ and $c:[0,T]\to\mathbb{R}$ be two continuous
functions. Put $\Phi_{t}=\exp\int_{0}^{t}a_{s}\differential s$, $\forall
t\in[0,T]$. Then the solution to

$\differential Z_{t}=(a_{t}Z_{t}+c_{t})\differential t+\sqrt{2}\differential
W_{t}$ (22)

is

$Z_{t}=\Phi_{t}\left(Z_{0}+\int_{0}^{t}\Phi_{s}^{-1}c_{s}\differential
s+\int_{0}^{t}\Phi_{s}^{-1}\sqrt{2}\differential W_{s}\right).$ (23)

Using Ito isometry, we get the following straightforward corollary.

######  Corollary C.2.

Under the setting of Lemma C.1, if $Z_{0}\sim\mathcal{N}(0,1)$, then

$\displaystyle\mathbb{E}[Z_{t}]=\Phi_{t}\int_{0}^{t}\Phi_{s}^{-1}c_{s}\differential
s,\quad\operatorname{Var}(Z_{t})=\Phi_{t}^{2}+2\Phi_{t}^{2}\int_{0}^{t}\Phi_{s}^{-2}\differential
s.$ (24)

We are now ready to give the proof of Proposition 2.2.

###### Proof of Proposition 2.2.

Without loss of generality, we can assume that $\beta_{t}\equiv 1$. Indeed,
putting $\hat{\bar{g}}_{t}(x_{t}):=g_{0}(e^{-t}x_{t})$ and defining
${\bar{Z}_{t}}^{(T)}$ by

$\differential{\bar{Z}_{t}}^{(T)}=\left[-{\bar{Z}_{t}}^{(T)}+2\nabla\log\hat{\bar{g}}_{T-t}({\bar{Z}_{t}}^{(T)})\right]\differential
t+\sqrt{2}\differential\tilde{B}_{t},\quad{\bar{Z}_{0}}^{(T)}\sim\mathcal{N}(0,1),$
(25)

we see that $Z_{0:t}^{(T)}$ has the same law as
$\bar{Z}^{\int_{0}^{T}\beta_{s}\differential
s}_{\int_{T-t}^{T}\beta_{s}\differential s}$.

We only consider the case $\sigma\neq 1$ here. The case $\sigma=1$ can be
treated using similar but simpler calculations; it can also be recovered by
letting $\sigma\to 1$ in the expressions below.

Since $\nabla\log g_{0}(x_{0})=\nabla\log\pi(x_{0})-\nabla\log
p_{0}(x_{0})=-(x-\mu)/\sigma^{2}+x$, Equation (9) has the form (22) with

$a_{t}=-\left[1+2e^{2(t-T)}(1/\sigma^{2}-1)\right],\quad
c_{t}=2e^{t-T}\mu/\sigma^{2}.$ (26)

Tedious but standard calculations yield

$\Phi_{t}=\exp{-t-e^{-2T}(1/\sigma^{2}-1)(e^{2t}-1)}$ (27)

and the integral of $t$-dependent terms in
$\int_{0}^{T}\Phi_{t}^{-1}c_{t}\differential t$ is

$\int_{0}^{T}\exp{2t+(1/\sigma^{2}-1)e^{2(t-T)}}\differential
t=\frac{1}{2e^{-2T}(1/\sigma^{2}-1)}\left[\exp{1/\sigma^{2}-1}-\exp{e^{-2T}(1/\sigma^{2}-1)}\right].$

Using Corollary C.2, we have

$\mathbb{E}[Z_{T}^{(T)}]=\frac{\mu}{1-\sigma^{2}}\left[1-\exp{(e^{-2T}-1)(1/\sigma^{2}-1)}\right]$
(28)

and $\operatorname{Var}(Z_{T}^{(T)})=v_{T,1}+v_{T,2}$, where

$\displaystyle v_{T,1}$
$\displaystyle=\frac{1}{2(1/\sigma^{2}-1)}\left[1-\exp{-2(1/\sigma^{2}-1)(1-e^{-2T})}\right]$
(29) $\displaystyle v_{T,2}$
$\displaystyle=\exp{-2T-2(1/\sigma^{2}-1)(1-e^{-2T})}.$ (30)

The lemma is proved. ∎

###  C.3 Proof of Propositions 3.1 and 3.2

The propositions follow from a direct application of Theorem B.4 to the
Feynman-Kac model

$\mathbb{Q}_{K}(y_{0:K})=\frac{1}{\mathcal{Z}}\mathbb{M}_{K}(y_{0:K})G_{0}(y_{0})\prod_{k=1}^{K}G_{k}(y_{k-1},y_{k})$

where we have the correspondence $y_{k}=x_{K-k}$ and

$\displaystyle\mathbb{M}_{K}(y_{0:K})$
$\displaystyle=\mathcal{N}(x_{K}|0,\operatorname{Id})\hat{\pi}(x_{K-1}|x_{K})\ldots\hat{\pi}(x_{0}|x_{1})$
$\displaystyle G_{0}(y_{0})$ $\displaystyle=\hat{g}_{K}(x_{K})$ $\displaystyle
G_{k}(y_{k-1},y_{k})$ $\displaystyle=\omega_{K-k}(x_{K-k},x_{K-k+1}).$

###  C.4 Proof of Proposition 3.3

We start by providing some intuition for the result. Repeat that $X_{k}$ is a
shorthand for $X_{k\delta}$ where $\delta$ is the discretization step size.
This convention only applies for $X_{k}$.

When $K$ is big ($\delta$ is small), we have, as established in (14),

$\displaystyle\pi(x_{k}|x_{k+1})$
$\displaystyle\approx\mathcal{N}(x_{k};\sqrt{1-\alpha_{k+1}}x_{k+1}+\alpha_{k+1}\nabla\log
g_{k+1}(x_{k+1}),\alpha_{k+1}I)$ $\displaystyle\hat{\pi}(x_{k}|x_{k+1})$
$\displaystyle=\mathcal{N}(x_{k};\sqrt{1-\alpha_{k+1}}x_{k+1}+\alpha_{k+1}\nabla\log\hat{g}_{k+1}(x_{k+1}),\alpha_{k+1}I).$

Using the analytic formula for the chi-squared divergence between two
Gaussians we get

$\displaystyle\chi^{2}(\pi(x_{k}|x_{k+1})||\hat{\pi}(x_{k}|x_{k+1}))\approx
e^{\alpha_{k+1}\norm{\nabla\log
g_{k+1}-\nabla\log\hat{g}_{k+1}}^{2}(x_{k+1})}-1\approx
2\delta\norm{\nabla\log g_{k+1}-\nabla\log\hat{g}_{k+1}}^{2}(x_{k+1}).$

Thus

$\displaystyle\zeta_{K}^{2}$
$\displaystyle=\chi^{2}(\pi_{K}||\mathcal{N}(0,\operatorname{Id}))+\sum_{k}\int\frac{\pi_{k+1}(x_{k+1})^{2}}{\hat{\pi}_{k+1}(x_{k+1})}\chi^{2}(\pi(x_{k}|x_{k+1})||\hat{\pi}(x_{k}|x_{k+1}))\differential
x_{k+1}$
$\displaystyle\approx\chi^{2}(\pi_{K}||\mathcal{N}(0,\operatorname{Id}))+\sum_{k}\int\frac{\pi_{k+1}(x_{k+1})^{2}}{\hat{\pi}_{k+1}(x_{k+1})}2\delta\norm{\nabla\log
g_{k+1}-\nabla\log\hat{g}_{k+1}}^{2}(x_{k+1})\differential x_{k+1}$
$\displaystyle\approx\chi^{2}(\pi_{T}||\mathcal{N}(0,\operatorname{Id}))+2\int_{0}^{T}\int_{\mathcal{X}}\frac{\pi_{t}(x)^{2}}{\hat{\pi}_{t}(x)}\norm{\nabla\log
g_{t}(x)-\nabla\log\hat{g}_{t}(x)}^{2}\differential x\differential t.$

####  C.4.1 Regularity conditions for Proposition 3.3

We assume that the sequence of distributions $\pi_{t}(\cdot)$ satisfy the
following properties.

######  Assumption C.3.

There exists $M_{1}>0$ such that $\norm{\nabla\log\pi_{t}(x_{t})}\leq
M_{1}(1+\norm{x_{t}})$.

######  Assumption C.4.

There exist $M_{2}>0$, $M_{3}>0$, $\alpha_{2}\geq 1$, and $\alpha_{3}\geq 1$
such that $\norm{\nabla^{2}\log\pi_{t}(x_{t})}\leq
M_{2}(1+\norm{x_{t}})^{\alpha_{2}}$ and
$\norm{\nabla^{3}\log\pi_{t}(x_{t})}\leq M_{3}(1+\norm{x_{t}})^{\alpha_{3}}$.

######  Assumption C.5.

There exist $\vartheta>0$ and $M_{\infty}<\infty$ such that
$\int\pi_{t}(x_{t})e^{\vartheta\norm{x_{t}}^{2}}\differential
x_{t}<M_{\infty}$, $\forall t$.

These assumptions are satisfied, for example, when the target distribution
$\pi_{0}$ is Gaussian.

######  Remark C.6.

Let $\varphi(x_{t})=\nabla\log\pi_{t}(x_{t})$. Then the notation
$\nabla^{3}\log\pi_{t}(x_{t})$ refers to the second order differential
$\varphi^{\prime\prime}(x_{t})$ which is a bilinear mapping from
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ to $\mathbb{R}^{d}$. For any multilinear
operator $H:\mathcal{X}_{1}\times\ldots\times\mathcal{X}_{n}\to\mathcal{Y}$,
we define

$\norm{H}_{\operatorname{op}}:=\sup_{x_{1},\ldots,x_{n}\neq
0}\frac{\norm{H(x_{1},\ldots,x_{n})}}{\norm{x_{1}}\ldots\norm{x_{n}}}.$

By writing $\norm{H}$ we implicitly refer to $\norm{H}_{\operatorname{op}}$.
In fact, the space of such operators is of finite dimensions in our cases of
interests, hence any two norms are bounded by a constant factor of each other.

The above assumptions only concern the differential of
$\nabla\log\pi_{t}(x_{t})$ with respect to $x$. The following lemma derives a
bound with respect to $t$.

######  Lemma C.7.

Under the above assumptions, there exist constants $\bar{M}_{1}$ and
$\bar{\alpha}_{1}$ such that

$\norm{\frac{\partial}{\partial
t}\nabla\log\pi_{t}(x_{t})}\leq\bar{M}_{1}(1+\norm{x_{t}})^{\bar{\alpha}_{1}},\forall
t.$

###### Proof.

Using the Fokker-Planck equation for the score (Lai et al., 2022), we write

$\partial_{t}\log g_{t}(x_{t})=\operatorname{div}_{x}\nabla\log
g_{t}(x_{t})+\nabla\log g_{t}(x_{t})\circ(\nabla\log g_{t}(x_{t})-x_{t}).$

For a fixed $t$, put $\varphi(x_{t})=\nabla\log g_{t}(x_{t})$ and
$\psi(x_{t})=\Tr(\varphi^{\prime}(x_{t}))+\varphi(x_{t})\circ(\varphi(x_{t})-x_{t})$.
Then $\partial_{t}\nabla\log g_{t}(x_{t})=\nabla\psi(x_{t})$. Viewing
$\psi^{\prime}(x_{t})$, $\varphi^{\prime}(x_{t})$, and
$\varphi^{\prime\prime}(x_{t})$ as elements of
$\mathcal{L}(\mathbb{R}^{d},\mathbb{R})$,
$\mathcal{L}(\mathbb{R}^{d},\mathbb{R}^{d})$, and
$\mathcal{L}(\mathbb{R}^{d},\mathcal{L}(\mathbb{R}^{d},\mathbb{R}^{d}))$
respectively; where $\mathcal{L}(A,B)$ is the space of linear operators from
$A$ to $B$; we can write

$\psi^{\prime}(x_{t})h=\Tr(\varphi^{\prime\prime}(x_{t})h)+(\varphi^{\prime}(x_{t})h)\circ(\varphi(x_{t})-x_{t})+\varphi(x_{t})\circ(\varphi^{\prime}(x_{t})h-h),\forall
h\in\mathbb{R}^{d}$

where $\circ$ stands for the usual scalar product between two vectors in
$\mathbb{R}^{d}$. There is a constant $C$ depending on the dimension such that
$\Tr(L)\leq C\norm{L}_{\operatorname{op}}$ for endomorphisms $L$. Thus

$\displaystyle\norm{\psi^{\prime}(x_{t})h}_{\operatorname{op}}$
$\displaystyle\leq
C\norm{\varphi^{\prime\prime}(x_{t})}_{\operatorname{op}}\norm{h}+\norm{\varphi^{\prime}(x_{t})}_{\operatorname{op}}\norm{\varphi(x_{t})-x_{t}}\norm{h}+\norm{\varphi(x_{t})}\norm{\varphi^{\prime}(x_{t})-\operatorname{Id}}_{\operatorname{op}}\norm{h}$
$\displaystyle\leq\bar{M}_{1}(1+\norm{x_{t}})^{\bar{\alpha}_{1}}\norm{h}$

for some $\bar{M}_{1}$ and $\bar{\alpha}_{1}$ by Assumptions C.3 and C.4. This
entails $\norm{\partial_{t}\nabla\log
g_{t}(x_{t})}=\norm{\psi^{\prime}(x_{t})}_{\operatorname{op}}\leq\bar{M}_{1}(1+\norm{x_{t}})^{\bar{\alpha}_{1}}$.
∎

####  C.4.2 Formal proof

To formalize the error of the heuristic approximations presented at the
beginning of Section C.4 we need the following technical lemma.

######  Lemma C.8.

Let $x_{0}$ and $v$ be two vectors in $\mathbb{R}^{d}$ and suppose that
$X^{\mathrm{A}}_{t}$ and $X^{\mathrm{B}}_{t}$ are respectively the unique
solutions of the SDEs:

$\displaystyle\mathrm{(A):}\quad\differential
X^{\mathrm{A}}_{t}=(-X^{\mathrm{A}}_{t}+2v)\differential
t+\sqrt{2}\differential W^{\mathrm{A}}_{t},\quad X_{0}=x_{0}$
$\displaystyle\mathrm{(B):}\quad\differential
X^{\mathrm{B}}_{t}=(-X^{\mathrm{B}}_{t}+2f_{t}(X^{\mathrm{B}}_{t}))\differential
t+\sqrt{2}\differential W^{\mathrm{B}}_{t},\quad X_{0}=x_{0}$

where there exist strictly positive constants $M_{1}$, $\bar{M}_{1}$, $M_{2}$,
and $M_{3}$; and strictly greater than $1$ constants $\bar{\alpha}_{1}$,
$\alpha_{2}$, and $\alpha_{3}$; such that the function $f_{t}(x_{t})$
satisfies $\norm{f_{t}(x_{t})}\leq M_{1}(1+\norm{x_{t}})$,
$\norm{\partial_{t}f_{t}(x_{t})}\leq\bar{M}_{1}(1+\norm{x_{t}})^{\bar{\alpha}_{1}}$,
and $\norm{\nabla^{i}f_{t}(x_{t})}\leq M_{i+1}(1+\norm{x_{t}})^{\alpha_{i+1}}$
for $i\in\left\\{1,2\right\\}$. (The notation $\nabla$ refers implicitly to
the gradient with respect to $x$.) Denote
$\mathbb{P}_{\mathrm{A}}(\differential x_{[0,T]})$ and
$\mathbb{P}_{\mathrm{B}}(\differential x_{[0:T]})$ respectively the path
measures associated with the solutions of (A) and (B). Then, there exist a
parameter $\widetilde{M}$ depending on all the aforementioned constants and a
parameter $\widetilde{M}_{1}$ depending only on $M_{1}$ such that for any
$t\leq 1/\widetilde{M}_{1}$, the chi-squared divergence of
$\mathbb{P}_{\mathrm{B}}(\differential x_{[0,t]})$ with respect to
$\mathbb{P}_{\mathrm{A}}(\differential x_{[0,t]})$ is finite, and

$\absolutevalue{\chi^{2}(\mathbb{P}_{\mathrm{B}}(\differential
x_{[0:t]})||\mathbb{P}_{\mathrm{A}}(\differential
x_{[0:t]}))-2t\norm{f_{0}(x_{0})-v}^{2}}\leq\widetilde{M}t^{2}e^{t\widetilde{M}_{1}(\norm{x_{0}}^{2}+\norm{v}^{2})}\left(1+\norm{x_{0}}+\norm{v}\right)^{4(1+\bar{\alpha}_{1}+\alpha_{2}+\alpha_{3})}$

###### Proof.

Put $\Delta_{t}(x_{t})=2f_{t}(x_{t})-2v$. By an application of Girsanov’s
theorem (Liptser & Shiryayev, 1977, Example 3, Section 6.2.3), we have

$D_{t}:=\frac{\differential\mathbb{P}_{\mathrm{B}}}{\differential\mathbb{P}_{\mathrm{A}}}(X_{[0,t]})=\exp{\int_{0}^{t}\frac{\langle\Delta_{s}(X_{s}),\differential
W^{\mathrm{A}}_{s}\rangle}{\sqrt{2}}-\frac{1}{4}\int_{0}^{t}\norm{\Delta_{s}(X_{s})}^{2}\differential
s}$

where, for a vector-valued process $V_{t}$, the notation $\int_{0}^{t}\langle
V_{s},\differential
W_{s}\rangle:=\sum_{i=1}^{d}\int_{0}^{t}V_{s}^{i}\differential W_{s}^{i}$. As
a preliminary step, we would like to bound
$\mathbb{E}_{\operatorname{A}}[D_{t}^{\alpha}(1+\norm{X_{t}}+\norm{v})^{n}]$
for some $\alpha,n\geq 1$. Here, $c(\cdot)$ denotes a constant whose value
might change from line to line and depends on the variables inside the round
bracket. We also drop the subscript/superscript A from
$\mathbb{E}_{\operatorname{A}}$ and $W^{\mathrm{A}}_{t}$ whenever there is no
risk of confusion. We have

$\begin{split}&\mathbb{E}[D_{t}^{\alpha}(1+\norm{X_{t}}+\norm{v})^{n}]=\mathbb{E}\left[\exp{\alpha\int_{0}^{t}\frac{\langle\Delta_{s},\differential
W_{s}\rangle}{\sqrt{2}}-\frac{\alpha}{4}\int_{0}^{t}\norm{\Delta_{s}}^{2}\differential
s}(1+\norm{X_{t}}+\norm{v})^{n}\right]\\\
&=\mathbb{E}\left[\exp{\alpha\int_{0}^{t}\frac{\langle\Delta_{s},\differential
W_{s}\rangle}{\sqrt{2}}-\frac{\alpha^{2}}{2}\int_{0}^{t}\norm{\Delta_{s}}^{2}\differential
s}\exp{\left(\frac{\alpha^{2}}{2}-\frac{\alpha}{4}\right)\int_{0}^{t}\norm{\Delta_{s}}^{2}\differential
s}(1+\norm{X_{t}}+\norm{v})^{n}\right]\\\ &\leq
e^{t(c(M_{1},\alpha)+c(\alpha)\norm{v}^{2})}\mathbb{E}\bigg{[}\exp{\alpha\int_{0}^{t}\frac{\langle\Delta_{s},\differential
W_{s}\rangle}{\sqrt{2}}-\frac{\alpha^{2}}{2}\int_{0}^{t}\norm{\Delta_{s}}^{2}\differential
s}\exp{c(M_{1},\alpha)\int_{0}^{t}\norm{X_{s}}^{2}\differential s}\times\\\
&\times(1+\norm{X_{t}}+\norm{v})^{n}\bigg{]}\text{ using
}\norm{\Delta_{s}}^{2}\leq c(M_{1})(1+\norm{x_{s}}^{2})+8\norm{v}^{2}\\\ &\leq
e^{t(c(M_{1},\alpha)+c(\alpha)\norm{v}^{2})}\mathbb{E}^{1/2}\left[\exp{\sqrt{2}\alpha\int_{0}^{t}\langle\Delta_{s},\differential
W_{s}\rangle-\alpha^{2}\int_{0}^{t}\norm{\Delta_{s}}^{2}\differential
s}\right]\times\\\
&\times\mathbb{E}^{1/4}\left[\exp{4c(M_{1},\alpha)\int_{0}^{t}\norm{X_{s}}^{2}\differential
s}\right]\mathbb{E}^{1/4}\left[(1+\norm{X_{t}}+\norm{v})^{n}\right]\end{split}$

using double Cauchy-Schwarz
$\mathbb{E}[XYZ]\leq\mathbb{E}^{1/2}[X^{2}]\mathbb{E}^{1/4}[Y^{4}]\mathbb{E}^{1/4}[Z^{4}]$.

In the last line of the above display, the first expectation is equal to $1$
by the same Girsanov argument as Liptser & Shiryayev (1977, Example 3, Section
6.2.3). To bound the second expectation, we note that under
$\mathbb{P}_{\operatorname{A}}$, we have
$X_{s}\sim\mathcal{N}(\sqrt{1-\lambda_{s}}(x_{0}-2v)+2v,\lambda_{s})$.
Elementary calculations yield the bound

$\mathbb{E}[e^{k\norm{X_{s}}^{2}}]=\left(\frac{1}{\sqrt{1-2k\lambda_{s}}}\right)^{d}\exp{\frac{k\norm{\sqrt{1-\lambda_{s}}(x_{0}-2v)+2v}^{2}}{1-2k\lambda_{s}}}\leq(\sqrt{2})^{d}e^{16k(\norm{x_{0}}^{2}+\norm{v}^{2})}$

for $0<k<1/4$. Write

$\begin{split}\mathbb{E}\left[\exp{4c(M_{1},\alpha)\int_{0}^{t}\norm{X_{s}}^{2}\differential
s}\right]\leq\frac{1}{t}\int_{0}^{t}\mathbb{E}\left[e^{4tc(M_{1},\alpha)\norm{X_{s}}^{2}}\right]\differential
s\leq c(1)e^{64tc(M_{1},\alpha)(\norm{x_{0}}^{2}+\norm{v}^{2})}\end{split}$

if $t\leq\frac{1}{16c(M_{1},\alpha)}$, using Jensen’s inequality and the above
bound. The third expectation can be bounded by
$c(n)(1+\norm{x_{0}}^{4n}+\norm{v_{0}}^{4n})$.

Putting everything together, we establish that there exist constants $c(n)$
and $c(M_{1},\alpha)$ such that

$\mathbb{E}[D_{t}^{\alpha}(1+\norm{X_{t}}+\norm{v})^{n}]\leq
c(n)e^{tc(M_{1},\alpha)(\norm{x_{0}}^{2}+\norm{v}^{2})}(1+\norm{x_{0}}^{4n}+\norm{v}^{4n}),\forall
t\leq\frac{1}{c(M_{1},\alpha)}.$ (31)

Now to study $\chi^{2}(\mathbb{P}_{\mathrm{B}}(\differential
x_{[0:t]})||\mathbb{P}_{\mathrm{A}}(\differential x_{[0:t]}))$, we apply Ito’s
formula to $D_{t}^{2}$ and get, under $\mathbb{P}_{\operatorname{A}}$

$D_{t}^{2}=1+\sqrt{2}\int_{0}^{t}D_{s}^{2}\langle\Delta_{s},\differential
W_{s}\rangle+\int_{0}^{t}\frac{D_{s}^{2}\norm{\Delta_{s}}^{2}}{2}\differential
t.$ (32)

Putting $\eta_{t}(x_{t})=\norm{\Delta_{t}(x_{t})}^{2}$ and
$\tilde{f}_{t}(x_{t})=2f_{t}(x_{t})-x_{t}$, we have

$D_{t}^{2}\eta_{t}=\eta_{0}+\sqrt{2}\int_{0}^{t}D_{s}^{2}\langle\eta_{s}\Delta_{s}+\nabla\eta_{s},\differential
W_{s}\rangle+\\\
+\int_{0}^{t}D_{s}^{2}\left\\{\eta_{s}\frac{\norm{\Delta_{s}}^{2}}{2}+\nabla\eta_{s}(\tilde{f}(X_{s})+2\Delta_{s})+\frac{\partial\eta}{\partial
s}+\Tr(\nabla^{2}\eta_{s})\right\\}\differential s.$ (33)

To study (32) and (33), we use (31) together with the following bounds which
are consequences of the lemma’s assumptions:

$\displaystyle\norm{\Delta_{s}(x_{s})}$ $\displaystyle\leq
c(M_{1})(1+\norm{x_{s}}+\norm{v})$ $\displaystyle\eta_{s}(x_{s})$
$\displaystyle\leq c(M_{1})(1+\norm{x_{s}}+\norm{v})^{2}$
$\displaystyle\norm{\nabla\eta_{s}(x_{s})}$ $\displaystyle\leq
c(M_{1},M_{2})(1+\norm{x_{s}}+\norm{v})^{1+\alpha_{2}}$
$\displaystyle\norm{\frac{\partial\eta}{\partial s}(s,X_{s})}$
$\displaystyle\leq c(M_{1},\bar{M}_{1})(1+\norm{x_{s}})^{1+\bar{\alpha}_{1}}$
$\displaystyle\norm{\Tr(\nabla^{2}\eta_{s})}$ $\displaystyle\leq
c(M_{1},M_{2},M_{3})(1+\norm{x_{s}}+\norm{v})^{1+\alpha_{2}+\alpha_{3}}.$

The last inequality is justified by the fact that
$\Tr(\nabla^{2}\eta_{s})\lesssim\norm{\nabla^{2}\eta_{s}}$, where by
considering $\nabla^{2}\eta_{s}(x_{s})$ as the second differential of
$\eta_{s}$ at $x_{s}$ (i.e. a bilinear form from
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ to $\mathbb{R}$), we have

$\frac{\partial^{2}\eta_{s}}{\partial
x^{2}}(x)[h,k]=\Delta_{s}(x)\circ\frac{\partial^{2}\Delta}{\partial
x^{2}}(x)[h,k]+(\frac{\partial\Delta}{\partial
x}(x)h)\circ(\frac{\partial\Delta}{\partial
x}(x)k),\forall(h,k)\in\mathbb{R}^{d}\times\mathbb{R}^{d}$

where $\circ$ stands for the usual scalar product between two vectors in
$\mathbb{R}^{d}$. These bounds show that the stochastic integrals (w.r.t.
$\differential W_{s}$) in (32) and (33) are true martingales (as opposed to
merely local martingales). Moreover, there exist a constant $\widetilde{M}$
depending on $M_{1}$, $M_{2}$, $M_{3}$, $\bar{M}_{1}$, $\bar{\alpha}_{1}$,
$\alpha_{2}$, and $\alpha_{3}$; and a constant $\widetilde{M}_{1}$ depending
on $M_{1}$ only such that

$\mathbb{E}\left[D_{s}^{2}\absolutevalue{\eta_{s}\frac{\norm{\Delta_{s}}^{2}}{2}+\nabla\eta_{s}(\tilde{f}(X_{s})+2\Delta_{s})+\frac{\partial\eta}{\partial
s}+\Tr(\nabla^{2}\eta_{s})}\right]\leq
4\widetilde{M}e^{t\widetilde{M}_{1}(\norm{x_{0}}^{2}+\norm{v}^{2})}\times\\\
\times(1+\norm{x_{0}}+\norm{v})^{4(1+\bar{\alpha}_{1}+\alpha_{2}+\alpha_{3})},\forall
s\leq t\leq\frac{1}{\widetilde{M}_{1}}.$

Taking expectation of both sides of (33) and rearranging yields

$\absolutevalue{\mathbb{E}[D_{t}^{2}\eta_{t}]-\eta_{0}}\leq
4\mathbf{t}\widetilde{M}e^{t\widetilde{M}_{1}(\norm{x_{0}}^{2}+\norm{v}^{2})}(1+\norm{x_{0}}+\norm{v})^{4(1+\bar{\alpha}_{1}+\alpha_{2}+\alpha_{3})},\forall
t\leq\frac{1}{\widetilde{M}_{1}}.$

Then we have, by taking expectation of both sides of (32):

$\begin{split}&\absolutevalue{\chi^{2}(\mathbb{P}_{\mathrm{B}}(\differential
x_{[0:t]})||\mathbb{P}_{\mathrm{A}}(\differential
x_{[0:t]}))-t\frac{\eta_{0}}{2}}=\absolutevalue{\mathbb{E}(D_{t}^{2})-1-t\frac{\eta_{0}}{2}}=\absolutevalue{\int_{0}^{t}\left(\mathbb{E}\left[\frac{D_{s}^{2}\eta_{s}}{2}\right]-\frac{\eta_{0}}{2}\right)\differential
s}\leq\int_{0}^{t}\absolutevalue{\mathbb{E}\left[\frac{D_{s}^{2}\eta_{s}}{2}\right]-\frac{\eta_{0}}{2}}\differential
s\\\
&\leq\int_{0}^{t}2s\widetilde{M}e^{t\widetilde{M}_{1}(\norm{x_{0}}^{2}+\norm{v}^{2})}(1+\norm{x_{0}}+\norm{v})^{4(1+\bar{\alpha}_{1}+\alpha_{2}+\alpha_{3})}\differential
s\\\
&=t^{2}\widetilde{M}e^{t\widetilde{M}_{1}(\norm{x_{0}}^{2}+\norm{v}^{2})}(1+\norm{x_{0}}+\norm{v})^{4(1+\bar{\alpha}_{1}+\alpha_{2}+\alpha_{3})}.\end{split}$

The proof is completed. ∎

We are now ready to give the proof of Proposition 3.3.

###### Proof.

To make the arguments clearer, we shall assume that the proposal distribution
is

$\hat{\pi}(x_{k}|x_{k+1})=\mathcal{N}(x_{k}|\sqrt{1-\alpha_{k+1}}x_{k+1}+2(1-\sqrt{1-\alpha_{k+1}})\nabla\log\hat{g}_{k+1}(x_{k+1}),\alpha_{k+1}I)$

which is slightly different from (16). As we will see, the proof also applies
to the original discretization with minimal changes. For $0<s<u<T$, the
distributions $\hat{\pi}(x_{s}|x_{u})$ and $\pi(x_{s}|x_{u})$ can be obtained
by respectively solving the following SDEs between times $T-u$ and $T-s$:

$\displaystyle\hat{\pi}(x_{s}|x_{u}):\quad\differential
Y^{\mathrm{A}}_{t}=(-Y^{\mathrm{A}}_{t}+2\nabla\log\hat{g}_{T-\mathbf{u}}(Y^{\mathrm{A}}_{t}))\differential
t+\sqrt{2}\differential W^{\mathrm{A}}_{t},\quad Y^{\mathrm{A}}_{T-u}=x_{u}$
$\displaystyle\pi(x_{s}|x_{u}):\quad\differential
Y^{\mathrm{B}}_{t}=(-Y^{\mathrm{B}}_{t}+2\nabla\log
g_{T-t}(Y^{\mathrm{B}}_{t}))\differential t+\sqrt{2}\differential
W^{\mathrm{B}}_{t},\quad Y^{\mathrm{B}}_{T-u}=x_{u}.$

The assumptions in Section C.4.1 and Lemma C.7 show that the conditions of
Lemma C.8 are satisfied for this pair of SDEs. Thus

$\absolutevalue{\chi^{2}(\mathbb{P}_{\mathrm{B}}(\differential
y_{[T-u,T-s]})||\mathbb{P}_{\mathrm{A}}(\differential
y_{[T-u,T-s]}))-2(u-s)\norm{\nabla\log
g_{u}(x_{u})-\nabla\log\hat{g}_{u}(x_{u})}^{2}}\leq\widetilde{M}(u-s)^{2}\times\\\
\times
e^{(u-s)\widetilde{M}_{1}(\norm{x_{u}}^{2}+\norm{\nabla\log\hat{g}_{u}(x_{u})}^{2})}(1+\norm{x_{u}}+\norm{\nabla\log\hat{g}_{u}(x_{u})})^{\alpha_{+}}$

for $\alpha_{+}=4(1+\bar{\alpha}_{1}+\alpha_{2}+\alpha_{3})$ and
$u-s\leq\frac{1}{\widetilde{M}_{1}}$. This, together with the data processing
inequality and the assumption on $|\nabla\log\hat{g}_{t}|$, implies

$\chi^{2}(\pi(x_{k}|x_{k+1})||\hat{\pi}(x_{k}|x_{k+1}))\leq
2\delta\norm{\nabla\log
g_{k+1}-\nabla\log\hat{g}_{k+1}}^{2}(x_{k+1})+\widetilde{M}\delta^{2}e^{\delta\widetilde{M}_{1}(1+2C_{2}^{2})(1+\norm{x_{k+1}})^{2}}\times\\\
\times(1+C_{2})^{\alpha_{+}}(1+\norm{x_{k+1}})^{\alpha_{+}},\forall\delta\leq\frac{1}{\widetilde{M}_{1}}.$

Thus, for a sufficiently fine discretization,

$\sum_{k}\int\frac{\pi_{k+1}(x_{k+1})^{2}}{\hat{\pi}_{k+1}(x_{k+1})}\chi^{2}(\pi(x_{k}|x_{k+1})||\hat{\pi}(x_{k}|x_{k+1}))\differential
x_{k+1}\leq\\\
\leq\sum_{k}\delta\int\frac{\pi_{k+1}(x_{k+1})^{2}}{\hat{\pi}_{k+1}(x_{k+1})}2\norm{\nabla\log
g_{k+1}-\nabla\log\hat{g}_{k+1}}^{2}(x_{k+1})\differential x_{k+1}+\\\
+\sum_{k}\int
C_{1}\pi_{k+1}(x_{k+1})\widetilde{M}\delta^{2}e^{\delta\widetilde{M}_{1}(1+2C_{2}^{2})(1+\norm{x_{k+1}}^{2})}(1+C_{2})^{\alpha_{+}}(1+\norm{x_{k+1}})^{\alpha_{+}}\differential
x_{k+1}.$ (34)

The first term is a Riemann sum and converges to
$\int_{0}^{T}\int\frac{\pi_{t}(x_{t})^{2}}{\hat{\pi}_{t}(x_{t})}2\norm{\nabla\log
g_{t}-\nabla\log\hat{g}_{t}}^{2}(x_{t})\differential x_{t}\differential t$. To
bound the second term, first note that

$\displaystyle\mathbb{E}_{\pi_{k+1}}[(1+\norm{X})^{2\alpha_{+}}]$
$\displaystyle\leq
2^{2\alpha_{+}-1}\mathbb{E}\left[1+\norm{X}^{2\alpha_{+}}\right]\leq
2^{2\alpha_{+}-1}\mathbb{E}\left[1+e^{\vartheta\norm{X}^{2}}\max\left(\frac{\lceil\alpha_{+}\rceil!}{\vartheta^{\lceil\alpha_{+}\rceil}},\frac{\lfloor\alpha_{+}-1\rfloor!}{\vartheta^{\lfloor\alpha_{+}-1\rfloor}}\right)\right]$
$\displaystyle\leq
2^{2\alpha_{+}-1}\left(1+M_{\infty}\max\left(\frac{\lceil\alpha_{+}\rceil!}{\vartheta^{\lceil\alpha_{+}\rceil}},\frac{\lfloor\alpha_{+}-1\rfloor!}{\vartheta^{\lfloor\alpha_{+}-1\rfloor}}\right)\right)$

where $M_{\theta}$ and $\vartheta$ appear in Assumption C.5. Thus, as long as
$\delta\widetilde{M}_{1}(1+2C_{2}^{2})\leq\vartheta/2$, it holds that

$\displaystyle\int\pi_{k+1}(x)e^{\delta\widetilde{M}_{1}(1+2C_{2}^{2})\norm{x}^{2}}(1+\norm{x})^{\alpha_{+}}\differential
x$
$\displaystyle\leq\mathbb{E}_{\pi_{k+1}}\left[e^{\vartheta\norm{X}^{2}/2}(1+\norm{X})^{\alpha_{+}}\right]$
$\displaystyle\leq\mathbb{E}^{1/2}\left[e^{\vartheta\norm{X}^{2}}\right]\mathbb{E}^{1/2}\left[(1+\norm{X})^{2\alpha_{+}}\right]\leq
C(M_{\infty},\alpha_{+},\vartheta)$

for some constant $C(M_{\infty},\alpha_{+})$ depending on $M_{\infty}$,
$\alpha_{+}$, and $\vartheta$. Hence the second sum of (34) tends to $0$ when
$\delta\to 0$. The proof is finished. ∎

###  C.5 Proof of Proposition 4.1

The denoising score matching identity is standard and recalled here for
convenience. We have

$\displaystyle\pi_{k}(x_{k})=\int p(x_{k}|x_{0})\pi_{0}(x_{0})\mathrm{d}x_{0}$

so by using the log derivative we obtain

$\displaystyle\nabla\log\pi_{k}(x_{k})$ $\displaystyle=\int\nabla\log
p(x_{k}|x_{0})\frac{\pi_{0}(x_{0})p(x_{k}|x_{0})}{\pi_{k}(x_{k})}\mathrm{d}x_{0}$
(35) $\displaystyle=\int\nabla\log
p(x_{k}|x_{0})\pi(x_{0}|x_{k})\mathrm{d}x_{0}.$ (36)

It can be easily verified that the interchange of differentiation and
integration here does not require any regularity assumption on
$\pi_{0}(x_{0})$ apart from differentiability.

To prove the novel score identity, we first note that, under the condition
$\int\norm{\nabla\pi(x_{0})}e^{-\eta\norm{x_{0}}^{2}}\differential
x_{0}<\infty,\forall\eta>0$, we have

$\int\nabla_{x_{0}}\log\pi(x_{0}|x_{k})\pi(x_{0}|x_{k})\differential x_{0}=0$
(37)

according to Lemma C.9. Combining this identity with (36), we have, for any
$\alpha\in\mathbb{R}$,

$\nabla\log\pi(x_{k})=\int\left[\nabla_{x_{k}}\log
p(x_{k}|x_{0})+\alpha\nabla_{x_{0}}\log\pi(x_{0}|x_{k})\right]\pi(x_{0}|x_{k})\differential
x_{0}.$

In particular:

  * •

For $\alpha=\frac{1}{\sqrt{1-\lambda_{k}}}$, we get

$\nabla\log\pi(x_{k})=\frac{1}{\sqrt{1-\lambda_{k}}}\int\nabla\log\pi(x_{0})\pi(x_{0}|x_{k})\differential
x_{0}$

which is the identity presented in Appendix C.1.3 (De Bortoli et al., 2021);

  * •

For $\alpha=\sqrt{1-\lambda_{k}}$, we get

$\nabla\log\pi(x_{k})=\int\left(\sqrt{1-\lambda_{k}}\nabla\log
g_{0}(x_{0})-x_{k}\right)\pi(x_{0}|x_{k})\differential x_{0}$

which is the identity we wanted to prove.

The verifications are straightforward by remarking that
$\nabla_{x_{0}}\log\pi(x_{0}|x_{k})=\nabla\log g_{0}(x_{0})+\nabla_{x_{0}}\log
p(x_{0}|x_{k})$. We also note that choosing $\alpha=0$ brings us back to the
classical score matching loss. Therefore, different values of $\alpha$ give
losses with different properties.

We finish this section with a technical lemma giving conditions for (37) to
hold.

######  Lemma C.9.

Let $f:\mathbb{R}^{d}\to\mathbb{R}$ be a probability density, i.e. $f\geq 0$
and $\int f(x)\differential x=1$. Suppose that $f$ is continuously
differentiable and $\int\norm{\nabla f(x)}\differential x<\infty$. Then
$\int\nabla f(x)\differential x=0$.

######  Remark C.10.

The condition $\int\norm{\nabla f(x)}\differential x<\infty$ is clearly the
minimum necessary for $\int\nabla f(x)\differential x=0$ to make sense. On the
other hand, we do not explicitly require that $f$ or $\nabla f$ vanishes at
infinity.

###### Proof.

Without loss of generality, we only prove that
$\int\partial_{1}f(x)\differential x=0$. Put $g(x_{1}):=\int
f(x_{1},x_{2:d})\differential x_{2:d}$. Fubini’s theorem then implies that
$\int_{\mathbb{R}}g(x_{1})\differential x_{1}=1$. We have

$\begin{split}\int_{\mathbb{R}^{d}}\partial_{1}f(x)\differential
x&=\lim_{M\to\infty}\int_{\mathbb{R}^{d-1}}\int_{-M}^{M}\partial_{1}f(x_{1},x_{2:d})\differential
x_{1}\differential
x_{2:d}=\lim_{M\to\infty}\int_{\mathbb{R}^{d-1}}f(M,x_{2:d})-f(-M,x_{2:d})\differential
x_{2:d}\\\ &=\lim_{M\to\infty}g(M)-g(-M).\end{split}$ (38)

Put $I(M):=\int_{0}^{\infty}\absolutevalue{g(M+x)-g(-M-x)}\differential x$. We
have

$I(M)\leq\int_{0}^{\infty}\absolutevalue{g(M+x)}+\absolutevalue{g(-M-x)}\differential
x,$

thus $\lim_{M\to\infty}I(M)=0$ by the integrability of $g$. Combining this
with Fatou’s lemma yields

$0=\liminf_{M\to\infty}I(M)\geq\int_{0}^{\infty}\liminf_{M\to\infty}\absolutevalue{g(M+x)-g(-M-x)}\differential
x=\int_{0}^{\infty}\absolutevalue{\int_{\mathbb{R}^{d}}\partial_{1}f(y)\differential
y}\differential x$

by (38). This means that $\int_{\mathbb{R}^{d}}\partial_{1}f(y)\differential
y=0$. ∎

###  C.6 Proof of Proposition 4.2

###### Proof.

Since we are in the Gaussian case with $d=1$, we have $\nabla\log
g_{0}(x_{0})=ax_{0}+b$ for some $a,b\in\mathbb{R}$. Therefore
$\ell_{\textup{NSM}}(\theta)$ and
$\operatorname{Var}(\hat{\ell}_{\textup{DSM}}(\theta))$ are trivially bounded.
To study $\ell_{\textup{DSM}}(\theta)$, we first note that

$\operatorname{Var}(X_{0}|X_{t})=\frac{\lambda_{t}\sigma^{2}}{\lambda_{t}+(1-\lambda_{t})\sigma^{2}}=:\rho_{t}.$

Write

$\begin{split}\mathbb{E}\left[\left.{\norm{s_{\theta}(t,X_{t})-\nabla\log
p(X_{t}|X_{0})}^{2}}\right|{X_{t}}\right]&\geq\operatorname{Var}\left(\left.{s_{\theta}(t,X_{t})-\nabla\log
p(X_{t}|X_{0})}\right|{X_{t}}\right)=\operatorname{Var}\left(\left.{-\frac{X_{t}-\sqrt{1-\lambda_{t}}X_{0}}{\lambda_{t}}}\right|{X_{t}}\right)\\\
&=\frac{1-\lambda_{t}}{\lambda_{t}^{2}}\rho_{t}\end{split}$

so

$\ell_{\textup{DSM}}(\theta)\geq\int_{0}^{T}\frac{1-\lambda_{t}}{\lambda_{t}^{2}}\rho_{t}\differential
t=\infty$

since $\rho_{t}\sim 2t$ as $t\to 0$. Concerning
$\hat{\nabla}\ell_{\textup{DSM}}(\theta)=2T(s_{\theta}(\tau,X_{\tau})-\nabla\log
p(X_{\tau}|X_{0}))\nabla_{\theta}s_{\theta}(\tau,X_{\tau})$, we have

$\begin{split}\mathbb{E}\left[\left.{\norm{\hat{\nabla}\ell_{\textup{DSM}}(\theta)}^{2}}\right|{X_{t},\tau=t}\right]&=4T^{2}\norm{\nabla_{\theta}s_{\theta}(t,X_{t})}^{2}\mathbb{E}\left[\left.{\norm{s_{\theta}(t,X_{t})-\nabla\log
p(X_{t}|X_{0})}^{2}}\right|{X_{t},\tau=t}\right]\\\ &\geq
4T^{2}\norm{\nabla_{\theta}s_{\theta}(t,X_{t})}^{2}\frac{1-\lambda_{t}}{\lambda_{t}^{2}}\rho_{t}\end{split}$

so

$\begin{split}\mathbb{E}\left[\norm{\hat{\nabla}\ell_{\textup{DSM}}(\theta)}^{2}\right]&=\frac{1}{T}\int_{0}^{T}\mathbb{E}\left[\left.{\norm{\hat{\nabla}\ell_{\textup{DSM}}(\theta)}^{2}}\right|{\tau=t}\right]\differential
t\geq
4T\int_{0}^{T}\mathbb{E}\left[\norm{\nabla_{\theta}s_{\theta}(t,X_{t})}^{2}\frac{1-\lambda_{t}}{\lambda_{t}^{2}}\rho_{t}\right]\differential
t\\\
&=4T\mathbb{E}\left[\int_{0}^{T}\norm{\nabla_{\theta}s_{\theta}(t,X_{t})}^{2}\frac{1-\lambda_{t}}{\lambda_{t}^{2}}\rho_{t}\differential
t\right].\end{split}$

The integral inside the last expectation is infinite whenever the event
$\norm{\nabla_{\theta}s_{\theta}(0,X_{0})}\neq 0$ holds, thanks to the
continuity of $\nabla_{\theta}s$ and the path $X_{[0,T]}$. Since
$\mathbb{E}\norm{\nabla_{\theta}s_{\theta}(t,X_{0})}^{2}>0$ by assumption,
that event has non-zero probability, which concludes the proof. ∎

##  Appendix D Experimental Details

In this section we give additional details and ablations relating to our
experimental results. We begin by providing details of the sampling tasks we
considered. We then provide details of our implementation and the baseline
methods. We finally provide additional ablation studies and results which
demonstrate the properties of our method.

###  D.1 Benchmarking targets

##### Gaussian

Here we consider the target $\pi(x)=\mathcal{N}(x;2.75,0.25^{2})$.

##### Mixture

This target was used in Arbel et al. (2021). It is an equally weighted mixture
of 6 bivariate Gaussian distributions with means
$\mu_{1}=(3.0,0.0),\mu_{2}=(-2.5,0.0),\mu_{3}=(2.0,3.0),\mu_{4}=(0.0,3.0),\mu_{5}=(0.0,-2.5),\mu_{6}=(3.0,2.0)$
and covariances $\Sigma_{1}=\Sigma_{2}=\big{(}\begin{smallmatrix}0.7&0.0\\\
0.0&0.05\end{smallmatrix}\big{)},\Sigma_{4}=\Sigma_{5}=\big{(}\begin{smallmatrix}0.05&0.0\\\
0.0&0.07\end{smallmatrix}\big{)},\Sigma_{3}=\Sigma_{6}=\big{(}\begin{smallmatrix}1.0&0.95\\\
0.95&1.0\end{smallmatrix}\big{)}$. This target is symmetric around $y=x$.

##### Funnel

This target was proposed by Neal (2003). Its density follows
$x_{0}\sim\mathcal{N}(0,\sigma_{f}^{2}),x_{1:9}|x_{0}\sim\mathcal{N}(0,\exp(x_{0})\mathbf{I})$,
with $\sigma_{f}=3$.

##### Logistic Regression

The Bayesian logistic regression model is defined by the prior distribution
$\theta\sim\mathcal{N}(0,\sigma^{2}\mathrm{I})$ and likelihood
$y|\theta,x\sim\text{Bernoulli}(\sigma(\theta^{T}x))$ where $\sigma$ is the
sigmoid function. We consider sampling the posterior $\theta|y,x$ on the
Ionosphere and Sonar datasets, which are of $35$ and $61$ dimensions
respectively.

##### Brownian Motion

In this task, we make noisy observations of a simple Brownian motion over 30
time steps. The model was introduced by Sountsov et al. (2020) and is defined
by the prior $\alpha_{\text{inn}}\sim\text{LogNormal}(0,2)$,
$\alpha_{\text{obs}}\sim\text{LogNormal}(0,2)$,
$x_{1}\sim\mathcal{N}(0,\alpha_{\text{inn}}^{2})$ and
$x_{i}\sim\mathcal{N}(x_{i-1},\alpha_{\text{inn}}^{2})$ for $i=2,...,30$. The
observation likelihood is given by
$y_{i}\sim\mathcal{N}(x_{i},\alpha_{\text{obs}}^{2})$ for $i=1,...,30$. The
goal is to sample the posterior distribution of
$\alpha_{\text{inn}},\alpha_{\text{obs}},x_{1},...,x_{30}|\\{y_{i}\\}_{i=1}^{10}\bigcup\\{y_{i}\\}_{i=21}^{30}$.
This task is in $32$ dimensions.

##### Log Gaussian Cox Process

The LGCP model (Møller et al., 1998) was developed for the analysis of spatial
data. The Poisson rate parameter $\lambda(x)$ is modelled on the grid using an
exponentially-transformed Gaussian process, and observations come from a
Poisson point process. The unnormalized posterior density is given directly by
$\gamma(x)=\mathcal{N}(x;\mu,K)\prod_{i\in[1:M]^{2}}\exp(x_{i}y_{i}-a\exp(x_{i}))$,
where $x_{i}$ are the points of a regular $M\times M$ grid. In our
experiments, we fit this model on the Pines forest dataset where $M=40$,
resulting in a problem in $1600$ dimensions.

###  D.2 Algorithmic details and hyperparameter settings

Here we give details of the algorithmic settings used in our experiments. We
first describe the considerations taken to ensure a fair comparison between
baselines, and then we detail exact hyperparameter settings.

Our method was implemented in Python using the libraries of JAX (Bradbury et
al., 2018), Haiku and Optax. Our implementation is available on
Github111https://github.com/angusphillips/particle_denoising_diffusion_sampler.
We used the open source code-bases of Arbel et al. (2021) to run the SMC and
CRAFT baselines and of Vargas et al. (2023) to run the PIS and DDS benchmarks,
both of which are also implemented in JAX.

In all experiments we used $2000$ particles to estimate the normalizing
constant.

##### Variational approximation

We used a variational approximation as the reference distribution for all
methods. We found that this was required for numerical stability of the
potential function $g_{t}(x_{t})$ in our method. We therefore used the same
variational approximation for all methods to ensure a fair comparison. Note
that PIS reverses a pinned brownian motion and thus the reference distribution
depends on the diffusion time span $T$ and the noise coefficient $\sigma$.
Since $\sigma$ affects the performance of the PIS algorithm itself, we tune
this parameter independently rather than setting this via the variational
approximation. The variational approximation was a mean-field variational
distribution i.e. a diagonal Gaussian distribution learnt by optimizing the
ELBO. We used $20,000$ optimisation steps ($50,000$ for the `Funnel` and
`Brownian` tasks) with the Adam optimizer (Kingma & Ba, 2015) and learning
rate $1e-3$. We did not use a variational approximation in the `Gaussian`
task.

##### Network architectures and optimizer settings

For the CRAFT baseline, we followed the flow network architectures and
optimizer settings given in Matthews et al. (2022), which are restated below
for completeness. For the diffusion-based methods (PDDS, DDS and PIS) we use
the same network architecture and optimizer settings for each method. The
neural network follows the PISGRAD network of Zhang & Chen (2022) with minor
adaptations. We use a sinusoidal embedding of 128 dimensions for the time
input. We use a 3-layer MLP with 64 hidden units per layer for the ‘smoothing’
network ($r_{\eta}(t)$ in our notation and $\mathrm{NN}_{2}(t)$ in Zhang &
Chen (2022)). For the main potential/score network ($\mathrm{N}_{\gamma}$ in
our notation and $\mathrm{NN}_{1}(t,x)$ in Zhang & Chen (2022)), we use a 2
layer MLP of 64 hidden units per layer to encode the 128-dimensional time
embedding. This is concatenated with the state input $x$ before passing
through a 3-layer MLP with 64 hidden units per layer, outputting a vector of
dimension $d$. In PDDS, we take the scalar product of this output with the
state input $x$ to approximate the potential function, while PIS and DDS use
the $d$-dimensional output to approximate the optimal control term. The
activation function is GeLU throughout. We train for $10,000$ iterations of
the Adam optimizer (Kingma & Ba, 2015) with batch size $300$ and a learning
rate of $1e-3$, which decays exponentially at a rate of $0.95$ every 50
iterations (with the exception of the `Funnel` task where we did not use any
learning rate decay).

The number of trainable parameters for each method and task can be found in
Table 1, along with training time and sampling time.

| Gaussian (16) | Mixture (16) | Funnel (32) | Brownian (16) | Ion (32) | Sonar (32) | LGCP (128)
---|---|---|---|---|---|---|---
PDDS1 |  37570 / 74 / 0.10 |  37764 / 81 / 0.13 |  39316 / 98 / 0.15 |  43584 / 79 / 0.19 |  44166 / 78 / 0.12 |  49210 / 77 / 0.12 |  347776 / 1404 / 1.86
CRAFT |  32 / 4 / 0.02 |  176608 / 26 / 0.09 |  6077440 / 178 / 0.18 |  1024 / 27 / 0.21 |  2240 / 38 / 0.09 |  3904 / 40 / 0.09 |  409600 / 3072 / 14.9
PIS |  37570 / 129 / 0.00 |  37764 / 167 / 0.01 |  39316 / 394 / 0.02 |  43854 / 176 / 0.01 |  44166 / 324 / 0.02 |  49210 / 326 / 0.01 |  347776 / 3931 / 0.64
DDS |  37570 / 136 / 0.00 |  37764 / 187 / 0.01 |  39316 / 381 / 0.02 |  43854 / 183 / 0.01 |  44166 / 338 / 0.01 |  49210 / 332 / 0.02 |  347776 / 3941 / 0.68
SMC | 0 / 0 / 0.02 | 0 / 0 / 0.07 | 0 / 0 / 0.17 | 0 / 0 / 0.20 | 0 / 0 / 0.09 | 0 / 0 / 0.09 | 0 / 0 / 14.6
Table 1: Number of trainable parameters / training time total (seconds) /
sampling time per 2000 particles (seconds). Timings are averaged over 3
training seeds.

##### Annealing and noise schedules

For the annealing based approaches (SMC and CRAFT) we used a geometric
annealing schedule with initial distribution the variational approximation as
described above. For the diffusion based approaches (PDDS, DDS and PIS) we
carefully considered the appropriate noise schedules for each method. Firstly,
we fix the diffusion time span at $T=1$ and adapt the discretization step size
depending on the number of steps $K$ of the experiment, i.e. $\delta=T/K$.
This choice is equivalent to the fixed discretization step and varying
diffusion time $T$ as considered by Vargas et al. (2023), up to the choice of
$\alpha_{\text{max}}$.

For PIS, the original work of Zhang & Chen (2022) used by default a uniform
noise schedule controlled by $\sigma$. Further, Vargas et al. (2023) were
unable to find a noise schedule which improved performance above the default
uniform noise schedule. As such as we stick with the uniform noise schedule
and tune the noise coefficient $\sigma$ by optimizing the ELBO objective over
a grid search up to a resolution of $0.02$.

For DDS, it was observed by Vargas et al. (2023) that controlling the
transition scale $\sqrt{\alpha_{k}}$ such that it goes smoothly to zero as
$k\to 0$ is critical to performance. To achieve this they choose a cosine-
based schedule
$\alpha_{k}^{1/2}=\alpha_{\text{max}}^{1/2}\cos^{2}\big{(}\frac{\pi}{2}\frac{1-k/K+s}{1+s}\big{)}$
for $s$ small ($0.008$ following Nichol & Dhariwal (2021)), which we term the
DDS cosine schedule. The parameter $\alpha_{\text{max}}$ is tuned such that
the noise at the final step in the reverse process is sufficiently small. We
found that this scheduler did indeed result in the best performance when
compared to a linear noise schedule ($\beta_{t}=\beta_{0}+\beta_{T}t/T$) or
the popular cosine schedule of Nichol & Dhariwal (2021)
($\lambda_{t}=1-\cos^{2}\big{(}\frac{\pi}{2}\frac{t/T+s}{1+s}\big{)}$). As
such we use the DDS cosine schedule and tune $\alpha_{\text{max}}$ by
optimizing the ELBO objective over a grid search to a resolution of $0.02$.

For our method PDDS, we obtained the best performance using the cosine
schedule of Nichol & Dhariwal (2021), which sets
$\lambda_{t}=1-\cos^{2}\big{(}\frac{\pi}{2}\frac{t/T+s}{1+s}\big{)}$ where we
recall that $\lambda_{t}=1-\exp(-2\int_{0}^{t}\beta_{s}ds)$, i.e. the variance
of the transition from $0$ to $t$. We provide an illustration of the benefits
of this schedule in the following section. In particular we found that the
alternative DDS cosine schedule did not improve performance and added the
complexity of tuning the parameter $\alpha_{\text{max}}$.

In summary, while each of the diffusion based approaches used different noise
schedulers, each was chosen to provide the best performance for the individual
approach and thus ensures a fair comparison.

####  D.2.1 SMC and CRAFT settings

We used 1 iteration of an HMC kernel with 10 leapfrog integrator steps as the
proposal distribution in the SMC and CRAFT baselines. We tuned the HMC step
sizes based on initial runs and obtained the step size schedules given below.
We performed simple resampling when the ESS dropped below $30\%$. We trained
CRAFT for $500$ iterations ($1000$ on `Funnel`) with a batch size of $2000$
and learning rate schedule detailed below. We also list the flow architecture
in each task. Our parameter settings differ to those in Matthews et al. (2022)
since we use the variational approximation, which results in larger MCMC step
sizes and smaller learning rates.

##### Gaussian

Step sizes $[0.7,0.7,0.5,0.4]$ linearly interpolated between times
$[0.0,0.25,0.5,1.0]$. CRAFT used a diagonal affine flow, with a learning rate
of $1e-2$.

##### Mixture

Step sizes $[0.5,0.5,0.5,0.3]$ linearly interpolated between times
$[0.0,0.25,0.5,1.0]$. CRAFT used a spline inverse autoregressive flow with 10
spline bins and a 3 layer autoregressive MLP of 30 hidden units per layer,
with a learning rate of $1e-3$.

##### Funnel

Step sizes $[1.0,0.9,0.8,0.7,0.6]$ linearly interpolated between times
$[0.0,0.25,0.5,0.75,1.0]$. CRAFT used an affine inverse autoregressive flow,
trained for 4000 iterations with a learning rate of $1e-3$.

##### Brownian

Step sizes $[0.8,0.8,0.7,0.6,0.5]$ linearly interpolated between times
$[0.0,0.25,0.5,0.75,1.0]$. CRAFT used a diagonal affine flow, with a learning
rate of $1e-3$.

##### Ion

Step sizes $[0.7,0.7,0.6,0.5,0.4]$ linearly interpolated between times
$[0.0,0.1,0.25,0.5,1.0]$. CRAFT used a diagonal affine flow, with a learning
rate of $1e-3$.

##### Sonar

Step sizes $[0.7,0.7,0.6,0.5,0.35]$ linearly interpolated between times
$[0.0,0.1,0.25,0.5,1.0]$. CRAFT used a diagonal affine flow, with a learning
rate of $1e-3$.

##### LGCP

Step sizes $[0.35,0.35,0.3,0.2]$ linearly interpolated between times
$[0.0,0.25,0.5,1.0]$. CRAFT used a diagonal affine flow, with a learning rate
of $1e-4$.

We used SMC with the above settings for $1000$ steps to estimate the ‘ground
truth’ normalizing constant on the Bayesian posterior targets.

####  D.2.2 DDS settings

Optimal settings for $\alpha_{\text{max}}$ are given in Table 2. Note that we
do not tune $\sigma$ as in Vargas et al. (2023) since we used a variational
approximation.

| Gaussian | Mixture | Brownian | Funnel | Ion | Sonar | LGCP
---|---|---|---|---|---|---|---
Steps | 1 | 1 | 1 | 2 | 2 | 2 | 8
1 | 0.86 | 0.28 | 0.76 | 0.60 | 0.68 | 0.68 | 0.74
2 | 0.86 | 0.28 | 0.76 | 0.68 | 0.80 | 0.82 | 0.62
4 | 1.00 | 0.36 | 0.84 | 0.68 | 0.74 | 0.78 | 0.60
8 | 0.96 | 0.52 | 0.80 | 0.60 | 0.64 | 0.64 | 0.44
16 | 0.82 | 0.54 | 0.72 | 0.64 | 0.52 | 0.50 | 0.26
Table 2: Optimal settings for $\alpha_{\text{max}}$. The number of steps for a
given entry is the base steps in the first row multiplied by the step
multiplier in the first column.

####  D.2.3 PIS settings

Optimal settings for $\sigma$ are given in Table 3. Note that we were unable
to obtain reasonable performance with PIS with only 1 step.

| Gaussian | Mixture | Brownian | Funnel | Ion | Sonar | LGCP
---|---|---|---|---|---|---|---
Steps | 1 | 1 | 1 | 2 | 2 | 2 | 8
1 | NA | NA | NA | 1.50 | 0.37 | 0.25 | 1.36
2 | 1.00 | 2.40 | 0.08 | 1.00 | 0.40 | 0.31 | 1.64
4 | 1.00 | 2.20 | 0.10 | 1.00 | 0.40 | 0.40 | 1.78
8 | 1.00 | 1.92 | 0.10 | 1.00 | 0.46 | 0.46 | 1.99
16 | 1.00 | 1.88 | 0.13 | 1.00 | 0.46 | 0.49 | 2.06
Table 3: Optimal settings for $\sigma$. The number of steps for a given entry
is the base steps in the first row multiplied by the step multiplier in the
first column. We were unable to tune PIS with one step size.

####  D.2.4 PDDS settings

No tuning of the cosine noise schedule was required. We performed systematic
resampling (Douc & Cappé, 2005) when the ESS dropped below $30\%$. We did not
perform any optional MCMC steps. In experiments where we used additional
iterations of potential approximation, we train each refinement for $10,000$
steps with a fresh instance of the learning rate schedule ($1e-3$ with
exponential decay at a rate of $0.95$ per 50 iterations). At each refinement
we initialise the potential approximation at it’s previous state, rather than
training from scratch at each refinement.

###  D.3 Ablation studies

##### Cosine scheduler

As stated above, we follow the cosine scheduler introduced by Nichol &
Dhariwal (2021). We found, as demonstrated in Figure 6 and Figure 7, that the
cosine schedule was effective in ensuring the forward SDE converges to the
target distribution while regularly spacing the ESS drops across the sampling
path.

Figure 6: Left: linear schedule, $\beta_{f}=8$. Middle: linear schedule:
$\beta_{f}=12$. Right: cosine schedule.

Figure 7: Left: linear schedule, $\beta_{f}=8$. Middle: linear schedule:
$\beta_{f}=12$. Right: cosine schedule.

##### Number of particles

In Figure 8 we show that the normalizing constant estimation error with the
naive potential approximation does decrease as the number of particles
increases. However, we could not eliminate the error before exceeding the
computer memory. We conclude that we must improve our initial potential
approximation to obtain feasible results, hence motivating the iterative
potential approximation scheme.

Figure 8: Samples and estimated normalizing constant from PDDS on
$\mathcal{N}(2.75,0.25^{2})$ target.

##### Guidance path of NN potential approximation

In Figure 9 we demonstrate that using a neural network to correct the naive
approximation has the effect of correcting the path of the guidance SDE. Here
we use a linear schedule for $\beta_{t}$ since the cosine schedule does not
allow for analytic roll-out of the SDE.

Figure 9: Guidance path using analytic potential function (left), naive
potential approximation (middle) and neural network potential approximation
(right).

###  D.4 Additional and uncurated results

In Figure 10 we display the normalising constant estimates on all tasks,
including the `Brownian` and `Ion` tasks which were omitted from the main text
for space. In Figure 11 we display the same results with outliers present.
Note that all methods are susceptible to erroneous over-estimation of the
normalising constant due to numerical errors an instability. It appears that
all methods are equally susceptible to this issue, with no single method
displaying more erroneous overestimation than the others. Results on the
Gaussian task are displayed separately in Figure 12.

Figure 10: Normalizing constant estimation results on all tasks. Outliers are
hidden for clarity. Each box consists of 2000 estimates, coming from 20
training seeds each with 100 evaluation seeds. Figure 11: Normalizing constant
estimation results on all tasks. Each box consists of 2000 estimates, coming
from 20 training seeds each with 100 evaluation seeds. Figure 12: Normalizing
constant estimation results on the Gaussian task. Each box consists of 2000
estimates, coming from 20 training seeds each with 100 evaluation seeds. Left:
outliers removed, right: uncurated.