# Learned Reference-based Diffusion Sampler for multi-modal distributions

Maxence Noble, Louis Grenioux, Marylou Gabrié, Alain Oliviero Durmus

► To cite this version:

Maxence Noble, Louis Grenioux, Marylou Gabrié, Alain Oliviero Durmus. Learned Reference-based Diffusion Sampler for multi-modal distributions. The Thirteenth International Conference on Learning Representations, Apr 2025, Singapour, France. hal-05164914

HAL Id: hal-05164914

https://hal.science/hal-05164914v1

Submitted on 16 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# LEARNED REFERENCE-BASED DIFFUSION SAMPLING FOR MULTI-MODAL DISTRIBUTIONS

Maxence Noble*, Louis Grenioux*, Marylou Gabrié & Alain Oliviero Durmus

CMAP, CNRS

École polytechnique, Institut Polytechnique de Paris
91120 Palaiseau, France

## ABSTRACT

Over the past few years, several approaches utilizing score-based diffusion have been proposed to sample from probability distributions, that is without having access to exact samples and relying solely on evaluations of unnormalized densities. The resulting samplers approximate the time-reversal of a noising diffusion process, bridging the target distribution to an easy-to-sample base distribution. In practice, the performance of these methods heavily depends on key hyperparameters that require ground truth samples to be accurately tuned. Our work aims to highlight and address this fundamental issue, focusing in particular on multi-modal distributions, which pose significant challenges for existing sampling methods. Building on existing approaches, we introduce *Learned Reference-based Diffusion Sampler* (LRDS), a methodology specifically designed to leverage prior knowledge on the location of the target modes in order to bypass the obstacle of hyperparameter tuning. LRDS proceeds in two steps by (i) learning a *reference* diffusion model on samples located in high-density space regions and tailored for multimodality, and (ii) using this reference model to foster the training of a diffusion-based sampler. We experimentally demonstrate that LRDS best exploits prior knowledge on the target distribution compared to competing algorithms on a variety of challenging distributions.

## 1 INTRODUCTION

We consider the problem of sampling from a probability density known up to a normalizing constant. More precisely, we consider a target distribution $\pi \in \mathcal{P}(\mathbb{R}^d)$ with probability density given by $x \mapsto \gamma(x)/\mathcal{Z}$, where $\gamma: \mathbb{R}^d \to \mathbb{R}_+$ can be pointwise evaluated and the normalizing constant $\mathcal{Z} = \int_{\mathbb{R}^d} \gamma(x)dx$ is intractable. This problem appears in a wide variety of applications such as Bayesian statistics (Liu & Liu, 2001; Kroese et al., 2011), statistical mechanics (Krauth, 2006) or molecular dynamics (Stoltz et al., 2010). In particular, we are interested in sampling from *multi-modal* distributions, i.e., distributions whose density admits multiple local maxima, called *modes*. Finding the modes of such distributions is a notoriously hard problem, yet, maybe surprisingly, even if the location of the modes is known, sampling $\pi$ remains a very challenging problem (Noé et al., 2019; Pompe et al., 2020; Grenioux et al., 2023). In this work, we aim to address this specific issue and will assume that we have access to the location of the modes as prior information on $\pi$. However, we do not assume to have access a *priori* to ground truth samples from $\pi$.

**Annealed MCMC.** Markov Chain Monte Carlo (MCMC) samplers are among the most popular approaches for sampling. In particular, gradient-based methods based on discretizations of Langevin or Hamiltonian dynamics (Roberts & Tweedie, 1996; Neal, 2012; Hoffman & Gelman, 2014) are guaranteed to be efficient for high-dimensional target distributions that are log-concave or satisfy functional inequalities (Dalalyan, 2017; Durmus & Moulines, 2017). However, when applied to multi-modal distributions, these MCMC methods suffer from high mixing time, that dramatically increases with the magnitude of the energy barriers between the modes, see Appendix I.3 for an illustration. To alleviate this issue, *annealing* strategies have been proposed. They consist in bridging an easy-to-sample *base* distribution $\pi^{\text{base}}$ to the *target* distribution $\pi$ via intermediate distributions

<footnote>*Both authors contributed equally.

Corresponding authors: {maxence.noble-bourillot, louis.grenioux}@polytechnique.edu</footnote>
$\{\pi_k\}_{k=0}^K$ with $K \ge 1$, $\pi_0 = \pi^{\text{base}}$ and $\pi_K = \pi$. This sequence should form a smooth path in the distribution space, so that sampling from $\pi_{k+1}$ starting from samples from $\pi_k$ is relatively simple. Typically, up to normalizing constants again, the densities of these annealed distributions are chosen as geometric interpolation between $\pi^{\text{base}}$ and $\gamma$, and are sampled either sequentially (Neal, 2001; Del Moral et al., 2006) or in parallel (Swendsen & Wang, 1986). However, these approaches may be prone to mode switching, i.e., the relative weight of the modes of the intermediate distributions may vary dramatically along the prescribed path (Woodard et al., 2009; Tawn et al., 2020; Syed et al., 2022). Hence, recovering the proportions of the different modes of $\pi$ is hard in practice.

Alternatively, the bridging distributions $\{\pi_k\}_{k=0}^K$ can be implicitly defined as the marginals of a diffusion process. A particular case of interest is the annealing path given by the time-reversal of a noising process, which is at the foundation of diffusion-based generative models (Sohl-Dickstein et al., 2015; Song et al., 2021; Ho et al., 2020), and has been proved to be robust to mode switching (Phillips et al., 2024). Following this approach, several diffusion-based samplers relying exclusively on Monte Carlo methods have recently shown promises for sampling from multi-modal distributions (Huang et al., 2024a;b; Grenioux et al., 2024). However, these non-parametric approaches require the repeated estimation of an intractable drift at each draw and at each annealing step.

**Annealed VI.** Due to their ability to amortize inference, Variational Inference (VI) techniques are a popular alternative to MCMC methods. Given a class of parameterized (variational) distributions that are easy to sample from, they aim to find the optimal parameters that minimize a fixed discrepancy metric with respect to $\pi$. Standard VI settings involve mean-field approximations (Wainwright & Jordan, 2008), mixture models (Arenz et al., 2023), and more recently normalizing flows (Rezende et al., 2014; Papamakarios et al., 2021), that typically optimize a reverse Kullback-Leibler (KL) objective. However, these methods struggle with multi-modal target distributions as they suffer from mode collapse (Jerfel et al., 2021; Blessing et al., 2024).

To tackle this issue, it has been proposed to follow the annealing paradigm by considering an extended target distribution with marginals corresponding to the $\{\pi_k\}_{k=0}^K$. This approach has been explored following explicit bridging paths, i.e., $\{\pi_k\}_{k=0}^K$ admit tractable unnormalized densities (Wu et al., 2020; Arbel et al., 2021; Geffner & Domke, 2021; Matthews et al., 2022; Doucet et al., 2022; Geffner & Domke, 2023; Vargas et al., 2024), but the resulting samplers may suffer from mode switching, as explained above. On the other hand, a recent class of annealed VI samplers using implicit diffusion-based paths has emerged (Tzen & Raginsky, 2019; Holdijk et al., 2023; Pavon, 2022; Zhang & Chen, 2022; Berner et al., 2023; Vargas et al., 2023a;b; Zhang & Chen, 2022; Berner et al., 2023; Vargas et al., 2023b; Akhound-Sadegh et al., 2024; Phillips et al., 2024) and seems well-suited for multi-modal target distributions. As we will further evidence in the present work, the promising results obtained with the latter methods are nevertheless mitigated by the need for careful tuning of their hyperparameters, which requires access to ground truth samples.

**Contributions** Within this context, we propose the *Learned Reference-based Diffusion Sampler* (LRDS), a variational diffusion-based sampler specifically designed for multi-modal distributions, which extends works from Zhang & Chen (2022); Vargas et al. (2023a); Richter et al. (2023):

* In the multi-modal setting, we highlight the sensitivity of previous variational diffusion-based methods with respect to their hyperparameters, which restrains their use in practice.

* To address this limitation, LRDS leverages the knowledge of the mode locations following two approaches: (a) GMM-LRDS relies on Gaussian Mixture Models and is well suited for a large variety of target distributions while being relatively lightweight, and (b) EBM-LRDS, a more computationally intensive scheme, takes advantage of Energy-Based Models for harder sampling problems.

* We show that GMM-LRDS and EBM-LRDS outperform previous diffusion-based samplers on challenging multi-modal settings.

**Notation.** For any measurable space $(\mathcal{X}, \mathcal{X})$, we denote by $\mathcal{P}(\mathcal{X})$ the space of probability measures defined on $(\mathcal{X}, \mathcal{X})$. Unless specified, if $\mathcal{X}$ is a topological space, $\mathcal{X}$ is defined as the Borel $\sigma$-field of $\mathcal{X}$. For any $T > 0$, we denote by $C_T = C([0, T], \mathbb{R}^d)$ the space of continuous functions from $[0, T]$ to $\mathbb{R}^d$ endowed with the uniform topology; hence $\mathcal{P}(C_T)$ is the set of continuous-time stochastic processes (or path measures) defined on $[0, T]$. For any $\mathbb{P} \in \mathcal{P}(C_T)$, we denote by $\mathbb{P}^R$ its time-reversal, defined such that if $(X_t)_{t \in [0, T]} \sim \mathbb{P}$, then $(X_{T-t})_{t \in [0, T]} \sim \mathbb{P}^R$. For any $t \in [0, T]$, we denote by $\mathbb{P}_t$ the marginal distribution of $\mathbb{P}$ at time $t$. The density of the Gaussian distribution with mean $\mathbf{m} \in \mathbb{R}^d$ and covariance $\Sigma \in \mathbb{R}^{d \times d}$ is denoted by $x \mapsto N(x; \mathbf{m}, \Sigma)$ and $\Delta_J$ refers to the $J$-dimensional simplex, i.e., $\Delta_J = \{(w_j)_{j=1}^J \in [0, 1]^J : \sum_{j=1}^J w_j = 1\}$ where $J \ge 1$.
Finally, for ease of understanding, we will denote *noising* diffusion processes by $(X_t)_{t \in [0,T]}$, driven by Brownian motion $(W_t)_{t \in [0,T]}$, while $(Y_t)_{t \in [0,T]}$ will refer to *denoising* diffusion processes, driven by Brownian motion $(B_t)_{t \in [0,T]}$. In this paper, $\Theta$ will refer to as the variational parameter space.

# 2 REFERENCE-BASED DIFFUSION SAMPLING

## 2.1 THEORETICAL FRAMEWORK

We first recall the variational diffusion-based framework formulated by Richter et al. (2023) in continuous time, that generalizes the approaches from Zhang & Chen (2022) and Vargas et al. (2023a).

**Time-reversed sampling process.** Consider a general *noising* diffusion process defined on $[0, T]$ by the linear SDE

$$dX_t = f(t)X_t dt + \sqrt{\beta(t)}dW_t, \quad X_0 \sim \pi, \qquad (1)$$

where the horizon $T > 0$ is fixed, $(W_t)_{t \in [0,T]}$ is a standard $d$-dimensional Brownian motion, $\beta: [0, T] \to (0, \infty)$ and $f: [0, T] \to \mathbb{R}$. We denote by $\mathbb{P}^*$ the path measure associated to the SDE (1). With appropriate choices of $f$ and $\beta$, this SDE admits explicit transition kernels and $X_T$ is approximately, or exactly, distributed according to an easy-to-sample *base* distribution $\pi^{\text{base}}$.

In the case where $\mathbb{P}_T^* = \pi^{\text{base}}$, we refer to this setting as an ‘exact’ noising scheme. One particular instance of this setting the Pinned Brownian Motion (PBM), considered by Tzen & Raginsky (2019); Zhang & Chen (2022); Vargas et al. (2023b), where $\pi^{\text{base}} = \delta_0$. Otherwise, when we only have $\mathbb{P}_T^* \approx \pi^{\text{base}}$, we define this setting as an ‘ergodic’ noising scheme. A well-known example is the Variance Preserving (VP) noising process for which $\pi^{\text{base}} = \mathbb{N}(0, \sigma^2 I_d)$ with $\sigma > 0$, proposed by Song et al. (2021) for score-based generative models and Vargas et al. (2023a) for sampling. Numerical experiments presented in this paper focus on these two schemes to provide fair comparison with previous methods, but the presented approach is applicable for an arbitrary noising scheme. We refer to Appendix C for more details.

Denote by $p_t^*$ the density of $\mathbb{P}_t^*$ w.r.t. the Lebesgue measure. Under mild assumptions on $\pi$, $f$ and $\beta$, see e.g., Cattiaux et al. (2023), the time-reversal of the noising process $\mathbb{P}^*$, i.e., the distribution of $(X_{T-t})_{t \in [0,T]}$ and denoted by $(\mathbb{P}^*)^R$, is associated to the SDE

$$dY_t = -f(T-t)Y_t dt + \beta(T-t)\nabla \log p_{T-t}^*(Y_t) dt + \sqrt{\beta(T-t)}dB_t, \quad Y_0 \sim \mathbb{P}_T^*, \qquad (2)$$

where $(B_t)_{t \in [0,T]}$ is another standard $d$-dimensional Brownian motion. By definition of the time-reversal, it holds that $Y_T \sim \pi$. Therefore, if we were able to simulate this diffusion process, we would obtain approximate samples from $\pi$. We therefore refer to $(Y_t)_{t \in [0,T]}$ as the *target process*. However, the scores $(\nabla \log p_t^*)_{t \in [0,T]}$ involved in the drift function of (2) are intractable in general. In addition, these scores cannot be estimated via usual score matching techniques used in generative modeling (Hyvärinen & Dayan, 2005; Vincent, 2011) since samples from $\pi$ are not available a *priori* in the setting at hand. This point has been addressed by Richter et al. (2023) who alternatively proposed a general variational approach on path measure space, which requires the definition of a reference process as described next.

**Variational reference-based approach.** This approach relies on a *reference process* that is solution of SDE (1) with initial condition $X_0 \sim \pi^{\text{ref}}$, where $\pi^{\text{ref}} \in \mathcal{P}(\mathbb{R}^d)$ is chosen such that the marginal scores of the associated path measure $\mathbb{P}^{\text{ref}}$ are *tractable*. Similarly to $\pi$, we assume that the probability density of $\pi^{\text{ref}}$ is known up to a multiplicative constant and is given by $x \mapsto \gamma^{\text{ref}}(x)/\mathcal{Z}^{\text{ref}}$, where $\gamma^{\text{ref}}: \mathbb{R}^d \to \mathbb{R}_+$ and $\mathcal{Z}^{\text{ref}} = \int_{\mathbb{R}^d} \gamma^{\text{ref}}(x)dx$ is *not necessarily tractable*. Finally, for any $t \in [0, T]$, we denote the marginal scores $s_t^{\text{ref}} = \nabla \log p_t^{\text{ref}}$, with $p_t^{\text{ref}}$ being the density of $\mathbb{P}_t^{\text{ref}}$ w.r.t. the Lebesgue measure. Based on the tractable reference scores, the SDE describing the target process (2) can be rewritten as

$$dY_t = -f(T-t)Y_t dt + \beta(T-t)\{s_{T-t}^{\text{ref}}(Y_t) + g_{T-t}(Y_t)\}dt + \sqrt{\beta(T-t)}dB_t, \quad Y_0 \sim \mathbb{P}_T^*, \qquad (3)$$

where $g_t = \nabla \log p_t^*/p_t^{\text{ref}}$ is now the only intractable term. This formulation of the time-reversed SDE has already been largely used in the diffusion model literature, especially for conditional generative models; see e.g., Dhariwal & Nichol (2021). In this specific context, $g_t$ is called a guidance term. From a probabilistic perspective, $g_t$ can be interpreted as a Doob-h transform control, see Appendix I.1. Exploiting the reference-based formulation given by (3), Richter et al. (2023) propose to estimate the guidance terms $(g_t)_{t \in [0,T]}$ rather than the target scores $(\nabla \log p_t^*)_{t \in [0,T]}$.
<table><caption>Table 1: Connection to prior works. Here, $E^\varphi : [0, T] \times \mathbb{R}^d \to \mathbb{R}$ refers to a neural network.</caption><thead><tr><th>Method</th><th>Noising</th><th>γ<sup>ref</sup>(x)</th><th>Z<sup>ref</sup></th><th>s<sub>t</sub><sup>ref</sup>(x)</th><th>Reference parameters</th></tr></thead><tbody><tr><td>PIS (Zhang & Chen, 2022)</td><td>PBM(σ)</td><td>N(x; 0, σ<sup>2</sup>TI<sub>d</sub>)</td><td>1</td><td>-x/σ<sup>2</sup>(T - t)</td><td>σ (tuned)</td></tr><tr><td>DDS (Vargas et al., 2023a)</td><td>VP(σ)</td><td>N(x; 0, σ<sup>2</sup>TI<sub>d</sub>)</td><td>1</td><td>-x/σ<sup>2</sup></td><td>σ (tuned)</td></tr><tr><td>GMM-LRDS (Section 3.2)</td><td>arbitrary</td><td>Σ<sub>j=1</sub><sup>J</sup> w<sub>j</sub>N(x; m<sub>j</sub>, Σ<sub>j</sub>)</td><td>1</td><td>analytical</td><td>{w<sub>j</sub>, m<sub>j</sub>, Σ<sub>j</sub>}<sub>j=1</sub><sup>J</sup> (learned)</td></tr><tr><td>EBM-LRDS (Section 3.3)</td><td>arbitrary</td><td>exp(-E<sup>φ</sup>(t = 0, x))</td><td>unknown</td><td>-∇<sub>x</sub>E<sup>φ</sup>(t, x)</td><td>φ (learned)</td></tr></tbody></table>

To this end, the authors consider the class of *variational* path measures ($\mathbb{P}^\theta$)$_{\theta \in \Theta} \subset \mathcal{P}(\mathcal{C}_T)$ where for any $\theta \in \Theta$, $\mathbb{P}^\theta$ is associated to the diffusion process

$$dY_t^\theta = -f(T-t)Y_t^\theta dt + \beta(T-t)\{s_{T-t}^\text{ref}(Y_t^\theta) + g_{T-t}^\theta(Y_t^\theta)\}dt + \sqrt{\beta(T-t)}dB_t, \quad Y_0^\theta \sim \pi^\text{base}, \quad (4)$$

where $g^\theta: [0, T] \times \mathbb{R}^d \to \mathbb{R}^d$ is typically a neural network. In particular, choosing $\pi^\text{ref}$ as a Gaussian distribution recovers exactly the previous VI methods *Path Integral Sampler* (PIS) (Zhang & Chen, 2022) and *Denoising Diffusion Sampler* (DDS) (Vargas et al., 2023a), see Table 1. Furthermore, Richter et al. (2023) propose to minimize the following LV-based objective

$$\mathcal{L}_{\text{LV}}(\theta) = \text{Var} \left[ \log \left\{ \frac{d\mathbb{P}^\theta}{d(\mathbb{P}^*)^R} (Y_{[0,T]}^\hat{\theta}) \right\} \right], \quad Y_{[0,T]}^\hat{\theta} \sim \mathbb{P}^\theta, \quad (5)$$

where $\hat{\theta}$ is a detached version of $\theta$, meaning that gradients with respect to $\theta$ will not see $\hat{\theta}$. In contrast to the reverse KL divergence (previously used in PIS and DDS), which is notably known to suffer from mode collapse (Midgley et al., 2023), the LV divergence (Nüsken & Richter, 2021) has the benefits to improve mode exploration, to avoid costly computations of the loss gradients and to have zero variance at the optimal solution (Richter et al., 2023, Section 2.3). Finally, the LV loss (5) can be further explicited, as detailed in the following proposition.

**Proposition 1.** Assume that $\mathbb{P}_T^* = \mathbb{P}_T^{\text{ref}} = \pi^{\text{base}}$ and there exists $\theta^* \in \Theta$ such that $g_t^{\theta^*} = g_t$. Then, the loss defined in (5) achieves optimal solution at $\theta^*$ and, setting $\varrho = \log(\gamma^{\text{ref}}/\gamma)$, it simplifies as

$$\mathcal{L}_{\text{LV}}(\theta) = \text{Var} \left[ \int_{0}^{T} \frac{\beta(T-t)}{2} \left\| g_{T-t}^{\theta}(Y_t^{\hat{\theta}}) \right\|^2 dt + \int_{0}^{T} \sqrt{\beta(T-t)} g_{T-t}^{\theta}(Y_t^{\hat{\theta}})^{\top} dB_t + \varrho(Y_t^{\hat{\theta}}) \right]. \quad (6)$$

This result is an adaptation of (Richter et al., 2023, Lemma 3.1), which proof is restated in Appendix D.1 for completeness. We also emphasize that the assumption made in Proposition 1 is only needed for ergodic noising schemes. In practice, Richter et al. (2023) only consider two specific versions of this continuous-time variational loss, where $\pi^\text{ref}$ is determined by the PIS and DDS settings. In the following, we will refer to these extensions of PIS and DDS as LV-PIS and LV-DDS. In contrast, we keep a general perspective, by letting $\pi^\text{ref}$ completely arbitrary, and describe our proposition of practical implementation of a Reference-based Diffusion Sampler (RDS) in this context.

## 2.2 REFERENCE-BASED DIFFUSION SAMPLING IN PRACTICE

**Discrete-time setting.** In practice, the parameterized process $(Y_t^{\hat{\theta}})_{t \in [0, T]}$ cannot be simulated exactly and SDE (4) can only be numerically solved with a small step-size using, for example, the Euler-Maruyama (EM) or Exponential Integration (EI) (Durmus & Moulines, 2015) discretization schemes, see Appendix B for additional details. Consider a time discretization of the interval $[0, T]$ given by an increasing sequence of timesteps $\{t_k\}_{k=0}^K$ such that $t_0 = 0, t_K = T$ and $K \ge 1$. Given an arbitrary reference distribution $\pi^\text{ref}$, we propose to approximate the continuous-time objective defined in (6) by the following discrete time RDS objective

$$\mathcal{L}_{\text{RDS}}(\theta) = \text{Var} \left[ \sum_{k=0}^{K-1} w_k g_{T-t_k}^{\theta}(Y_k)^{\top} \left\{ g_{T-t_k}^{\hat{\theta}}(Y_k) - \frac{1}{2} g_{T-t_k}^{\theta}(Y_k) \right\} + \sum_{k=0}^{K-1} \sqrt{w_k} g_{T-t_k}^{\theta}(Y_k)^{\top} Z_k + \varrho(Y_k) \right], \quad (7)$$

where $\{Z_k\}_{k=0}^{K-1}$ are independently distributed according to $N(0, I_d)$, $\{Y_k\}_{k=0}^K$ is recursively defined by $Y_0 \sim \pi^{\text{base}}$ and for any $k \in \{0, \dots, K-1\}$,

$$Y_{k+1} = a_k Y_k + b_k \{s_{T-t_k}^{\text{ref}}(Y_k) + g_{T-t_k}^{\hat{\theta}}(Y_k)\} + \sqrt{c_k} Z_k, \quad (8)$$
with $\hat{\theta}$ being a detached version of $\theta$ and $\{w_k, a_k, b_k, c_k\}_{k=0}^{K-1}$ being tractable coefficients that depend on the choice of discretization and noising schemes. In particular, if $K$ is sufficiently large, $\{Y_k\}_{k=0}^K$ has approximately the same distribution as $\{Y_{t_k}^{\hat{\theta}}\}_{k=0}^K$. We refer to Appendix D.2 for a detailed explanation on the computation of our loss and explicit values of the coefficients. We emphasize that this variational objective can be implemented with any reference distribution. In comparison, the practical loss functions<fnref n="1" /> exhibited by Richter et al. (2023) are only designed for LV-PIS and LV-DDS, and may be seen as particular instances of (7).

**Training and sampling with RDS.** The training stage of RDS consists in running $N$ iterations of Stochastic Gradient Descent (SGD) to minimize the discrete-time objective (7), for some $N \ge 1$. Like other variational diffusion-based objectives, the RDS loss is simulation-based, meaning that, the whole process $\{Y_k\}_{k=0}^K$ defined in (8) must be simulated at each SGD step. The complete training procedure of RDS is summarized in Algorithm 2, stated in Appendix A. After training, the final parameter $\theta^*$ is fixed, and approximate samples from $\pi$ are obtained by simulating $\{Y_k\}_{k=0}^K$ from $g^{\theta^*}$, using (8), and taking $Y_K$.

# 3 LEARNED RDS FOR MULTI-MODAL DISTRIBUTIONS

In this section, we highlight the crucial role played by the reference distribution $\pi^{\text{ref}}$ in RDS when targeting multi-modal distributions. We first show in Section 3.1 the limitations of previous variational reference-based methods and give intuition on how $\pi^{\text{ref}}$ should be chosen. Based on this observation, we introduce *Learned Reference-based Diffusion Sampler* (LRDS), a novel diffusion-based sampler defined as a practical instance of RDS, where $\pi^{\text{ref}}$ is learned from prior knowledge on $\pi$, namely the location of their modes. We present two approaches based on Gaussian Mixture Models (GMM-LRDS), see Section 3.2, and Energy-Based Models (EBM-LRDS), see Section 3.3.

## 3.1 ON THE IMPORTANCE OF THE REFERENCE DISTRIBUTION

Choosing a Gaussian distribution for $\pi^{\text{ref}}$ ensures that the entire reference process is Gaussian. Therefore, the scores $(s_l^{\text{ref}})_{t \in [0,T]}$ can be analytically computed. We will refer to this case as Gaussian RDS (G-RDS). In this setting, previous variational reference-based methods (PIS, DDS and their log-variance version) focused on using isotropic centered Gaussian distributions as reference distributions (see Table 1), i.e., $\pi^{\text{ref}} = N(0, \sigma^2 I_d)$ where the hyperparameter is $\sigma \in (0, \infty)$. This parameter is fixed before running the optimization over $\theta$ and a priori disconnected from the target distribution. To assess the sensitivity of the algorithms with respect to the hyperparameter $\sigma$, we run LV-DDS and LV-PIS on a simple multi-modal distribution – a 16-dimensional bi-modal Gaussian mixture with two distinct weights – and examine the accuracy of estimation of these weights (see Figure 1, left). We observe that the performance is highly dependent on the value of $\sigma$, with an optimal accuracy achieved when the isotropic variance is directly related to statistics of the target distribution; see Appendix I.2 for details. However, it is hard to estimate this quantity without access to ground truth samples.

**Towards a multi-modal reference distribution.** The previous experiment suggests that choosing a reference distribution that is close to the target distribution improves the performance of the VI sampler. To validate this intuition, we target the same Gaussian mixture as above using RDS with PBM and VP noising schemes, setting $\pi^{\text{ref}}$ as a Gaussian mixture with the same components as the target with their *relative weights* as hyperparameter (see Figure 1, right). We find that this design of $\pi^{\text{ref}}$, with a well adjusted Gaussian mixture reference, enables RDS to robustly recover the information of the relative weight of the target modes. These observations motivate us to adopt the following paradigm for RDS: given a multi-modal target distribution $\pi$, we aim to design $\pi^{\text{ref}}$ close to $\pi$, sharing the same multi-modality characteristics. Below, we present *Learned Reference-based Diffusion Sampler* (LRDS), a complete sampling methodology relying on the RDS framework that leverages prior knowledge on the target modes to learn a well-suited reference distribution.

<footnote n="1">¹Although the authors do not provide any expression of a discrete-time version of (6) in their paper, a discrete-time loss is implemented for LV-PIS and LV-DDS in their official codebase available at https://github.com/juliusberner/sde_sampler.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: **Illustration of the decisive role of the reference distribution.** Here, we target a 16-dimensional Gaussian mixture with two modes, that have respective weights $w = 2/3$ and $1 - w = 1/3$, and display the estimation error of $w$ with different methods. (Left): Results for LV-PIS and LV-DDS when varying the value of their hyperparameter $\sigma$ (which directly determines $\pi^{\text{ref}}$ as shown in Table 1). The green dotted line represents the optimal variance for Gaussian approximation of $\pi$, see Appendix I.2 for related computations. (Right) Results for RDS in PBM and VP settings when setting $\pi^{\text{ref}}$ as a Gaussian mixture with the same modes as $\pi$, but $w$ is replaced by $w_{\text{ref}}$. Details on the design of this experiment are given in Appendix H.1.</figcaption></figure>

**LRDS pipeline.** Our algorithm LRDS proceeds in three main steps:

(a) **Obtain reference samples.** Given a standard MCMC sampler which targets $\pi$, such as the *Metropolis-Adjusted Langevin Algorithm* (MALA) (Roberts & Tweedie, 1996), we simulate multiple Markov chains that are initialized in the target mode locations. We refer to the obtained samples as the *reference samples* and denote by $\hat{\pi}^{\text{ref}}$ the corresponding empirical distribution. As recalled in Appendix I.3, in presence of high energy barriers between the modes, $\hat{\pi}^{\text{ref}}$ fails at estimating the global energy landscape. Moreover, $\hat{\pi}^{\text{ref}}$ does not admit a density w.r.t. the Lebesgue measure.

(b) **Define the reference process.** In this stage, we aim at (i) learning $\pi^{\text{ref}}$ as a continuous approximation of $\hat{\pi}^{\text{ref}}$, and (ii) computing the marginal scores at times $\{T-t_k\}_{k=0}^{K-1}$ of the induced reference process $\mathbb{P}^{\text{ref}}$, i.e., the noising process defined by SDE (1), where $X_0 \sim \pi^{\text{ref}}$. In practice, we set $\pi^{\text{ref}}$ either as a Gaussian Mixture Model, which relies on a light parameterization, or as an Energy-Based Model, which is more expressive at the cost of higher computational budget.

(c) **Run the RDS variational optimization.** Once the reference distribution is learned and the reference scores are computed, we run the RDS procedure described in Section 2.2.

In the next sections, we describe two practical implementations of stage (b) in LRDS using different parameterizations of the reference distribution.

## 3.2 GAUSSIAN MIXTURE MODEL LRDS

We first propose to set $\pi^{\text{ref}}$ as a Gaussian mixture to integrate multi-modality in RDS. In particular, we have $\mathcal{Z}^{\text{ref}} = 1$ and any reference density is parameterized as $\gamma^{\text{ref}}(\cdot) = \sum_{j=1}^{J} w_j N(\cdot; \mathbf{m}_j, \mathbf{\Sigma}_j)$ where $J \ge 1$ is the number of reference modes, $\{w_j\}_{j=1}^{J} \in \Delta_J$ are the reference weights, $\mathbf{m}_j \in \mathbb{R}^d$ and $\mathbf{\Sigma}_j \in \mathbb{R}^{d \times d}$ are respectively the mean and the covariance of the j-th reference mode. In this setting, each marginal of $\mathbb{P}^{\text{ref}}$ is also a Gaussian mixture whose parameters can be simply deduced from $\gamma^{\text{ref}}$, making their score fully tractable, see Appendix B. Therefore, the only difficulty at stage (b) lies in the estimation of the reference parameters $\varphi = \{w_j, \mathbf{m}_j, \mathbf{\Sigma}_j\}_{j=1}^{J}$ for some $J \ge 1$ fixed in advance. If $J > 1$, we refer to this version of the algorithm as *Gaussian Mixture Model LRDS (GMM-LRDS)*. In this case, we learn $\varphi$ by running the Expectation Minimization (EM) algorithm (Dempster et al., 1977) with samples from $\hat{\pi}^{\text{ref}}$. If $J = 1$, LRDS reduces to *Gaussian LRDS (G-LRDS)*, and we learn parameters $(\mathbf{m}_1, \mathbf{\Sigma}_1)$ by Maximum Likelihood (ML) estimation. In practice, we observe that setting $J$ equal or larger to the number of target modes can lead to better performance. We summarize the whole GMM-LRDS methodology in Algorithm 3, provided in Appendix A, and illustrate in Figure 2 the behaviours of GMM-LRDS and G-LRDS for a bi-modal target distribution. Our results show that choosing a Gaussian reference distribution leads to mode collapse, while defining $\pi^{\text{ref}}$ as in GMM-LDRS enables RDS to recover the relative weights accurately.
<figure><img src="image_3.png" /><figcaption>Figure 2: **Comparison between G-LRDS and GMM-LRDS.** Here, the target distribution is the same 16-dimensional Gaussian mixture as in Figure 1 (top left), see Appendix H.1 for more details. For illustration purpose, projections along the first two coordinates are used. In each cell, the value of 'Mode weight' refers to the effective weight of the left mode. Reference samples (bottom left) are obtained by running MALA sampler initialized in both target modes: each color (orange/purple) depicts one MALA Markov chain. In particular, none of them mixes between the modes. Running G-LRDS with an ML-estimated Gaussian reference (top middle) leads to mode collapse (bottom middle). Conversely, GMM-LRDS with an EM-estimated Gaussian mixture reference (top right) appropriately recovers the target distribution and the true mode proportions (bottom right).</figcaption></figure>

3.3 ENERGY-BASED MODEL LRDS

Although GMM-LRDS efficiently introduces multi-modality in the reference distribution while con-
serving the tractability of the reference scores, we anticipate that in some cases, a Gaussian mixture
may not provide a correct approximation of $\hat{\pi}^{\text{ref}}$, regardless of the number of reference modes $J$, see Appendix I.4 for illustrations. To provide more flexibility, we propose a second version of LRDS where the reference distribution is parameterized by an *Energy-Based Model* (EBM), i.e., $\gamma^{\text{ref}}(x) = \exp(-E^{\varphi}(x))$ where $E^{\varphi}: \mathbb{R}^d \rightarrow \mathbb{R}$ is a neural network with parameters $\varphi$ and the normalizing constant $Z^{\text{ref}} = \int_{\mathbb{R}^d} \exp(-E^{\varphi}(x))dx$ is *intractable*. In contrast to GMM-LRDS, the corresponding reference scores cannot be computed analytically anymore; hence, we suggest to *estimate* them with neural networks. In the following, we will denote $\pi^{\text{ref}}$ by $p^{\varphi}$ to insist on its parametric nature.

At first sight, it seems natural to learn the density $\gamma^{\text{ref}}$ and the scores $(s_t^{\text{ref}})_{t \in [0,T]}$ independently from each other. Indeed, following the diffusion model literature (Song et al., 2021; Karras et al., 2022), the reference scores could be easily learned using a *Score Matching* (SM) objective with samples from $\hat{\pi}^{\text{ref}}$ (Hyvärinen & Dayan, 2005; Vincent, 2011). On the other hand, one could learn $p^{\varphi}$ by maximizing the standard *Maximum Likelihood* (ML) objective $\mathcal{L}^{\text{ML}}: \varphi \mapsto \mathbb{E}[\log p^{\varphi}(X)]$, where $X \sim \hat{\pi}^{\text{ref}}$, with corresponding gradient given by $\varphi \mapsto \mathbb{E}[\nabla_{\varphi} E^{\varphi}(X^{-})] - \mathbb{E}[\nabla_{\varphi} E^{\varphi}(X^{+})]$, where $X^{-} \sim p^{\varphi}$ and $X^{+} \sim \hat{\pi}^{\text{ref}}$. Nevertheless, computing gradients of objective $\mathcal{L}^{\text{ML}}$ requires to sample from $p^{\varphi}$ (*negative sampling*), which is a well known hurdle in EBM training since $p^{\varphi}$ is expected to be as multi-modal as $\hat{\pi}^{\text{ref}}$, and thus $\pi$.

Therefore, we rather suggest to learn a path of parametric distributions $(p_t^\varphi)_{t \in [0,T]}$ as the marginal distributions of $(\hat{X}_t^{\text{ref}})_{t \in [0,T]}$, defined as the noising process induced by SDE (1), starting at $\hat{X}_0^{\text{ref}} \sim \hat{\pi}^{\text{ref}}$. To do so, we set $(p_t^\varphi)_{t \in [0,T]}$ as a *multi-level EBM*, i.e., $p_t^\varphi(x) = \exp(-E^\varphi(t, x))/Z_t^\varphi$, where $E^\varphi: [0, T] \times \mathbb{R}^d \to \mathbb{R}$ is now a *time-dependent* neural network with parameters $\varphi$ and the normalizing constants $Z_t^\varphi = \int_{\mathbb{R}^d} \exp(-E^\varphi(t, x))dx$ are still intractable. In the RDS framework, this amounts to consider $\gamma^{\text{ref}}(x) = \exp(-E^\varphi(0, x))$ and $s_t^{\text{ref}}(x) = -\nabla_x E^\varphi(t, x)$. To learn this multi-level EBM, we propose to maximize the integrated ML objective $\varphi \mapsto \int_0^T \mathbb{E}[\log p_t^\varphi(\hat{X}_t^{\text{ref}})]dt$. In this case, we can leverage the correlations between the single-level EBMs to alleviate their individual issue of negative sampling, a strategy that has already been investigated in several works (Gao et al., 2021; Zhu et al., 2024; Zhang et al., 2023).
Here, since $(\hat{X}_t^{\text{ref}})_{t \in [0,T]}$ defines a path of increasingly simpler distributions, annealed MCMC samplers can be conveniently used to sample from the multi-level EBM densities. This allows us to kill two birds in one stone: (i) we obtain negative samples for each single-level EBM, which is needed to compute the gradient of the ML objective, and (ii) we overcome the individual sampling issues at every level thanks to annealing. This method is completely detailed in Algorithm 11 of Appendix F together with previous literature on multi-level EBMs. Morever, Appendix E provides details on annealed MCMC samplers.

This version of LRDS called *Energy-Based Model LRDS* (EBM-LRDS) is summarized in Algorithm 4 of Appendix A. In Figure 3, we illustrate the superiority of EBM-LRDS over GMM-LRDS for a target distribution $\pi$ that exhibits complex geometry. While the GMM fails at capturing the local energy landscape of $\pi$, which results in a poor performance of GMM-LRDS, EBM-LRDS captures well the target distribution since the reference EBM recovers the complex geometry of the target.

<figure><img src="image_4.png" /><figcaption>Figure 3: **Comparison between GMM-LRDS and EBM-LRDS in a multi-modal setting.** Here, we target the 2-dimensional Rings distribution, which has 3 *unbalanced* modes represented by the rings. **(Left)**: Target density (top) and exact samples (bottom). **(Middle)**: 16-component GMM reference distribution (top) and resulting GMM-LRDS samples (bottom). **(Right)**: EBM reference distribution (top) and resulting EBM-LRDS samples (bottom).</figcaption></figure>

# 4 RELATED WORKS

In the past years, numerous parametric methods have been proposed to approximate a diffusion process $(Y_t)_{t \in [0,T]} \sim \mathbb{P}$ induced by an SDE with intractable drift, such as (2), bridging an easy-to-sample distribution $\pi^{\text{base}}$ to a target distribution $\pi$. Based on deep learning techniques, they suggest to learn a path measure $\mathbb{P}^{\theta}$, induced by a *neural* SDE that admits a neural network (parameterized by $\theta \in \Theta$ to be optimized) as drift function. Below, we review two main classes of such algorithms.

**Variational diffusion-based methods.** To optimize $\theta$, several approaches have adopted a variational formulation on path measure space, which consists in minimizing a divergence-based loss with samples from $\mathbb{P}^{\theta}$. We adopt the same perspective in the present work. For instance, a recent line of works has defined $\mathbb{P}$ as the time-reversal of a noising diffusion process, such as PIS (Zhang & Chen, 2022), DDS (Vargas et al., 2023a) and *Time-reversed Diffusion Sampler* (DIS) (Berner et al., 2023). In these methods, the neural network used in $\mathbb{P}^{\theta}$ is crucially required to be *target-informed*, i.e., it is parameterized with the score of the target distribution. For a large number of variational training steps, this may result in a costly procedure, since each training step requires a full simulation of $\mathbb{P}^{\theta}$. While those algorithms were originally implemented with the largely used reverse KL divergence, Richter et al. (2023) demonstrate the benefits of using the Log-Variance (LV) divergence to avoid mode collapse in the variational optimization stage. However, even combined with a LV-based loss, these variational methods still require a target-informed parameterization. On the other hand, Vargas et al. (2024) propose a different perspective by defining $\mathbb{P}$ as a controlled version of an *annealed Langevin* diffusion, that is expected to follow a prescribed path of tractable marginal densities. Here, both $\mathbb{P}$ and $\mathbb{P}^{\theta}$ are induced by neural SDEs. The authors present two versions of their algorithm, *Controlled Monte Carlo Diffusion* (CMCD), using either reverse KL or LV divergence. For clear comparison with the RDS discrete-time setting presented in Section 2.2, we describe all of these approaches (PIS, DDS, DIS, CMCD) under the discrete time scope in Appendix D.2.

**Adaptive diffusion-based approaches.** To alleviate the computational difficulties of divergence-based losses, recent works have proposed to learn $\mathbb{P}^{\theta}$ with an adaptive procedure (Phillips et al., 2024; Akhound-Sadegh et al., 2024). These methods iterate two steps which consist of (a) obtaining approximate target samples by sampling from $\mathbb{P}^{\theta}$ and (b) optimizing $\theta$ using those samples via learning techniques usually restricted to generative modeling. In particular, Akhound-Sadegh et al. (2024) present *Iterated Denoising Energy Matching* (iDEM), where $\mathbb{P}$ corresponds to the time-reversal of a Variance-Exploding noising diffusion process (Song et al., 2021). In their setting, stage (a) is conducted by running the SDE induced by $\mathbb{P}^{\theta}$ and stage (b) relies on a novel energy-matching loss which directly depends on $\pi$. On the other hand, the *Particle Denoising Diffusion Sampler*
(PDDS) (Phillips et al., 2024) (a) introduces a SMC-based scheme when sampling from $P^θ$, see Appendix E for more details, and (b) implements an extension of the Target Score Matching loss (De Bortoli et al., 2024). When it comes to practice, we find that both of these algorithms face significant limitations: while PDDS performance is highly sensitive to the choice of the initial $P^θ$, the iDEM energy-matching loss suffers from very high variance.

# 5 NUMERICAL EVALUATION OF RDS

To validate our approach, we compare GMM-LRDS and EBM-LRDS on a variety of multi-modal distributions against the following annealed methods: (a) **annealed MCMC methods** – *Sequential Monte Carlo* (SMC) (Del Moral et al., 2006) and *Replica Exchange* (RE) (Swendsen & Wang, 1986); (b) **variational diffusion-based methods**, implemented with the LV loss – LV-PIS (Zhang & Chen, 2022), LV-DDS (Vargas et al., 2023a), LV-DIS (Berner et al., 2023) and LV-CMCD (Vargas et al., 2024); (c) **adaptive diffusion-based approaches** – iDEM (Akhound-Sadegh et al., 2024) and PDDS (Phillips et al., 2024). To assess the performance of each sampler in multi-modal settings, we will evaluate how well the obtained samples are able to recover the weights of the target modes<fnref n="2" />. The details of each target distribution can be found in Appendix H.1.

**General experimental setting.** To ensure fair comparison with previous approaches, we make sure that all competing methods are as informed as LRDS of prior knowledge on the target distribution. More specifically: we set $\pi^{\text{base}}$ as a Gaussian approximation of $\hat{\pi}^{\text{ref}}$ in SMC, RE or CMCD; we choose $\sigma$ based on a Gaussian isotropic approximation of $\hat{\pi}^{\text{ref}}$ for PIS, DDS or DIS (see Section 3.1); we standardize the target distributions using the empirical mean and variances of $\hat{\pi}^{\text{ref}}$ for iDEM and PDDS. Additionally, we pre-fill iDEM’s training buffer with samples from $\hat{\pi}^{\text{ref}}$. Note that all competing methods use the score of the target distribution, either in their training procedure when computing the loss (iDEM, PDDS) or in the sampling procedure: through a target-informed neural network parameterization in PIS, DDS and DIS; through the target-informed base drift in CMCD; through the MCMC gradient steps in SMC or RE. In contrast, LRDS only requires evaluations of the target density, which makes it an interesting alternative in settings where the score of $\pi$ is expensive to compute. We refer to Appendix H.3 for complete details of the implementation of these methods.

**High-dimensional Gaussian mixtures.** We first consider a synthetic but challenging setting, where $\pi$ is a bi-modal Gaussian mixture whose modes are $N(-1_d, \Sigma_1)$, with weight $w_1 = 2/3$, and $N(1_d, \Sigma_2)$, with weight $w_2 = 1/3$, where $1_d$ is the $d$-dimensional vector with all components equal to 1, and for $i \in \{1, 2\}$, $\Sigma_i \in \mathbb{R}^{d \times d}$ is a diagonal positive matrix with conditioning number equal to 100. For each method, the ground truth weight $w_1$ is estimated by a Monte Carlo estimator $\hat{w}_1$. We report the estimation error $|w_1 - \hat{w}_1|$ for increasing values of $d$ in Table 2. Note that mode collapse occurs when $\hat{w}_1 \in \{0, 1\}$. We observe that GMM-LRDS outperforms competing methods in all the considered dimensions. In Appendix I.5.1, we complete these mode weight results with probability metrics (see Appendix B), present further experiments with bi-modal Gaussian mixtures with lower condition numbers and dimension, which show that competing algorithms are able to perform on par with LRDS in these simpler settings, and conduct ablation studies on this two-mode type of target (by lightly perturbing the reference distribution or increasing the distance between the modes). To further assess the performance and robustness of LRDS, we also consider the extension to more than two modes in Appendix I.5.2, where LRDS still outperforms its competitors.

**Field system $\phi^4$ from statistical mechanics.** Next, we sample from the 1D $\phi^4$ model, previously studied (Gabrié et al., 2022; Grenioux et al., 2024). At the chosen temperature, the distribution has two well distinct modes with respective weights $w_-$ and $w_+$ such that the relative weight $w_-/w_+$ can be adjusted through a ‘local-field’ parameter $h$. We discretize this continuous model with a grid size of 32 (i.e., $d = 32$). For each method, we compute a Monte Carlo estimation of $w_-/w_+$ and compare the results with a Laplace approximation (0-th and 2-nd orders), see Appendix H.1 for the computations. In this setting, all competing approaches suffer from mode collapse while GMM-LRDS is close to the ground truth as shown by Figure 4. Complete results showing the failure of the competing methods are provided in Appendix I.5.4. Given the satisfying results obtained with the lightweight GMM-LRDS sampler, we did not run EBM-LRDS for this target distribution.

<footnote n="2">²In Appendix G, we give expressions of additional variational metrics proposed by Blessing et al. (2024) to quantify mode collapse in the specific case of RDS.</footnote>
<table><caption>Table 2: Absolute mode weight estimation error for a bi-modal Gaussian mixture with growing *d*, averaged over 16 sampling runs. Bold font indicates best result, orange cells refer to settings with uninformative mode weight estimation (i.e., uniform mixture), red cells denote mode collapse. N/A denotes settings with numerical issues.</caption><thead><tr><th>Algorithm</th><th>d = 16 ↓</th><th>d = 32 ↓</th><th>d = 64 ↓</th></tr></thead><tbody><tr><td>SMC</td><td>11.4%±9.1%</td><td>15.8%±8.5%</td><td>15.2%±7.5%</td></tr><tr><td>RE</td><td>16.5%±1.3%</td><td>15.9%±1.4%</td><td>17.0%±1.4%</td></tr><tr><td>LV-PIS</td><td>6.0%±3.4%</td><td>33.2%±0.1%</td><td>33.0%±0.1%</td></tr><tr><td>LV-DDS</td><td>11.8%±9.3%</td><td>31.5%±2.9%</td><td>33.1%±0.1%</td></tr><tr><td>LV-DIS</td><td>16.4%±0.5%</td><td>16.5%±0.4%</td><td>16.8%±0.6%</td></tr><tr><td>LV-CMCD</td><td>32.2%±15.4%</td><td>50.1%±8.8%</td><td>16.3%±10.6%</td></tr><tr><td>iDEM</td><td>33.3%±0.0%</td><td>66.7%±0.0%</td><td>11.7%±0.4%</td></tr><tr><td>PDDS</td><td>0.8%±0.6%</td><td>66.7%±0.0%</td><td>N/A</td></tr><tr><td>GMM-LRDS</td><td>1.7%±0.6%</td><td>2.7%±0.8%</td><td>4.1%±0.6%</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: Estimation of the relative weight of $\phi^4$ modes with increasing $h$, averaged over 16 sampling runs.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Samples obtained for Rings distribution. Reasonable results could not be obtained with PDDS due to numerical issues.</figcaption></figure>

**Compactly supported multi-modal distributions.** Then, we aim to sample from 2-dimensional multi-modal distributions with complex geometries. We consider (a) *Rings distribution*, which has 3 ring-shaped modes, and (b) *Checkerboard distribution*, which has 8 square-shaped modes. In both cases, the modes are not evenly weighted. In this setting, we consider $J = 64$ components for GMM-LRDS and leverage the Replica Exchange algorithm as backbone annealed MCMC sampler in the EBM-LRDS training algorithm. Apart from adaptive methods, the structure of Rings modes are well recovered by all approaches, see Figure 5. However, we observe that non diffusion-based approaches fail to recover the ground truth mode weights. On the other hand, LRDS is the only sampling method to obtain samples that are close to exact for the Checkerboard distribution, see Figure 6, while being able to correctly estimate the ground truth mode weights.

<figure><img src="image_7.png" /><figcaption>Figure 6: Samples obtained for Checkerboard distribution. All sampling methods except LRDS fail to provide accurate samples.</figcaption></figure>

**Bayesian logistic regression models.** Lastly, we provide in Appendix I.5.5 sampling results on standard Bayesian posterior distributions, which are not however explicitly multi-modal.

# 6 DISCUSSION

Building on the class of recently developed annealed VI methods, this paper presents the *Learned Reference-based Diffusion Sampler* (LRDS), which specifically addresses the challenging case of multi-modal target distributions. In essence, LRDS aims at learning a reference process – based on GMMs or multi-level EBMs – adapted to the target distribution by using samples obtained via local MCMC samplers initialized in the modes. Our numerical experiments show that GMM-LRDS accurately recovers the global information of the relative weights of the modes in several tens of dimensions, unlike competing methods, and that EBM-LRDS can help tackling distributions with non-Gaussian properties such as sharp supports. However this advantage of LRDS comes at the computational cost of the necessary pre-training of the reference process model. Concerning future work, we note in particular that EBM-LRDS is a promising tool for real-world sampling tasks on non-euclidean spaces which may benefit from the flexible definition of the reference process as an EBM. An interesting test case of this kind that shall be considered is the sampling of the Boltzmann distribution of proteins in internal coordinates.
## ACKNOWLEDGEMENTS

We would like to thank Lorenz Richter, Julius Berner, Denis Blessing, Junhua Chen, Francisco Vargas and Yazid Janati for many stimulating discussions. LG and MG acknowledge funding from Hi! Paris. AD and MN would like to thank the Isaac Newton Institute for Mathematical Sciences and the Alan Turing Institute for support and hospitality during the programme *Diffusions in machine learning: Foundations, generative models and non-convex optimisation* when work on this paper was undertaken. Part of the work of AD is funded by the European Union (ERC, Ocean, 101071601). Views and opinions expressed here are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. This work was partially performed using HPC resources from GENCI-IDRIS (AD011015234).

## REFERENCES

Tara Akhound-Sadegh, Jarrid Rector-Brooks, Joey Bose, Sarthak Mittal, Pablo Lemos, Cheng-Hao Liu, Marcin Sendera, Siamak Ravanbakhsh, Gauthier Gidel, Yoshua Bengio, Nikolay Malkin, and Alexander Tong. Iterated denoising energy matching for sampling from boltzmann densities. In *Proceedings of the 41st International Conference on Machine Learning*, volume 235 of *Proceedings of Machine Learning Research*, pp. 760–786. PMLR, 21–27 Jul 2024. URL https://proceedings.mlr.press/v235/akhound-sadegh24a.html.

Michael Arbel, Alex Matthews, and Arnaud Doucet. Annealed flow transport Monte Carlo. In *International Conference on Machine Learning*, pp. 318–330. PMLR, 2021.

Oleg Arenz, Philipp Dahlinger, Zihan Ye, Michael Volpp, and Gerhard Neumann. A unified perspective on natural gradient variational inference with gaussian mixture models. *Transactions on Machine Learning Research*, 2023. ISSN 2835-8856. URL https://openreview.net/forum?id=tLBjsX4tjs.

Julius Berner, Lorenz Richter, and Karen Ullrich. An optimal control perspective on diffusion-based generative modeling. *Transactions on Machine Learning Research*, 2023. URL https://openreview.net/forum?id=oYIjw37pTP.

Denis Blessing, Xiaogang Jia, Johannes Esslinger, Francisco Vargas, and Gerhard Neumann. Beyond ELBOs: A large-scale evaluation of variational methods for sampling. In *Proceedings of the 41st International Conference on Machine Learning*, volume 235 of *Proceedings of Machine Learning Research*, pp. 4205–4229. PMLR, 21–27 Jul 2024. URL https://proceedings.mlr.press/v235/blessing24a.html.

Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng. *Handbook of Markov Chain Monte Carlo*. Chapman and Hall/CRC, May 2011. ISBN 9780429138508. doi: 10.1201/b10905. URL http://dx.doi.org/10.1201/b10905.

Patrick Cattiaux, Giovanni Conforti, Ivan Gentil, and Christian Léonard. Time reversal of diffusion processes under a finite entropy condition. *Annales de l'Institut Henri Poincaré, Probabilités et Statistiques*, 59(4):1844–1881, 2023. doi: 10.1214/22-AIHP1320. URL https://doi.org/10.1214/22-AIHP1320.

Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. 26, 2013. URL https://proceedings.neurips.cc/paper_files/paper/2013/file/af21d0c97db2e27e13572cbf59eb343d-Paper.pdf.

Yin Dai, Yuan Gao, Jian Huang, Yuling Jiao, Lican Kang, and Jin Liu. Lipschitz transport maps via the Föllmer flow. *arXiv preprint arXiv:2309.03490*, 2023.

Arnak S Dalalyan. Theoretical guarantees for approximate sampling from smooth and log-concave densities. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 79(3):651–676, 2017.

Valentin De Bortoli, Michael Hutchinson, Peter Wirnsberger, and Arnaud Doucet. Target score matching. *arXiv preprint arXiv:2402.08667*, 2024.
Pierre Del Moral, Arnaud Doucet, and Ajay Jasra. Sequential Monte Carlo samplers. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 68(3):411–436, 2006.

Arthur P Dempster, Nan M Laird, and Donald B Rubin. Maximum likelihood from incomplete data via the em algorithm. *Journal of the royal statistical society: series B (methodological)*, 39(1): 1–22, 1977.

Prafulla Dhariwal and Alexander Quinn Nichol. Diffusion models beat GANs on image synthesis. In *Advances in Neural Information Processing Systems*, 2021. URL https://openreview.net/forum?id=AAWuCvzaVt.

Arnaud Doucet, Will Grathwohl, Alexander G Matthews, and Heiko Strathmann. Score-based diffusion meets annealed importance sampling. *Advances in Neural Information Processing Systems*, 35:21482–21494, 2022.

Simon Duane, A.D. Kennedy, Brian J. Pendleton, and Duncan Roweth. Hybrid monte carlo. *Physics Letters B*, 195(2):216–222, 1987. ISSN 0370-2693. doi: https://doi.org/10.1016/0370-2693(87)91197-X. URL https://www.sciencedirect.com/science/article/pii/037026938791197X.

Alain Durmus and Éric Moulines. Quantitative bounds of convergence for geometrically ergodic Markov chain in the Wasserstein distance with application to the Metropolis Adjusted Langevin Algorithm. *Statistics and Computing*, 25:5–19, 2015.

Alain Durmus and Eric Moulines. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. 2017.

Marylou Gabrié, Grant M Rotskoff, and Eric Vanden-Eijnden. Adaptive Monte Carlo augmented with normalizing flows. *Proceedings of the National Academy of Sciences*, 119(10):e2109420119, 2022.

Marylou Gabrié, Grant M. Rotskoff, and Eric Vanden-Eijnden. Adaptive Monte Carlo augmented with normalizing flows. *Proceedings of the National Academy of Sciences*, 119(10): e2109420119, 2022. doi: 10.1073/pnas.2109420119. URL https://www.pnas.org/doi/abs/10.1073/pnas.2109420119.

Ruiqi Gao, Yang Song, Ben Poole, Ying Nian Wu, and Diederik P Kingma. Learning energy-based models by diffusion recovery likelihood. In *International Conference on Learning Representations*, 2021. URL https://openreview.net/forum?id=v.1Soh8QUNC.

Tomas Geffner and Justin Domke. Mcmc variational inference via uncorrected hamiltonian annealing. *Advances in Neural Information Processing Systems*, 34:639–651, 2021.

Tomas Geffner and Justin Domke. Langevin diffusion variational inference. In *International Conference on Artificial Intelligence and Statistics*, pp. 576–593. PMLR, 2023.

Louis Grenioux, Alain Oliviero Durmus, Eric Moulines, and Marylou Gabrié. On sampling with approximate transport maps. In *Proceedings of the 40th International Conference on Machine Learning*, volume 202 of *Proceedings of Machine Learning Research*, pp. 11698–11733. PMLR, 23–29 Jul 2023. URL https://proceedings.mlr.press/v202/grenioux23a.html.

Louis Grenioux, Maxence Noble, Marylou Gabrié, and Alain Oliviero Durmus. Stochastic localization via iterative posterior sampling. In *Proceedings of the 41st International Conference on Machine Learning*, volume 235 of *Proceedings of Machine Learning Research*, pp. 16337–16376. PMLR, 21–27 Jul 2024. URL https://proceedings.mlr.press/v235/grenioux24a.html.

Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. *The Journal of Machine Learning Research*, 13(1):723–773, 2012.

Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. *Advances in neural information processing systems*, 33:6840–6851, 2020.
Matthew D. Hoffman and Andrew Gelman. The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo. *Journal of Machine Learning Research*, 15(47):1593–1623, 2014. URL http://jmlr.org/papers/v15/hoffman14a.html.

Lars Holdijk, Yuanqi Du, Ferry Hooft, Priyank Jaini, Bernd Ensing, and Max Welling. PIPS: Path integral stochastic optimal control for path sampling in molecular dynamics, 2023. URL https://openreview.net/forum?id=TnIZfXSFJAh.

Xunpeng Huang, Hanze Dong, Yifan HAO, Yian Ma, and Tong Zhang. Reverse Diffusion Monte Carlo. In *The Twelfth International Conference on Learning Representations*, 2024a. URL https://openreview.net/forum?id=kIPEyMSdFV.

Xunpeng Huang, Difan Zou, Hanze Dong, Yi-An Ma, and Tong Zhang. Faster sampling without isoperimetry via diffusion-based monte carlo. In *Proceedings of Thirty Seventh Conference on Learning Theory*, volume 247 of *Proceedings of Machine Learning Research*, pp. 2438–2493. PMLR, 30 Jun–03 Jul 2024b. URL https://proceedings.mlr.press/v247/huang24a.html.

Aapo Hyvärinen and Peter Dayan. Estimation of non-normalized statistical models by score matching. *Journal of Machine Learning Research*, 6(4), 2005.

Ghassen Jerfel, Serena Wang, Clara Wong-Fannjiang, Katherine A. Heller, Yian Ma, and Michael I. Jordan. Variational refinement for importance sampling using the forward Kullback-Leibler divergence. In *Proceedings of the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence*, pp. 1819–1829. PMLR, December 2021. URL https://proceedings.mlr.press/v161/jerfel21a.html. ISSN: 2640-3498.

Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion-based generative models. *Advances in Neural Information Processing Systems*, 35:26565–26577, 2022.

Beomsu Kim and Jong Chul Ye. Denoising MCMC for accelerating diffusion-based generative models. In *Proceedings of the 40th International Conference on Machine Learning*, volume 202 of *Proceedings of Machine Learning Research*, pp. 16955–16977. PMLR, 23–29 Jul 2023. URL https://proceedings.mlr.press/v202/kim23z.html.

Werner Krauth. Statistical mechanics: algorithms and computations. *OUP Oxford*, 13, 2006.

Dirk P. Kroese, Thomas Taimre, and Zdravko I. Botev. *Handbook of Monte Carlo methods*. Number 706 in Wiley series in probability and statistics. Wiley, Hoboken, 2011. ISBN 978-0-470-17793-8.

Christian Léonard. Some properties of path measures. *Séminaire de Probabilités XLVI*, pp. 207–230, 2014.

Jun S Liu and Jun S Liu. *Monte Carlo strategies in scientific computing*, volume 10. Springer, 2001.

Alex Matthews, Michael Arbel, Danilo Jimenez Rezende, and Arnaud Doucet. Continual repeated annealed flow transport monte carlo. In *International Conference on Machine Learning*, pp. 15196–15219. PMLR, 2022.

Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller. Equation of State Calculations by Fast Computing Machines. *The Journal of Chemical Physics*, 21(6):1087–1092, 06 1953. ISSN 0021-9606. doi: 10.1063/1.1699114. URL https://doi.org/10.1063/1.1699114.

Laurence Illing Midgley, Vincent Stimper, Gregor N. C. Simm, Bernhard Schölkopf, and José Miguel Hernández-Lobato. Flow annealed importance sampling bootstrap. In *The Eleventh International Conference on Learning Representations*, 2023. URL https://openreview.net/forum?id=XCTVFJwS9LJ.

Radford M Neal. Annealed importance sampling. *Statistics and computing*, 11:125–139, 2001.

Radford M Neal. Mcmc using hamiltonian dynamics. *arXiv preprint arXiv:1206.1901*, 2012.
Alexander Quinn Nichol and Prafulla Dhariwal. Improved denoising diffusion probabilistic models. In *International Conference on Machine Learning*, pp. 8162–8171. PMLR, 2021.

Frank Noé, Simon Olsson, Jonas Köhler, and Hao Wu. Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning. *Science*, 365(6457):eaaw1147, September 2019. ISSN 0036-8075, 1095-9203. doi: 10.1126/science.aaw1147. URL https://www.science.org/doi/10.1126/science.aaw1147.

Nikolas Nüsken and Lorenz Richter. Solving high-dimensional hamilton–jacobi–bellman pdes using neural networks: perspectives from the theory of controlled diffusions and measures on path space. *Partial differential equations and applications*, 2(4):48, 2021.

George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. *Journal of Machine Learning Research*, 22(57):1–64, 2021.

Michele Pavon. On local entropy, stochastic control, and deep neural networks. *IEEE Control Systems Letters*, 7:437–441, 2022.

Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit-learn: Machine learning in python. *Journal of machine learning research*, 12(Oct):2825–2830, 2011.

Gabriel Peyré, Marco Cuturi, et al. Computational optimal transport: With applications to data science. *Foundations and Trends® in Machine Learning*, 11(5-6):355–607, 2019.

Angus Phillips, Hai-Dang Dau, Michael John Hutchinson, Valentin De Bortoli, George Deligiannidis, and Arnaud Doucet. Particle denoising diffusion sampler. In *Proceedings of the 41st International Conference on Machine Learning*, volume 235 of *Proceedings of Machine Learning Research*, pp. 40688–40724. PMLR, 21–27 Jul 2024. URL https://proceedings.mlr.press/v235/phillips24a.html.

Emilia Pompe, Chris Holmes, and Krzysztof Łatuszyński. A framework for adaptive MCMC targeting multimodal distributions. *The Annals of Statistics*, 48(5):2930–2952, October 2020. ISSN 0090-5364, 2168-8966. doi: 10.1214/19-AOS1916. URL https://projecteuclid.org/journals/annals-of-statistics/volume-48/issue-5/A-framework-for-adaptive-MCMC-targeting-multimodal-distributions/10.1214/19-AOS1916.full. Publisher: Institute of Mathematical Statistics.

Teodora Reu, Francisco Vargas, Anna Kerekes, and Michael M. Bronstein. To smooth a cloud or to pin it down: Expressiveness guarantees and insights on score matching in denoising diffusion models. In *The 40th Conference on Uncertainty in Artificial Intelligence*, 2024. URL https://openreview.net/forum?id=BV2STKHKE6.

D J Rezende, S Mohamed, and D Wierstra. Stochastic backpropagation and approximate inference in deep generative models. *Proceedings of The 31st International Conference in Machine Learning*, Beijing China, 32:1278–, 2014. URL http://jmlr.org/proceedings/papers/v32/rezende14.html%5Cnpapers3://publication/uuid/F2747569-7719-4EAC-A5A7-9ECA9D6A8FE6. arXiv:1401.4082v3 ISBN: 9781634393973.

Lorenz Richter, Julius Berner, and Guan-Horng Liu. Improved sampling via learned diffusions. In *ICML Workshop on New Frontiers in Learning, Control, and Dynamical Systems*, 2023. URL https://openreview.net/forum?id=uLgYD7ie0O.

Gareth O. Roberts and Richard L. Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. *Bernoulli*, 2(4):341–363, 1996. ISSN 13507265. URL http://www.jstor.org/stable/3318418.

Tim Salimans and Jonathan Ho. Should EBMs model the energy or the score?, 2021. URL https://openreview.net/forum?id=9AS-TF2jRNb.
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In *International conference on machine learning*, pp. 2256–2265. PMLR, 2015.

Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. In *Advances in Neural Information Processing Systems*, volume 32. Curran Associates, Inc., 2019. URL https://proceedings.neurips.cc/paper_files/paper/2019/file/3001ef257407d5a371a96dcd947c7d93-Paper.pdf.

Yang Song and Stefano Ermon. Improved techniques for training score-based generative models. In *Advances in Neural Information Processing Systems*, volume 33, pp. 12438–12448. Curran Associates, Inc., 2020. URL https://proceedings.neurips.cc/paper_files/paper/2020/file/92c3b916311a5517d9290576e3ea37ad-Paper.pdf.

Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. In *The Ninth International Conference on Learning Representations*, 2021. URL https://openreview.net/forum?id=PdTIG12RRHS.

Gabriel Stoltz, Mathias Rousset, et al. *Free energy computations: A mathematical perspective.* World Scientific, 2010.

Robert H. Swendsen and Jian-Sheng Wang. Replica Monte Carlo Simulation of Spin-Glasses. *Physical Review Letters*, 57(21):2607-2609, November 1986. doi: 10.1103/PhysRevLett.57.2607. URL https://link.aps.org/doi/10.1103/PhysRevLett.57.2607. Publisher: American Physical Society.

Saifuddin Syed, Alexandre Bouchard-Côté, George Deligiannidis, and Arnaud Doucet. Non-reversible parallel tempering: a scalable highly parallel mcmc scheme. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 84(2):321–350, 2022.

Nicholas G Tawn, Gareth O Roberts, and Jeffrey S Rosenthal. Weight-preserving simulated tempering. *Statistics and Computing*, 30(1):27–41, 2020.

Belinda Tzen and Maxim Raginsky. Theoretical guarantees for sampling and inference in generative models with latent diffusions. In *Conference on Learning Theory*, pp. 3084–3114. PMLR, 2019.

Francisco Vargas, Will Sussman Grathwohl, and Arnaud Doucet. Denoising diffusion samplers. In *The Eleventh International Conference on Learning Representations*, 2023a. URL https://openreview.net/forum?id=8pvnfTABulf.

Francisco Vargas, Andrius Ovsianas, David Fernandes, Mark Girolami, Neil D Lawrence, and Nikolas Nüsken. Bayesian learning via neural Schrödinger–Föllmer flows. *Statistics and Computing*, 33(1):3, 2023b.

Francisco Vargas, Shreyas Padhy, Denis Blessing, and Nikolas Nüsken. Transport meets variational inference: Controlled monte carlo diffusions. In *The Twelfth International Conference on Learning Representations*, 2024. URL https://openreview.net/forum?id=PP1rudnxiW.

Pascal Vincent. A connection between score matching and denoising autoencoders. *Neural computation*, 23(7):1661–1674, 2011.

Martin J. Wainwright and Michael I. Jordan. Graphical Models, Exponential Families, and Variational Inference. *Foundations and Trends® in Machine Learning*, 1:1–305, 2008. doi: 10.1561/2200000001. URL http://discovery.ucl.ac.uk/185880/.

Dawn Woodard, Scott Schmidler, and Mark Huber. Sufficient conditions for torpid mixing of parallel and simulated tempering. 2009.

Hao Wu, Jonas Köhler, and Frank Noé. Stochastic normalizing flows. *Advances in Neural Information Processing Systems*, 33:5933–5944, 2020.

Qinsheng Zhang and Yongxin Chen. Path Integral Sampler: a stochastic control approach for sampling. In *The Tenth International Conference on Learning Representations*, 2022.
Xinwei Zhang, Zhiqiang Tan, and Zhijian Ou. Persistently trained, diffusion-assisted energy-based models. *Stat*, 12(1):e625, 2023. doi: https://doi.org/10.1002/sta4.625. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sta4.625.

Yaxuan Zhu, Jianwen Xie, Ying Nian Wu, and Ruiqi Gao. Learning energy-based models by co-operative diffusion recovery likelihood. In *The Twelfth International Conference on Learning Representations*, 2024. URL https://openreview.net/forum?id=AyzkDpuqcl.
# ORGANIZATION OF THE SUPPLEMENTARY

The appendix is organized as follows. In Appendix A, we provide the pseudo-codes of every RDS instance presented in the main part of this paper (standard RDS, GMM-LRDS & EBM-LRDS). Appendix B summarizes general facts that will be useful for proofs and corresponding computations. In Appendix C, we describe the noising diffusion schemes considered in this work, namely the Pinned Brownian Motion and the Variance-Preserving setting, and propose novel schemes. For RDS and related variational diffusion-based methods, we dispense details on their continuous-time and corresponding discrete-time variational objectives in Appendix D. We respectively provide an overview of existing annealed MCMC methods and multi-level Energy-Based models in Appendix E and Appendix F. In Appendix G, we provide expressions of additional variational metrics proposed by Blessing et al. (2024) in the RDS setting. Details on our experimental settings, our implementation of RDS and other competing methods are given in Appendix H. Finally, we present additional theoretical and experimental results in Appendix I.

Our codebase is available at https://github.com/h2o64/sde_sampler_lrds.

**Notation.** Let $\mathbb{P} \in \mathcal{P}(C_T)$ and $t \in [0, T]$. For any $x_t \in \mathbb{R}^d$, we denote by $\mathbb{P}_{|t}^{x_t} \in \mathcal{P}(C_T)$ the path measure associated to the stochastic process $(X_t)_{t \in [0,T]} \sim \mathbb{P}$ conditioned on $X_t = x_t$. Furthermore, for any $\pi \in \mathcal{P}(\mathbb{R}^d)$, $\pi \otimes \mathbb{P}_{|t}^{x_t} \in C_T$ stands for the path measure $\int_{\mathbb{R}^d} \mathbb{P}_{|t}^{x_t} d\pi(x_t)$. For any $K \ge 1$, we denote the set of joint distributions $\mathcal{P}((\mathbb{R}^d)^{K+1})$ by $\mathcal{P}^{(K+1)}$. Finally, we adopt the following notation to design sample batches: for any $\bar{L} \ge \underline{L} \ge 0$ and $\bar{K} \ge \underline{K} \ge 0$, we denote by $X_{\underline{K}:\bar{K}}^{\underline{L}:\bar{L}}$ the sequence of samples $\{X_k^\ell\}_{k=\underline{K},\ell=\underline{L}}^{\bar{K},\bar{L}}$.

## A PSEUDO-CODES OF RDS-BASED ALGORITHMS

We respectively give sampling procedures and training procedures of a general version of RDS in Algorithm 1 and Algorithm 2. Relying on this, we derive the complete training schemes for GMM-LRDS (Algorithm 3) and EBM-LRDS (Algorithm 4).

**Algorithm 1: Reference-based Diffusion Sampling (RDS) : sampling stage**

**Input:** Time discretization $\{t_k\}_{k=0}^K$ of $[0, T]$, reference scores $\{s_{T-t_k}^{\text{ref}}\}_{k=0}^{K-1}$, neural network $g^\theta: [0, T] \times \mathbb{R}^d \to \mathbb{R}^d$, variational coefficients $\{a_k, b_k, c_k\}_{k=0}^{K-1}$ defined in Appendix D.2

▷ Initialization

$$Y_0 \sim \pi^{\text{base}}$$

$$(Z_k)_{k=0}^{K-1} \stackrel{\text{i.i.d.}}{\sim} N(0, I_d)$$

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

▷ Compute the $k$-th diffusion step

$$Y_{k+1} = a_k Y_k + b_k \{s_{T-t_k}^{\text{ref}}(Y_k) + g_{T-t_k}^{\hat{\theta}}(Y_k)\} + \sqrt{c_k} Z_k$$

**Output:** Discrete time process $Y_{0:K}$ approximating $(Y_{t_k}^{\hat{\theta}})_{k=0}^K$

**Algorithm 2: Reference-based Diffusion Sampling (RDS) : training stage**

**Input:** Time discretization $\{t_k\}_{k=0}^K$ of $[0, T]$, target density $\gamma$, reference density $\gamma^{\text{ref}}$ and scores $\{s_{T-t_k}^{\text{ref}}\}_{k=0}^{K-1}$, neural network $g^\theta: [0, T] \times \mathbb{R}^d \to \mathbb{R}^d$ initialized such that $g^{\theta_0} = 0$, number of training iterations $N$, batch size $B$, variational coefficients $\{w_k, a_k, b_k, c_k\}_{k=0}^{K-1}$ defined in Appendix D.2

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

▷ Simulate $B$ trajectories of the process $(Y_t^{\hat{\theta}})_{t \in [0,T]}$

$$Y_{0:K}^{1:B} \stackrel{\text{i.i.d.}}{\sim} \text{SamplingRDS}(g^{\theta_n}), \text{see Algorithm 1}$$

▷ Apply a stochastic gradient descent step on $\theta_n$

Compute a MC estimator $\hat{L}_{\text{RDS}}(\theta_n)$ of the loss $\mathcal{L}_{\text{RDS}}(\theta_n)$ defined in (7) with samples $Y_{0:K}^{1:B}$

Compute the gradient $\nabla_{\theta} \hat{L}_{\text{RDS}}(\theta_n)$ and update $\theta_n$ to $\theta_{n+1}$ with Adam optimizer

**Output:** Learned guidance $g^{\theta_N}$
**Algorithm 3: Gaussian Mixture Model LRDS (GMM-LRDS) : training stage**

**Input:** Time discretization {$t_k$}$_{k=0}^K$ of $[0, T]$, target density $\gamma$, location of the target modes {$\mathbf{x}_i$}$_{i=1}^I$, number of Markov chains per target mode $M \ge 1$, size of Markov chains $N_{tot} \ge 1$, effective size of Markov chains $N_{eff} \in [1, N_{tot}]$, number of reference modes $J \ge I$

▷ (a) Obtain reference samples

For each $i \in \{1, \dots, I\}$, build via MALA M Markov chains {$X_{0:N_{tot}}^{i,m}$}$_{m=1}^M$ of size $N_{tot}$ starting at $\mathbf{x}_i$

Keep the last $N_{eff}$ samples of each Markov chain to define $\hat{\pi}^{\text{ref}} \approx \{X_{N_{tot}-N_{eff}+1:N_{tot}}^{i,m}\}_{i=1,m=1}^{I,M}$

▷ (b) Define the reference process

if $J = 1$ then

Fit on $\hat{\pi}^{\text{ref}}$ a Gaussian model with parameterized density $\gamma^{\varphi}$ (Maximum Likelihood estimation)

else

Fit on $\hat{\pi}^{\text{ref}}$ a J-component Gaussian Mixture Model with parameterized density $\gamma^{\varphi}$ (EM algorithm)
Set $\gamma^{\text{ref}} = \gamma^{\varphi}$ and compute $\{s_{T-t_k}^{\text{ref}}\}_{k=0}^{K-1}$ analytically, see Appendix B

▷ (c) Run the RDS variational optimization

Run Algorithm 2

**Output:** Trained GMM-LRDS sampler

**Algorithm 4: Energy-Based Model LRDS (EBM-LRDS) : training stage**

**Input:** Time discretization {$t_k$}$_{k=0}^K$ of $[0, T]$, target density $\gamma$, location of the target modes {$\mathbf{x}_i$}$_{i=1}^I$, number of Markov chains per target mode $M \ge 1$, size of Markov chains $N_{tot} \ge 1$, effective size of Markov chains $N_{eff} \in [1, N_{tot}]$

▷ (a) Obtain reference samples

For each $i \in \{1, \dots, I\}$, build via MALA M Markov chains {$X_{0:N_{tot}}^{i,m}$}$_{m=1}^M$ of size $N_{tot}$ starting at $\mathbf{x}_i$

Keep the last $N_{eff}$ samples of each Markov chain to define $\hat{\pi}^{\text{ref}} \approx \{X_{N_{tot}-N_{eff}+1:N_{tot}}^{i,m}\}_{i=1,m=1}^{I,M}$

▷ (b) Define the reference process

Based on $\hat{\pi}^{\text{ref}}$, fit a multi-level EBM $E^{\varphi}$ using Algorithm 11

Set $\gamma^{\text{ref}}(x) = \exp(-E^{\varphi}(0, x))$ and $s_{T-t_k}^{\text{ref}}(x) = -\nabla_x E^{\varphi}(T-t_k, x)$ for any $k \in \{0, \dots, K-1\}$

▷ (c) Run the RDS variational optimization

Run Algorithm 2

**Output:** Trained EBM-RDS sampler

B PRELIMINARIES

**Linear SDE integration.** We first dispense a useful lemma to compute exact integration in SDEs with linear drift.

**A0 (Integrability conditions on f and β).** Coefficients $f: [0, T] \to \mathbb{R}$ and $\beta: [0, T] \to (0, \infty)$ are such that (a) $f$ is integrable on $(0, T)$ and (b) $\beta$ is integrable on $(0, T)$.

**Lemma 2.** Let $T > 0$ and $b \in \mathbb{R}^d$. Consider the SDE defined on $[0, T]$ by $dY_t = f(t)(Y_t + b)dt + \sqrt{\beta(t)}dB_t$, where coefficients $f$ and $\beta$ verify **A0**. Then, for any pair of time-steps $(s, t)$ such that $T \ge t > s \ge 0$, the conditional distribution of $Y_t$ given $Y_s = y_s \in \mathbb{R}^d$, denoted by $p_{t|s}(\cdot|y_s)$, verifies

$$p_{t|s}(\cdot|y_s) = N \left( \exp\left(\int_s^t f(u)du\right)y_s + \left(\exp\left(\int_s^t f(u)du\right) - 1\right)b, \int_s^t \beta(u) \exp\left(2\int_u^t f(r)dr\right)du I_d \right).$$

*Proof.* Assume **A0**. Define the function $\zeta: t \in [0, T] \to \exp(-\int_0^t f(u)du)$ and consider the stochastic process $(Z_t)_{t \in [0, T]}$ defined by $Z_t = \zeta(t)Y_t$ for any $t \in [0, T]$. By $\hat{\text{Ito}}$'s formula, we have $dZ_t = f(t)\zeta(t)bdt + \zeta(t)\sqrt{\beta(t)}dB_t = -\dot{\zeta}(t)bdt + \zeta(t)\sqrt{\beta(t)}dB_t$. Therefore, for any time-steps $(s, t)$ such that $T \ge t > s \ge 0$, we have

$$\zeta(t)Y_t - \zeta(s)Y_s = \{\zeta(s) - \zeta(t)\}b + \int_s^t \zeta(u)\sqrt{\beta(u)}dB_u,$$

and then

$$Y_t = \exp\left(\int_s^t f(u)du\right)Y_s + \left(\exp\left(\int_s^t f(u)du\right) - 1\right)b + \int_s^t \sqrt{\beta(u)} \exp\left(\int_u^t f(r)dr\right)dB_u,$$

which gives the result using $\hat{\text{Ito}}$'s isometry and the fact that $Y_s$ is independent from $(B_t - B_s)_{t \in [s, T]}$.

<figure><img src="image_8.png" /></figure>
**General SDE integration.** Consider a general SDE of the form $dY_t = f(t)(Y_t + g_t(Y_t))dt + \sqrt{\beta(t)}dB_t$ defined on $[0, T]$, where $f, \beta$ satisfy A0 and $g: [0, T] \times \mathbb{R}^d \to \mathbb{R}^d$ is a black-box function. Let $(s, t)$ be a pair of time-steps such that $T \ge t > s \ge 0$. We aim to compute the conditional distribution of $Y_t$ given $Y_s = y_s \in \mathbb{R}^d$. While this distribution is intractable for general function $g$, it can be approximated by a tractable Gaussian kernel if $t$ is close enough to $s$. Below, we restate two common first-order methods that produce such approximation following the result of Lemma 2.

(a) Euler-Maruyama (EM) scheme, which consists of exactly integrating on time interval $[s, t]$ the SDE defined by $dY_u = f(s)(y_s+g_s(y_s))du + \sqrt{\beta(s)}dB_u$ with constant drift and constant volatility.

(b) Exponential Integration (EI) scheme (Durmus & Moulines, 2015), more precise than the EM scheme, which consists of exactly integrating on time interval $[s, t]$ the SDE defined by $dY_u = f(u)(Y_u+g_u(y_u))du + \sqrt{\beta(u)}dB_u$ with linear drift, where the intractable term is frozen at $g_s(y_s)$.

**Computation of scores for Gaussian mixtures.** Below, we provide a useful lemma to compute the marginal scores of a noising diffusion process applied to a general Gaussian mixture.

**Lemma 3.** Let $\pi \in \mathcal{P}(\mathbb{R}^d)$. Consider the noising process $(X_t)_{t \in [0,T]}$ defined by SDE (2) where $X_0 \sim \pi$. Then, for any $t \in (0, T)$, the conditional distribution of $X_t$ given $X_0 = x_0 \in \mathbb{R}^d$ is defined by the Gaussian kernel

$$q_{t|0}(\cdot|x_0) = N(S(t)x_0, S(t)^2\sigma^2(t)I_d), \quad (9)$$

where $S(t) = \exp(\int_0^t f(u)du)$ and $\sigma^2(t) = \int_0^t \beta(u)/S(u)^2du$.

In particular, if $\pi$ is a Gaussian distribution i.e., $\pi = N(\mathbf{m}, \Sigma)$ with mean $\mathbf{m} \in \mathbb{R}^d$ and covariance $\Sigma \in \mathbb{R}^{d \times d}$, then for any $t \in [0, T]$, the marginal distribution of $X_t$ is Gaussian, defined by

$$p_t = N(S(t)\mathbf{m}, S(t)^2\Sigma + S(t)^2\sigma^2(t)I_d). \quad (10)$$

On the other hand, if $\pi$ is a Gaussian mixture i.e., $\pi = \sum_{j=1}^J w_j N(\mathbf{m}_j, \Sigma_j)$ where $J \ge 1$ is the number of components, $\{w_j\}_{j=1}^J \in \Delta_J$ are the component weights, $\mathbf{m}_j \in \mathbb{R}^d$ and $\Sigma_j \in \mathbb{R}^{d \times d}$ are respectively the mean and the covariance of the j-th component, for any $t \in [0, T]$, the marginal distribution of $X_t$ is a Gaussian mixture, defined by

$$p_t = \sum_{j=1}^{J} w_j N(S(t)\mathbf{m}_j, S(t)^2\Sigma_j + S(t)^2\sigma^2(t)I_d). \quad (11)$$

Note that differentiating (11) provides closed-form solutions of the score of the stochastic process.

*Proof.* The proof of (9) directly follows from Lemma 2. Then, (10) and (11) are obtained by applying the affine mapping defined by the transition kernel (9) to the target distribution $\pi$. $\square$

In particular, (11) demonstrates that any noising diffusion scheme applied to a Gaussian mixture does not show any mode switching.

**Integral probability Metrics.** Finally, we dispense the definition of standard probability metrics, that will be used to assess the performance of the considered samplers, when having access to ground truth samples. Consider two distributions $\mu \in \mathcal{P}(\mathbb{R}^d)$ and $\nu \in \mathcal{P}(\mathbb{R}^d)$ that we aim to compare.

We recall that the 2-Wasserstein distance between $\mu$ and $\nu$ is given by

$$W_2(\mu, \nu) = \inf\{\int_{\mathbb{R}^d \times \mathbb{R}^d} \|x_1 - x_0\|^2 d\pi(x_0, x_1) : \pi \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d), \pi_0 = \mu, \pi_1 = \nu\}^{1/2},$$

where $\pi_i$ denotes the i-th marginal of $\pi$ for $i \in \{0, 1\}$. In practice, we rather turn to its entropy regularized version (Peyré et al., 2019)

$$W_{2,\varepsilon}(\mu, \nu) = \inf\{\int_{\mathbb{R}^d \times \mathbb{R}^d} \|x_1 - x_0\|^2 d\pi(x_0, x_1) - \mathcal{H}(\pi) : \pi \in \mathcal{P}(\mathbb{R}^d \times \mathbb{R}^d), \pi_0 = \mu, \pi_1 = \nu\}^{1/2},$$

where $\varepsilon > 0$ is a regularization hyper-parameter and $\mathcal{H}(\pi) = -\int_{\mathbb{R}^d \times \mathbb{R}^d} \log \pi(x_0, x_1)d\pi(x_0, x_1)$ refers to as the entropy of $\pi$. When having access to samples from $\mu$ and $\nu$, we compute a statistical
estimation of $W_{2,ε}(μ, ν)$ via Sinkhorn algorithm (Cuturi, 2013) with $ε = 10^{-3}$, based on the GPU-oriented implementation from https://github.com/fwilliams/scalable-pytorch-sinkhorn.

We also consider the Maximum Mean Discrepancy (MMD) (Gretton et al., 2012), which quantifies the dissimilarity between $μ$ and $ν$ by comparing their mean embeddings in a reproducing kernel Hilbert space $H_k$ with kernel $k$, as followed

$$MMD(μ, ν) = \sup_{f \in H_K : \|f\|_{H_K} \le 1} \mathbb{E}[f(X) - f(Y)], \quad X \sim μ, Y \sim ν.$$

In particular, $MMD(μ, ν) = 0$ if and only if $μ = ν$ almost surely. When having access to samples $\{x_i\}_{i=1}^n$ from $μ$ and $\{y_i\}_{i=1}^m$ from $ν$, we compute an unbiased statistical estimation of $MMD(μ, ν)$ given by

$$\sqrt{\frac{2}{n(n-1)} \sum_{1 \le i < j \le n} k(x_i, x_j) + \frac{2}{m(m-1)} \sum_{1 \le i < j \le m} k(y_i, y_j) - \frac{2}{mn} \sum_{i=1}^{n} \sum_{j=1}^{m} k(x_i, y_j)},$$

where $k(x, y) = \exp(-\|x - y\|^2 / (2α))$ is the squared exponential kernel whose positive bandwidth is fixed as the median of the squared pairwise distances computed on the joint set $\{x_i\}_{i=1}^n \cup \{y_j\}_{j=1}^m$ (Gretton et al., 2012).

Finally, we recall that the Kolmogorov-Smirnov (KS) distance between $μ$ and $ν$ is defined by

$$KS(μ, ν) = \sup_{x \in \mathbb{R}^d} |F_μ(x) - F_ν(x)|,$$

where $F_μ$ (respectively $F_ν$) denotes the cumulative distribution function of $μ$ (respectively $ν$). When having access to samples from $μ$ and $ν$, we implement a statistical estimation of the sliced KS distance (Grenioux et al., 2023) using 128 random projections.

## C DETAILS ON NOISING DIFFUSION PROCESSES

In this section, we consider a noising diffusion process $(X_t)_{t \in [0, T]}$ induced by SDE (1) for an arbitrary target distribution $π \in \mathcal{P}(\mathbb{R}^d)$ and $T > 0$, and denote by $\mathbb{P}^* \in \mathcal{P}(C_T)$ the corresponding path measure. Below, we provide detailed results on this diffusion process, depending on the choice of $f$ and $β$: we first provide details on the general setting in Appendix C.1; then, we fully describe the Pinned Brownian Motion and the Variance-Preserving scheme in Appendix C.2 and Appendix C.3, that are used in our experiments; finally, we present novel noising schemes in Appendix C.4.

### C.1 GENERAL NOISING SCHEME

Here, we consider the most general form of SDE (1).

**Lemma 4.** Let $π \in \mathcal{P}(\mathbb{R}^d)$. Assume that $f$ and $β$ both verify A0. Then, for any pair of time-steps $(s, t)$ such that $T \ge t > s \ge 0$, the conditional distribution of $X_t$ given $X_s = x_s \in \mathbb{R}^d$ is defined by the Gaussian kernel

$$q_{t|s}(\cdot|x_s) = N(\exp(\log S(t) - \log S(s))x_s, S(t)^2\{\sigma^2(t) - \sigma^2(s)\}I_d), \quad (12)$$

where $S(t) = \exp(\int_0^t f(u)du)$ and $\sigma^2(t) = \int_0^t \beta(u)/S(u)^2du$.

*Proof.* The proof of (12) directly follows from Lemma 2. □

### C.2 PINNED BROWNIAN MOTION

For any $t \in [0, T]$, define $α(t) = \int_0^t \beta(t)dt$ and consider the case where $f(t) = -\frac{β(t)}{α(T)-α(t)}$. Then, SDE (1) can be rewritten as

$$dX_t = -\frac{β(t)X_t}{α(T) - α(t)}dt + \sqrt{β(t)}dW_t, \quad X_0 \sim \pi.$$
This noising scheme, known as the *Pinned Brownian Motion* (PBM) and previously considered by Tzen & Raginsky (2019); Zhang & Chen (2022); Vargas et al. (2023b), can be obtained by applying a Doob's h-transform on the uncontrolled SDE $dX_t = \sqrt{\beta(t)}dW_t$, associated to path measure $\mathbb{Q}$, to hit 0 at time $T$. Interestingly, it can be shown that the resulting path measure $(\mathbb{P}^*)^R$ is solution to the following stochastic optimal control problem

$$ \arg \min \{ \mathrm{KL}(\mathbb{P} | \mathbb{Q}) : \mathbb{P} \in \mathcal{P}(\mathcal{C}_T), \mathbb{P}_T = \pi \}, $$

which is often referred to as a half Schrödinger Bridge problem.

**Lemma 5.** Let $\pi \in \mathcal{P}(\mathbb{R}^d)$. Assume that $f(t) = -\frac{\beta(t)}{\alpha(T)-\alpha(t)}$ and $\beta$ verifies A0. Then, for any pair of time-steps $(s, t)$ such that $T \ge t > s \ge 0$, the conditional distribution of $X_t$ given $X_s = x_s \in \mathbb{R}^d$ is defined by the Gaussian kernel

$$ q_{t|s}(\cdot|x_s) = N\left(\frac{\alpha(T) - \alpha(t)}{\alpha(T) - \alpha(s)}x_s, \frac{(\alpha(T) - \alpha(t))(\alpha(t) - \alpha(s))}{\alpha(T) - \alpha(s)}I_d\right). \quad (13) $$

Since $p_T^*(x) = \int_{\mathbb{R}^d} q_{T|0}(x|x_0)d\pi(x_0)$, it results that $\mathbb{P}_T^* = \delta_0$.

*Proof.* The proof of (13) directly follows from Lemma 2. □

Following Lemma 5, the PBM is an 'exact' noising scheme and we have $\pi^{\text{base}} = \delta_0$ in this setting. Moreover, under mild assumptions on $\pi$, the time-reversed SDE (2) writes as

$$ dY_t = \beta(T-t) \left\{ \frac{Y_t}{\alpha(T) - \alpha(T-t)} + \nabla \log p_{T-t}^*(Y_t) \right\} dt + \sqrt{\beta(T-t)} dB_t, \quad Y_0 = 0. $$

**On the choice of the $\beta$-schedule.** Previous works have considered constant schedule $\beta(t) = \sigma^2$ (and therefore, $\alpha(t) = \sigma^2 t$), where $\sigma > 0$ can be arbitrarily chosen, see e.g., Zhang & Chen (2022); Richter et al. (2023). We also follow this setting in practice.

### C.3 VARIANCE-PRESERVING DIFFUSION

For any $t \in [0, T]$, denote $\alpha(t) = \int_0^t \beta(t)dt$ and consider the case where $f(t) = -\beta(t)/2\sigma^2$ with $\sigma > 0$ and $\beta$ being such that $\int_0^T \beta(s)ds \gg 1$. Then, SDE (1) can be rewritten as

$$ dX_t = -\frac{\beta(t)X_t}{2\sigma^2}dt + \sqrt{\beta(t)}dW_t, \quad X_0 \sim \pi. $$

This noising scheme, known as the *Variance-Preserving* (VP) scheme (Song et al., 2021) and previously considered by Vargas et al. (2023a), is an Ornstein-Uhlenbeck diffusion process and is largely used in score-based generative models.

**Lemma 6.** Let $\pi \in \mathcal{P}(\mathbb{R}^d)$ and $\sigma > 0$. Assume that $f(t) = -\beta(t)/2\sigma^2$ and $\beta$ verifies A0 such that $\int_0^T \beta(s)ds \gg 1$. Then, for any pair of time-steps $(s, t)$ such that $T \ge t > s \ge 0$, the conditional distribution of $X_t$ given $X_s = x_s \in \mathbb{R}^d$ is defined by the Gaussian kernel

$$ q_{t|s}(\cdot|x_s) = N\left(\sqrt{1-\lambda_{s,t}}x_s, \sigma^2\lambda_{s,t}I_d\right), \quad (14) $$

where $\lambda_{s,t} = 1 - \exp(\alpha(s) - \alpha(t))$. Since $p_T^*(x) = \int_{\mathbb{R}^d} q_{T|0}(x|x_0)d\pi(x_0)$, it results that $\mathbb{P}_T^* \approx N(0, \sigma^2 I_d)$.

*Proof.* The proof of (14) directly follows from Lemma 2, along with fact that $\lambda_{0,T} \approx 1$. □

Following Lemma 6, the VP scheme is an 'ergodic' noising scheme, converging exponentially fast to the Gaussian distribution $N(0, \sigma^2 I_d)$; therefore, we have $\pi^{\text{base}} = N(0, \sigma^2 I_d)$ in this setting. Moreover, under mild assumptions on $\pi$, the time-reversed SDE (2) writes as

$$ dY_t = \frac{\beta(T-t)}{2} \{Y_t + 2\sigma^2 \nabla \log p_{T-t}^*(Y_t)\} dt + \sigma \sqrt{\beta(T-t)} dB_t, \quad Y_0 \sim \mathbb{P}_T^*, $$

after simple linear time reparameterization.
**On the choice of the β-schedule.** Previous works have considered a linear schedule $\beta(t) = \beta_{\min}(1 - t/T) + \beta_{\max}(t/T)$ where $\beta_{\min} = 0.1$, $\beta_{\max} \in \{10, 20\}$ and $T = 1$, see e.g., Song et al. (2021); Richter et al. (2023); Reu et al. (2024) or cosine parameterization (Nichol & Dhariwal, 2021; Vargas et al., 2023a), which has been proved to perform better in generative modeling. In our sampling experiments, we did not observe any significant difference between these two settings. Hence, we fix the linear $\beta$-schedule to be the default setting for our numerics (except DDS, originally implemented with the cosine schedule), and let $\sigma$ be arbitrarily chosen.

### C.4 ADDITIONAL EXAMPLES OF NOISING DIFFUSION PROCESSES

Below, we present original noising schemes, but did not consider them in our experiments.

**Pinned Ornstein-Uhlenbeck Motion.** We first propose to consider the case where $f(t) = -\frac{\beta(t)}{2} \coth(\int_t^T \frac{\beta(u)}{2} du)$. Then, SDE (1) can be rewritten as

$$dX_t = -\frac{\beta(t)}{2} \coth\left(\int_t^T \frac{\beta(u)}{2} du\right) X_t dt + \sqrt{\beta(t)} dW_t, X_0 \sim \pi.$$

Analogously to the PBM, one can show that this noising scheme can be obtained by applying a Doob's h-transform on the Ornstein-Uhlenbeck process induced by the SDE $d\tilde{X}_t = -\frac{\beta(t)}{2}\tilde{X}_t dt + \sqrt{\beta(t)}dW_t$ to hit 0 at time $T$. We refer to this noising diffusion scheme as the *Pinned Ornstein-Uhlenbeck Motion* (POUM).

**Lemma 7.** Let $\pi \in \mathcal{P}(\mathbb{R}^d)$. Assume that $f(t) = -\frac{\beta(t)}{2} \coth(\int_t^T \frac{\beta(u)}{2} du)$ and $\beta$ verifies A0. Then, for any pair of time-steps $(s, t)$ such that $T \ge t > s \ge 0$, the conditional distribution of $X_t$ given $X_s = x_s \in \mathbb{R}^d$ is defined by the Gaussian kernel $q_{t|s}(\cdot|x_s)$ given by

$$N\left(\frac{\sinh\left(\int_t^T \frac{\beta(u)}{2} du\right)}{\sinh\left(\int_s^T \frac{\beta(u)}{2} du\right)} x_s, 2 \sinh\left(\int_t^T \frac{\beta(u)}{2} du\right)^2 \left\{\coth\left(\int_s^T \frac{\beta(u)}{2} du\right) - \coth\left(\int_t^T \frac{\beta(u)}{2} du\right)\right\} I_d\right). \quad (15)$$

Since $p_T^*(x) = \int_{\mathbb{R}^d} q_{T|0}(x|x_0)d\pi(x_0)$, it results that $P_T^* = \delta_0$.

*Proof.* The proof of (15) directly follows from Lemma 2. □

Following Lemma 7, the POUM scheme is an 'exact' noising scheme and we have $\pi^{\text{base}} = \delta_0$ in this setting. Moreover, under mild assumptions on $\pi$, the time-reversed SDE (2) writes as

$$dY_t = \frac{\beta(T-t)}{2} \{\coth\left(\int_{T-t}^T \frac{\beta(u)}{2} du\right) Y_t + 2\nabla \log p_{T-t}^*(Y_t)\} dt + \sqrt{\beta(T-t)} dB_t, Y_0 = 0.$$

**Gaussianized Pinned Brownian Motion.** Let $\zeta: [0, T] \to (0, \infty)$ be a square integrable function on $(0, T)$. Denote $\alpha(t) = \int_0^t \zeta(u) du$. Inspired by Dai et al. (2023), we propose to consider $f(t) = -\frac{\zeta(t)}{\alpha(T)-\alpha(t)}$ and $\beta(t) = \frac{2\zeta(t)}{\alpha(T)-\alpha(t)}$. Then, SDE (1) can be rewritten as

$$dX_t = -\frac{\zeta(t)X_t}{\alpha(T) - \alpha(t)} dt + \sqrt{\frac{2\zeta(t)}{\alpha(T) - \alpha(t)}} dW_t, X_0 \sim \pi.$$

**Lemma 8.** Let $\pi \in \mathcal{P}(\mathbb{R}^d)$. Assume that $f(t) = -\frac{\zeta(t)}{\alpha(T)-\alpha(t)}$ and $\beta(t) = \frac{2\zeta(t)}{\alpha(T)-\alpha(t)}$, where $\zeta: [0, T] \to (0, \infty)$ is a square integrable function on $(0, T)$. Then, for any pair of time-steps $(s, t)$ such that $T \ge t > s \ge 0$, the conditional distribution of $X_t$ given $X_s = x_s \in \mathbb{R}^d$ is defined by the Gaussian kernel

$$q_{t|s}(\cdot|x_s) = N\left(\frac{\alpha(T) - \alpha(t)}{\alpha(T) - \alpha(s)} x_s, \left(1 - \frac{(\alpha(T) - \alpha(t))^2}{(\alpha(T) - \alpha(s))^2}\right) I_d\right). \quad (16)$$

Since $p_T^*(x) = \int_{\mathbb{R}^d} q_{T|0}(x|x_0)d\pi(x_0)$, it results that $P_T^* = N(0, I_d)$.

*Proof.* The proof of (16) directly follows from Lemma 2. □
Following Lemma 7, the noising scheme given above is 'exact' with base distribution given by $\pi^{\text{base}} = N(0, I_d)$, this setting being quite unusual in diffusion model community. Based on this, we refer to this noising diffusion process as the Gaussianized Pinned Brownian Motion (GPBM). Moreover, under mild assumptions on $\pi$, the time-reverse SDE (2) writes as

$$dY_t = \frac{\zeta(T-t)}{\alpha(T) - \alpha(T-t)} \{Y_t + 2\nabla \log p_{T-t}^*(Y_t)\} dt + \sqrt{\frac{2\zeta(T-t)}{\alpha(T) - \alpha(T-t)}} dB_t, Y_0 \sim N(0, I_d).$$

# D VARIATIONAL OBJECTIVES FOR DIFFUSION-BASED METHODS

In this section, we first provide details on the continuous-time framework of variational diffusion-based methods (including RDS) in Appendix D.1. Then, we describe in Appendix D.2 the RDS discrete-time formulation, related to the variational loss $\mathcal{L}_{\text{RDS}}$ given in (7), and compare this framework to discrete-time formulations of previous variational methods.

In the rest of the section, we consider a target distribution $\pi \in \mathcal{P}(\mathbb{R}^d)$ and recall that $\mathbb{P}^*$ denotes the path measure associated to the noising SDE (1) initialized at $X_0 \sim \pi$.

## D.1 CONTINUOUS TIME FORMULATION

### D.1.1 RDS SETTING

Let $\pi^{\text{ref}} \in \mathcal{P}(\mathbb{R}^d)$. We recall that $\mathbb{P}^{\text{ref}}$ denotes the path measure associated to the noising SDE (1) starting at $X_0 \sim \pi^{\text{ref}}$. First note that $\mathbb{P}^* = \pi \otimes \mathbb{P}_{[0,T]}^{\text{ref}}$. Therefore we have the following relation $\mathbb{P}^* = \frac{d\pi}{d\pi^{\text{ref}}} \cdot \mathbb{P}^{\text{ref}}$, i.e., $\mathbb{P}^*$ is absolutely continuous with respect to $\mathbb{P}^{\text{ref}}$, with Radon-Nikodym derivative given by $(d\mathbb{P}^*/d\mathbb{P}^{\text{ref}})(X_{[0,T]}) = (d\pi/d\pi^{\text{ref}})(X_0)$, for any diffusion process $(X_t)_{t \in [0,T]}$ with forward time direction.

Below, we first provide the proof of Proposition 1.

*Proof of Proposition 1.* Assume that $\mathbb{P}_T^* = \mathbb{P}_T^{\text{ref}} = \pi^{\text{base}}$ and there exists $\theta^* \in \Theta$ such that $g_t^{\theta^*} = g_t$. Under those two assumptions, it is clear that $(Y_t^{\theta^*})_{t \in [0,T]} \sim (\mathbb{P}^*)^R$; therefore, we have $\mathcal{L}_{\text{LV}}(\theta^*) = 0$, i.e., $\theta^*$ achieves optimal solution. We now detail how to obtain (6).

Consider the detached diffusion process $(Y_t^{\hat{\theta}})_{t \in [0,T]} \sim \mathbb{P}^{\hat{\theta}}$. Using the fact that $\mathbb{P}^* = \frac{d\pi}{d\pi^{\text{ref}}} \cdot \mathbb{P}^{\text{ref}}$, the log-ratio in the objective (5) can be computed as

$$\log \frac{d\mathbb{P}^{\theta}}{d(\mathbb{P}^{*})^{R}}(Y_{[0,T]}^{\hat{\theta}}) = \log \frac{d\mathbb{P}^{\theta}}{d(\mathbb{P}^{\text{ref}})^{R}}(Y_{[0,T]}^{\hat{\theta}}) + \log \frac{\gamma^{\text{ref}}}{\gamma}(Y_{T}^{\hat{\theta}}) - \log \mathcal{Z}^{\text{ref}} + \log \mathcal{Z}. \quad (17)$$

Moreover, by applying Girsanov's theorem to the path measures $(\mathbb{P}^{\text{ref}})^R$ and $\mathbb{P}^{\theta}$ with the assumption $\mathbb{P}^{\text{ref}} = \pi^{\text{base}}$, we have

$$\log \frac{d\mathbb{P}^{\theta}}{d(\mathbb{P}^{\text{ref}})^R}(Y_{[0,T]}^{\hat{\theta}}) = \int_{0}^{T} \frac{\beta(T-t)}{2} \|g_{T-t}^{\theta}(Y_t^{\hat{\theta}})\|^2 dt + \int_{0}^{T} \sqrt{\beta(T-t)} g_{T-t}^{\theta}(Y_t^{\hat{\theta}})^{\top} dB_t. \quad (18)$$

Since the Radon-Nikodym derivative between path measures and their respective time-reversals is the same, we obtain the explicit LV objective combine the previous results to obtain the LV objective (6) by combining (17) and (18). $\square$

In contrast, previous reference-based works from Zhang & Chen (2022); Vargas et al. (2023a) relied on the reverse KL continuous-time objective

$$\mathcal{L}_{\text{KL}}^{c}(\theta) = \text{KL}(\mathbb{P}^{\theta} | (\mathbb{P}^{*})^{R}) = \mathbb{E} \left[ \log \frac{d\mathbb{P}^{\theta}}{d(\mathbb{P}^{*})^{R}} (Y_{[0,T]}^{\theta}) \right], \quad Y_{[0,T]}^{\theta} \sim \mathbb{P}^{\theta} \quad (19)$$

where the expectation is taken w.r.t. to the parameterized process, not the detached one. In this case too, the objective may be simplified as presented next.
**Proposition 9.** Assume there exists $\theta^* \in \Theta$ such that $g_t^{\theta^*} = g_t$. Then, the loss defined in (19) achieves optimal solution at $\theta^*$ and, setting $\rho = \log(\gamma^{\text{ref}}/\gamma)$, it simplifies as

$$ \mathcal{L}_{\text{KL}}^{c}(\theta) = \mathbb{E} \left[ \int_{0}^{T} \frac{\beta(T-t)}{2} \| g_{T-t}^{\theta}(Y_{t}^{\theta}) \|^2 dt + \rho(Y_{T}^{\theta}) \right], \quad (20) $$

where the equality holds, up to constants independent of $\theta$.

*Proof.* For any $\theta \in \Theta$, we have by the KL chain rule

$$ \mathcal{L}_{\text{KL}}^{c}(\theta) = \text{KL}(\pi^{\text{base}} | \mathbb{P}_{T}^{*}) + \mathbb{E}_{\pi^{\text{base}}} \left[ \text{KL}(\mathbb{P}_{0}^{\theta} | (\mathbb{P}^{*})_{0}^{R}) \right], $$

where the first term does not depend on $\theta$. Then, it is clear that $\mathcal{L}_{\text{KL}}^{c}$ is minimized when $\theta = \theta^*$. The result of (20) directly follows from (17) and (18); in particular, the additional constants are $-\log \mathcal{Z}^{\text{ref}}$ and $\log \mathcal{Z}$. $\square$

Based on the simplification (20), one may recognize a standard formulation of stochastic control-affine problem with terminal cost $\rho$, where the uncontrolled process is given by $\mathbb{P}^{\text{ref}}$.

### D.1.2 DIS SETTING

In their paper, Berner et al. (2023) seek to approximate the time-reversal of $\mathbb{P}^*$ for the VP noising scheme, see Appendix C.3. To do so, they consider the class of variational path measures $(\mathbb{P}^{\theta})_{\theta \in \Theta} \subset \mathcal{P}(C_T)$ where for any $\theta \in \Theta$, $\mathbb{P}^{\theta}$ is associated to the SDE

$$ dY_t^\theta = \frac{\beta(T-t)}{2} Y_t^\theta dt + \sigma^2 \beta(T-t) s_{T-t}^\theta (Y_t^\theta) dt + \sigma \sqrt{\beta(T-t)} dB_t, \quad Y_0^\theta \sim \text{N}(0, \sigma^2 I_d), \quad (21) $$

where the $\beta$-schedule is detailed in Appendix C.3, and $s^\theta: [0, T] \times \mathbb{R}^d$ is a neural network that aims at approximating the marginal scores of $\mathbb{P}^*$, i.e., $s^\theta(t, \cdot) \approx \nabla \log p_t^*$. The resulting variational method, DIS, then consists in minimizing the KL-based objective $\theta \mapsto \text{KL}(\mathbb{P}^\theta | (\mathbb{P}^*)^R)$, later extended by Richter et al. (2023) to the LV-based setting. We refer to (Richter et al., 2023, Section 3.2.) for detailed expressions of these continuous-time losses.

### D.1.3 CMCD SETTING

Alternatively, Vargas et al. (2024) propose to approximate a diffusion process with marginals prescribed by a curve of distributions $(\pi_t)_{t \in [0, T]}$, such that $\pi_0 = \pi^{\text{base}}$ and $\pi_T = \pi$, with unnnormalized densities that are assumed to be tractable. To do so, they consider the class of pairs of variational path measures $((\mathbb{P}^\theta, \mathbb{Q}^\theta))_{\theta \in \Theta} \subset \mathcal{P}(C_T)^2$ where, for any $\theta \in \Theta$, $\mathbb{Q}^\theta$ is associated to the SDE

$$ dX_t^\theta = \frac{\sigma^2}{2} (\nabla \log \pi_{T-t}(X_t^\theta) - h_{T-t}^\theta(X_t^\theta)) dt + \sigma dW_t, \quad X_0^\theta \sim \pi. \quad (22) $$

and $\mathbb{P}^\theta$ is associated to the SDE

$$ dY_t^\theta = \frac{\sigma^2}{2} (\nabla \log \pi_t(Y_t^\theta) + h_t^\theta(Y_t^\theta)) dt + \sigma dB_t, \quad Y_0^\theta \sim \pi^{\text{base}}, \quad (23) $$

with $h^\theta: [0, T] \times \mathbb{R}^d \to \mathbb{R}^d$ being a neural network. The resulting variational method, CMCD, consists in minimizing the KL-based objective $\theta \mapsto \text{KL}(\mathbb{P}^\theta | (\mathbb{Q}^\theta)^R)$ or the LV-based extension. We refer to the original paper for detailed expressions of these continuous-time losses.

## D.2 DISCRETE TIME FORMULATION

Following the setting of Section 2.2, we consider a time discretization of the interval $[0, T]$ given by an increasing sequence of timesteps $\{t_k\}_{k=0}^K$ such that $t_0 = 0^3$, $t_K = T$ and $K \ge 1$. In the rest of the section, we denote $\delta_k = t_{k+1} - t_k$.

<footnote n="3">³In the only case of PBM, we consider $t_0 > 0$ to avoid numerical integration error (typically $t_0 = 10^{-4}$).</footnote>
This section is organized as follows. We first explain in Appendix D.3 the discretization setting that we adopt for RDS, before deriving the RDS variational objective obtained by EM integration (Appendix D.3.1), EI integration in the PBM case (Appendix D.3.2) and EI integration in the VP case (Appendix D.3.3). Finally, for comparison purpose, we provide a discrete time approach to DIS in Appendix D.4 and CMCD in Appendix D.5. All of the presented settings are implemented in our code.

## D.3 GENERAL RDS SETTING

Let $\theta \in \Theta$ and let $\pi^{\text{ref}} \in \mathcal{P}(\mathbb{R}^d)$ be an arbitrary reference distribution. Consider the variational diffusion process $(Y_t^\theta)_{t \in [0, T]}$ defined by SDE (4). We recall that under mild assumptions on $\pi^{\text{ref}}$, see e.g., Cattiaux et al. (2023), the time-reversal of $\mathbb{P}^{\text{ref}}$ is associated to the SDE

$$dY_t^{\text{ref}} = -f(T-t)Y_t^{\text{ref}}dt + \beta(T-t)s_{T-t}^{\text{ref}}(Y_t^{\text{ref}})dt + \sqrt{\beta(T-t)}dB_t, \quad Y_0^{\text{ref}} \sim \mathbb{P}_T^{\text{ref}}. \quad (24)$$

Assume that we are given transition kernels $(A, y_k) \mapsto p_{k+1|k}^{\theta}(A|y_k)$ and $(A, y_k) \mapsto p_{k+1|k}^{\text{ref}}(A|y_k)$ respectively approximating $\mathbb{P}(Y_{t_{k+1}}^{\theta} \in A | Y_{t_k}^{\theta} = y_k)$ and $\mathbb{P}(Y_{t_{k+1}}^{\text{ref}} \in A | Y_{t_k}^{\text{ref}} = y_k)$. Then, for $K$ sufficiently large, the path measures $\mathbb{P}^{\theta}$ and $(\mathbb{P}^{\text{ref}})^R$ may respectively be approximated by the joint distributions $p_{0:K}^{\theta} \in \mathcal{P}^{(K+1)}$ and $p_{0:K}^{\text{ref}} \in \mathcal{P}^{(K+1)}$, defined by $p_{0:K}^{\theta} = \pi^{\text{base}} \prod_{k=0}^{K-1} p_{k+1|k}^{\theta}$ and $p_{0:K}^{\text{ref}} = \pi^{\text{base}} \prod_{k=0}^{K-1} p_{k+1|k}^{\text{ref}}$. Using the relation $\mathbb{P}^* = \frac{d\pi}{d\pi^{\text{ref}}} \cdot \mathbb{P}^{\text{ref}}$, $(\mathbb{P}^*)^R$ may also be approximated by the joint distribution $p_{0:K}^{\star} = \frac{d\pi}{d\pi^{\text{ref}}} \cdot p_{0:K}^{\text{ref}} \in \mathcal{P}^{(K+1)}$.

Using this correspondence between path measures and joint distributions, we propose to approximate $\mathcal{L}_{\text{LV}}$ defined in (5) by the discrete time LV-based objective

$$\mathcal{L}_{\text{RDS}}(\theta) = \text{Var} \left[ \log \frac{dp_{0:K}^{\theta}(Y_{0:K}^{\theta})}{dp_{0:K}^{\star}} \right], \quad Y_{0:K}^{\theta} \sim p_{0:K}^{\theta}. \quad (25)$$

In a similar manner, $\mathcal{L}_{\text{KL}}^c$ defined in (19) may be approximated by the discrete time KL-based objective

$$\mathcal{L}_{\text{KL}}(\theta) = \mathbb{E} \left[ \log \frac{dp_{0:K}^{\theta}}{dp_{0:K}^{\star}}(Y_{0:K}^{\theta}) \right], \quad Y_{0:K}^{\theta} \sim p_{0:K}^{\theta}. \quad (26)$$

By definition of the joint distributions, and similarly to the continuous-time relation, the log-ratio in (25) and (26) can be simplified as

$$\log \frac{dp_{0:K}^{\theta}}{dp_{0:K}^{\star}}(Y_{0:K}^{\bar{\theta}}) = \log \frac{dp_{0:K}^{\theta}}{dp_{0:K}^{\text{ref}}}(Y_{0:K}^{\bar{\theta}}) + \log \frac{\gamma^{\text{ref}}}{\gamma}(Y_{K}^{\bar{\theta}}) + \log \mathcal{Z} - \log \mathcal{Z}^{\text{ref}}. \quad (27)$$

where $\bar{\theta} \in \{\theta, \hat{\theta}\}$. In the following sections, building upon (27), we derive the exact expression of the variational loss (25) for each combination of noising and discretization schemes. We emphasize that these derivations are in particular applicable to PIS and DDS settings by taking $\pi^{\text{ref}}$ as given in Table 1.

### D.3.1 EM-BASED RDS

In this section, we consider an arbitrary noising scheme that fits within the framework detailed in Appendix C.1. For any $k \in \{0, \dots, K-1\}$, define the coefficients

$$w_k^{\text{EM}} = \beta(T - t_k)\delta_k,$$

$$a_k^{\text{EM}} = 1 - f(T - t_k)\delta_k,$$

$$b_k^{\text{EM}} = \beta(T - t_k)\delta_k,$$

$$c_k^{\text{EM}} = \beta(T - t_k)\delta_k.$$

To approximate SDEs (24) and (4) on time interval $[t_k, t_{k+1}]$, we may define the following transition kernels

$$p_{k+1|k}^{\text{ref}}(\cdot|y_k) = N(a_k^{\text{EM}} y_k + b_k^{\text{EM}} s_{T-t_k}^{\text{ref}}(y_k), c_k^{\text{EM}} I_d), \quad (28)$$

$$p_{k+1|k}^{\theta}(\cdot|y_k) = N(a_k^{\text{EM}} y_k + b_k^{\text{EM}} \{s_{T-t_k}^{\text{ref}}(y_k) + g_{T-t_k}^{\theta}(y_k)\}, c_k^{\text{EM}} I_d), \quad (29)$$
that can be obtained by applying the EM integration scheme, see Appendix B. In this case, the log-ratio $\log dp^\theta/dp^{\text{ref}}$ can be computed as follows.

**Lemma 10.** Assume that the joint distributions $p^\theta$ and $p^{\text{ref}}$ are respectively induced by the transition kernels defined in (28) and (29). Then, for any $y_{0:K} \in (\mathbb{R}^d)^{K+1}$, we have

$$ \log \frac{dp^\theta}{dp^{\text{ref}}}(y_{0:K}) = \sum_{k=0}^{K-1} \frac{w_k^{\text{EM}}}{2} \|g_{T-t_k}^\theta(y_k)\|^2 + \sum_{k=0}^{K-1} \sqrt{w_k^{\text{EM}}} g_{T-t_k}^\theta(y_k)^\top \epsilon_k, \\ \text{where } \epsilon_k = \frac{1}{\sqrt{c_k^{\text{EM}}}} (y_{k+1} - a_k^{\text{EM}} y_k - b_k^{\text{EM}} \{s_{T-t_k}^{\text{ref}}(y_k) + g_{T-t_k}^\theta(y_k)\}). \quad (30) $$

*Proof.* In this proof, for ease of reading, we will respectively denote $a_k^{\text{EM}}$, $b_k^{\text{EM}}$, $c_k^{\text{EM}}$ and $w_k^{\text{EM}}$ by $a_k$, $b_k$, $c_k$ and $w_k$. Let $y_{0:K} \in (\mathbb{R}^d)^{K+1}$. We have

$$ \log \frac{dp^\theta}{dp^{\text{ref}}}(y_{0:K}) = \sum_{k=0}^{K-1} \log \frac{dp_{k+1|k}^\theta(y_{k+1}|y_k)}{dp_{k+1|k}^{\text{ref}}(y_{k+1}|y_k)} \\ = \sum_{k=0}^{K-1} \left\{ -\frac{1}{2} \|\epsilon_k\|^2 + \frac{1}{2} \|\epsilon_k + (b_k/\sqrt{c_k})g_{T-t_k}^\theta(y_k)\|^2 \right\} \\ = \sum_{k=0}^{K-1} \frac{b_k^2}{2c_k} \|g_{T-t_k}^\theta(y_k)\|^2 + \sum_{k=0}^{K-1} \frac{b_k}{\sqrt{c_k}} g_{T-t_k}^\theta(y_k)^\top \epsilon_k, $$

where $\epsilon_k$ is defined by (30). The result of Lemma 10 directly follows by noting that $w_k = b_k^2/c_k$. $\square$

By combining Lemma 10 with (27), we deduce the expression of the RDS loss in this setting.

**Corollary 11 (Expression of RDS loss – EM setting).** When using the EM scheme to integrate the SDE (4), the RDS loss simplifies as

$$ \mathcal{L}_{\text{RDS}}(\theta) = \text{Var} \left[ \sum_{k=0}^{K-1} w_k^{\text{EM}} g_{T-t_k}^{\theta}(Y_k)^{\top} \{g_{T-t_k}^{\hat{\theta}}(Y_k) - \frac{1}{2} g_{T-t_k}^{\theta}(Y_k)\} + \sum_{k=0}^{K-1} \sqrt{w_k^{\text{EM}}} g_{T-t_k}^{\theta}(Y_k)^{\top} Z_k + \log \frac{\gamma^{\text{ref}}}{\gamma}(Y_K) \right], $$

where $\{Z_k\}_{k=0}^{K-1}$ are independently distributed according to $N(0, I_d)$, and $\{Y_k\}_{k=0}^K$ is recursively defined by $Y_0 \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$,

$$ Y_{k+1} = a_k^{\text{EM}} Y_k + b_k^{\text{EM}} \{s_{T-t_k}^{\text{ref}}(Y_k) + g_{T-t_k}^{\hat{\theta}}(Y_k)\} + \sqrt{c_k^{\text{EM}}} Z_k, \quad (31) $$

with $\hat{\theta}$ being a detached version of $\theta$.

*Proof.* Consider the discrete time process $\{Y_k\}_{k=0}^K$ defined by recursion (31). With the same notation as in Lemma 10, we have for any $k \in \{0, \dots, K-1\}$

$$ \epsilon_k = Z_k + \sqrt{w_k^{\text{EM}}} \{g_{T-t_k}^{\hat{\theta}}(Y_k) - g_{T-t_k}^{\theta}(Y_k)\}, $$

recalling that $w_k^{\text{EM}} = (b_k^{\text{EM}})^2/c_k^{\text{EM}}$. Using (27), we finally obtain the expression of the RDS loss. $\square$

For completeness, we also provide a similar result for the KL-based objective given in (26).

**Corollary 12 (Expression of reverse KL loss – EM setting).** When using the EM scheme to integrate the SDE (4), up to additional constants independent of $\theta$, the discrete time KL loss simplifies as

$$ \mathcal{L}_{\text{KL}}(\theta) = \mathbb{E} \left[ \sum_{k=0}^{K-1} \frac{w_k^{\text{EM}}}{2} \|g_{T-t_k}^{\theta}(Y_k^{\theta})\|^2 + \log \frac{\gamma^{\text{ref}}}{\gamma}(Y_K^{\theta}) \right], $$
where {$Z_k$}$_{k=0}^{K-1}$ are independently distributed according to N(0, I$_d$), and {$Y_k^\theta$}$_{k=0}^K$ is recursively defined by $Y_0^\theta \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$

$$Y_{k+1}^{\theta} = a_{k}^{\text{EM}} Y_{k}^{\theta} + b_{k}^{\text{EM}} \{s_{T-t_k}^{\text{ref}}(Y_k^{\theta}) + g_{T-t_k}^{\theta}(Y_k^{\theta})\} + \sqrt{c_{k}^{\text{EM}}} Z_k.$$

As observed by Vargas et al. (2023a), there may be a significant bias raising from the use of the Euler-Maruyama discretization when using diffusion-based sampling methods. This motivates us to consider EI discretization in the RDS setting.

### D.3.2 EI-BASED RDS (PBM)

We first provide the expression of the RDS loss in the case of the PBM noising scheme, which is detailed in Appendix C.2. We recall the notation $\alpha(t) = \int_0^t \beta(u)du$. For any $k \in \{0, \dots, K-1\}$, define the coefficients

$$\begin{aligned}
w_k^{\text{PBM}} &= \frac{(\alpha(T) - \alpha(T - t_k))(\alpha(T - t_k) - \alpha(T - t_{k+1}))}{\alpha(T) - \alpha(T - t_{k+1})}, \\
a_k^{\text{PBM}} &= \frac{\alpha(T) - \alpha(T - t_{k+1})}{\alpha(T) - \alpha(T - t_k)}, \\
b_k^{\text{PBM}} &= \alpha(T - t_k) - \alpha(T - t_{k+1}), \\
c_k^{\text{PBM}} &= \frac{(\alpha(T) - \alpha(T - t_{k+1}))(\alpha(T - t_k) - \alpha(T - t_{k+1}))}{\alpha(T) - \alpha(T - t_k)}.
\end{aligned}$$

To approximate SDEs (24) and (4) on time interval $[t_k, t_{k+1}]$, we may define the following transition kernels

$$p_{k+1|k}^{\text{ref}}(\cdot|y_k) = N(a_k^{\text{PBM}} y_k + b_k^{\text{PBM}} s_{T-t_k}^{\text{ref}}(y_k), c_k^{\text{PBM}} I_d), \quad (32)$$

$$p_{k+1|k}^{\theta}(\cdot|y_k) = N(a_k^{\text{PBM}} y_k + b_k^{\text{PBM}} \{s_{T-t_k}^{\text{ref}}(y_k) + g_{T-t_k}^{\theta}(y_k)\}, c_k^{\text{PBM}} I_d), \quad (33)$$

that can be obtained by applying the EI scheme, see Appendix B. In this case, the log ratio $\log dp^\theta/dp^{\text{ref}}$ can be computed as follows.

**Lemma 13.** Assume that the joint distributions $p^\theta$ and $p^{\text{ref}}$ are respectively induced by the transition kernels defined in (32) and (33). Then for any $y_{0:K} \in (\mathbb{R}^d)^{K+1}$, we have

$$\log \frac{dp^\theta}{dp^{\text{ref}}}(y_{0:K}) = \sum_{k=0}^{K-1} \frac{w_k^{\text{PBM}}}{2} \|g_{T-t_k}^\theta(y_k)\|^2 + \sum_{k=0}^{K-1} \sqrt{w_k^{\text{PBM}}} g_{T-t_k}^\theta(y_k)^\top \epsilon_k, \\ \text{where } \epsilon_k = \frac{1}{\sqrt{c_k^{\text{PBM}}} \{y_{k+1} - a_k^{\text{PBM}} y_k - b_k^{\text{PBM}} \{s_{T-t_k}^{\text{ref}}(y_k) + g_{T-t_k}^\theta(y_k)\}\}}.$$

*Proof.* The proof is exactly the same as Lemma 10, noting that we still have $w_k^{\text{PBM}} = (b_k^{\text{PBM}})^2/c_k^{\text{PBM}}$. $\square$

By combining Lemma 13 with (27), we deduce the expression of the RDS loss in this setting. The proof of the following is the same as in Corollary 11.

**Corollary 14 (Expression of RDS loss - EI setting - PBM scheme).** When using the EI scheme to integrate SDE (4) implemented with PBM, the RDS loss simplifies as

$$\mathcal{L}_{\text{RDS}}(\theta) = \text{Var} \left[ \sum_{k=0}^{K-1} w_k^{\text{PBM}} g_{T-t_k}^{\theta}(Y_k)^{\top} \{g_{T-t_k}^{\hat{\theta}}(Y_k) - \frac{1}{2} g_{T-t_k}^{\theta}(Y_k)\} \right. \\ \left. + \sum_{k=0}^{K-1} \sqrt{w_k^{\text{PBM}}} g_{T-t_k}^{\theta}(Y_k)^{\top} Z_k + \log \frac{\gamma^{\text{ref}}}{\gamma}(Y_K) \right],$$
where {$Z_k$}$_{k=0}^{K-1}$ are independently distributed according to N(0, I$_d$), and {$Y_k$}$_{k=0}^{K}$ is recursively defined by $Y_0 \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$

$$Y_{k+1} = a_k^{\text{PBM}} Y_k + b_k^{\text{PBM}} \{s_{T-t_k}^{\text{ref}}(Y_k) + g_{T-t_k}^{\hat{\theta}}(Y_k)\} + \sqrt{c_k^{\text{PBM}}} Z_k,$$

with $\hat{\theta}$ being a detached version of $\theta$.

For completeness, we also provide a similar result for the KL-based objective given in (26).

**Corollary 15** (Expression of reverse KL loss – EI setting – PBM noising). When using the EI scheme to integrate SDE (4) implemented with PBM, up to additional constants independent of $\theta$, the discrete time KL loss simplifies as

$$\mathcal{L}_{KL}(\theta) = \mathbb{E} \left[ \sum_{k=0}^{K-1} \frac{w_k^{\text{PBM}}}{2} \|g_{T-t_k}^{\theta}(Y_k^{\theta})\|^2 + \log \frac{\gamma^{\text{ref}}}{\gamma} (Y_K^{\theta}) \right],$$

where {$Z_k$}$_{k=0}^{K-1}$ are independently distributed according to N(0, I$_d$), and {$Y_k^{\theta}$}$_{k=0}^{K}$ is recursively defined by $Y_0^{\theta} \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$

$$Y_{k+1}^{\theta} = a_k^{\text{PBM}} Y_k^{\theta} + b_k^{\text{PBM}} \{s_{T-t_k}^{\text{ref}}(Y_k^{\theta}) + g_{T-t_k}^{\theta}(Y_k^{\theta})\} + \sqrt{c_k^{\text{PBM}}} Z_k.$$

### D.3.3 EI-BASED RDS (VP)

Here, we provide the expression of the RDS loss in the case of the VP noising scheme, which is detailed in Appendix C.3. We recall the notation $\alpha(t) = \int_0^t \beta(u)du$. For any $k \in \{0, \dots, K-1\}$, define $\lambda_k = \exp((\alpha(T - t_k) - \alpha(T - t_{k+1})) - 1$ and the coefficients

$$w_k^{\text{VP}} = \frac{4\sigma^2(\sqrt{1+\lambda_k}-1)^2}{\lambda_k} = 4\sigma^2 \tanh \left( \frac{\alpha(T-t_k)-\alpha(T-t_{k+1})}{4} \right)$$

$$a_k^{\text{VP}} = \sqrt{1 + \lambda_k}$$

$$b_k^{\text{VP}} = 2\sigma^2\{\sqrt{1 + \lambda_k} - 1\}$$

$$c_k^{\text{VP}} = \sigma^2 \lambda_k$$

To approximate SDEs (24) and (4) on time interval $[t_k, t_{k+1}]$, we may define the following transition kernels

$$p_{k+1|k}^{\text{ref}}(\cdot|y_k) = N(a_k^{\text{VP}} y_k + b_k^{\text{VP}} s_{T-t_k}^{\text{ref}}(y_k), c_k^{\text{VP}} I_d), \quad (34)$$

$$p_{k+1|k}^{\theta}(\cdot|y_k) = N(a_k^{\text{VP}} y_k + b_k^{\text{VP}} \{s_{T-t_k}^{\text{ref}}(y_k) + g_{T-t_k}^{\theta}(y_k)\}, c_k^{\text{VP}} I_d),$$

that can be obtained by applying the EI scheme, see Appendix B. In this case, the log ratio $\log dp^{\theta}/dp^{\text{ref}}$ can be computed as follows.

**Lemma 16.** Assume that the joint distributions $p^{\theta}$ and $p^{\text{ref}}$ are respectively induced by the transition kernels defined in (34) and (34). Then, for any $y_{0:K} \in (\mathbb{R}^d)^{K+1}$, we have

$$\log \frac{dp^{\theta}}{dp^{\text{ref}}}(y_{0:K}) = \sum_{k=0}^{K-1} \frac{w_k^{\text{VP}}}{2} \|g_{T-t_k}^{\theta}(y_k)\|^2 + \sum_{k=0}^{K-1} \sqrt{w_k^{\text{VP}}} g_{T-t_k}^{\theta}(y_k)^{\top} \epsilon_k$$
$$\text{where } \epsilon_k = \frac{1}{\sqrt{c_k^{\text{VP}}}} (y_{k+1} - a_k^{\text{VP}} y_k - b_k^{\text{VP}} \{s_{T-t_k}^{\text{ref}}(y_k) + g_{T-t_k}^{\theta}(y_k)\})$$

*Proof.* The proof is exactly the same as Lemma 10, noting that we still have $w_k^{\text{VP}} = (b_k^{\text{VP}})^2/c_k^{\text{VP}}$. $\square$

By combining Lemma 16 with (27), we deduce the expression of the RDS loss in this setting. The proof of the following is the same as in Corollary 11.
**Corollary 17 (Expression of RDS loss - EI setting - VP scheme).** When using the EI scheme to integrate SDE (4) implemented with VP, the RDS loss simplifies as

$$ \mathcal{L}_{RDS}(\theta) = \mathrm{Var} \left[ \sum_{k=0}^{K-1} w_k^{\mathrm{VP}} g_{T-t_k}^{\theta}(Y_k)^{\top} \{g_{T-t_k}^{\hat{\theta}}(Y_k) - \frac{1}{2} g_{T-t_k}^{\theta}(Y_k)\} + \sum_{k=0}^{K-1} \sqrt{w_k^{\mathrm{VP}}} g_{T-t_k}^{\theta}(Y_k)^{\top} Z_k + \log \frac{\gamma^{\mathrm{ref}}}{\gamma}(Y_K) \right], $$

where $\{Z_k\}_{k=0}^{K-1}$ are independently distributed according to $N(0, I_d)$, and $\{Y_k\}_{k=0}^{K}$ is recursively defined by $Y_0 \sim \pi^{\mathrm{base}}$, and for any $k \in \{0, \dots, K-1\}$

$$ Y_{k+1} = a_k^{\mathrm{VP}} Y_k + b_k^{\mathrm{VP}} \{s_{T-t_k}^{\mathrm{ref}}(Y_k) + g_{T-t_k}^{\hat{\theta}}(Y_k)\} + \sqrt{c_k^{\mathrm{VP}}} Z_k, $$

with $\hat{\theta}$ being a detached version of $\theta$.

For completeness, we also provide a similar result for the KL-based objective given in (26).

**Corollary 18 (Expression of reverse KL loss - EI setting - VP noising).** When using the EI scheme to integrate SDE (4) implemented with VP, up to additional constants independent of $\theta$, the discrete time KL loss simplifies as

$$ \mathcal{L}_{KL}(\theta) = \mathbb{E} \left[ \sum_{k=0}^{K-1} \frac{w_k^{\mathrm{VP}}}{2} \|g_{T-t_k}^{\theta}(Y_k^{\theta})\|^2 + \log \frac{\gamma^{\mathrm{ref}}}{\gamma}(Y_K^{\theta}) \right], $$

where $\{Z_k\}_{k=0}^{K-1}$ are independently distributed according to $N(0, I_d)$, and $\{Y_k^{\theta}\}_{k=0}^{K}$ is recursively defined by $Y_0^{\theta} \sim \pi^{\mathrm{base}}$, and for any $k \in \{0, \dots, K-1\}$,

$$ Y_{k+1}^{\theta} = a_k^{\mathrm{VP}} Y_k^{\theta} + b_k^{\mathrm{VP}} \{s_{T-t_k}^{\mathrm{ref}}(Y_k^{\theta}) + g_{T-t_k}^{\theta}(Y_k^{\theta})\} + \sqrt{c_k^{\mathrm{VP}}} Z_k. $$

## D.4 DIS SETTING

Here, we propose a discrete time alternative to the continuous time objective proposed by Berner et al. (2023) and Richter et al. (2023). Let $\theta \in \Theta$. Consider the variational diffusion process $(Y_t^{\theta})_{t \in [0, T]}$ defined by SDE (21).

Assume that we are given transition kernels $(A, y_k) \mapsto p_{k+1|k}^{\theta}(A|y_k)$ and $(A, y_{k+1}) \mapsto p_{k|k+1}^{\theta}(A|y_{k+1})$ respectively approximating $\mathbb{P}(Y_{t_{k+1}}^{\theta} \in A | Y_{t_k}^{\theta} = y_k)$ and $\mathbb{P}(X_{t_k} \in A | X_{t_{k+1}} = y_{k+1})$. Then, for $K$ sufficiently large, the path measures $\mathbb{P}^{\theta}$ and $(\mathbb{P}^{\star})^{R}$ may be approximated by the joint distributions $p_{0:K}^{\theta} \in \mathcal{P}^{(K+1)}$ and $p_{0:K}^{\star} \in \mathcal{P}^{(K+1)}$ defined by $p_{0:K}^{\theta} = \pi^{\mathrm{base}} \prod_{k=0}^{K-1} p_{k+1|k}^{\theta}$ and $p_{0:K}^{\star} = \pi \prod_{k=0}^{K-1} p_{k|k+1}^{\star}$.

Following the RDS methodology, we propose to approximate the KL-based continuous-time DIS objective by

$$ \mathcal{L}_{\mathrm{KL}}^{\mathrm{DIS}}(\theta) = \mathbb{E} \left[ \log \frac{\mathrm{d}p_{0:K}^{\theta}}{\mathrm{d}p_{0:K}^{\star}}(Y_{0:K}^{\theta}) \right], \quad Y_{0:K}^{\theta} \sim p_{0:K}^{\theta}, \qquad (35) $$

and the LV-based continuous-time DIS objective by

$$ \mathcal{L}_{\mathrm{LV}}^{\mathrm{DIS}}(\theta) = \mathrm{Var} \left[ \log \frac{\mathrm{d}p_{0:K}^{\theta}}{\mathrm{d}p_{0:K}^{\star}}(Y_{0:K}^{\theta}) \right], \quad Y_{0:K}^{\theta} \sim p_{0:K}^{\theta}. \qquad (36) $$

To approximate SDEs (21) and (1) on time interval $[t_k, t_{k+1}]$ in the case of the VP noising scheme, we may define the following transition kernels

$$ p_{k|k+1}^{\star}(\cdot|y_{k+1}) = \mathrm{N}(\sqrt{1 - \lambda_{T-t_{k+1},T-t_k}} y_{k+1}, \sigma^2 \lambda_{T-t_{k+1},T-t_k} I_d), \qquad (37) $$

$$ p_{k+1|k}^{\theta}(\cdot|y_k) = \mathrm{N}(a_k^{\mathrm{VP}} y_k + b_k^{\mathrm{VP}} + s_{T-t_k}^{\theta}(y_k), c_k^{\mathrm{VP}} I_d), \qquad (38) $$
where $\lambda_{s,t}$ is defined in Lemma 6 with $s = T - t_{k+1}$ and $t = T - t_k$, and $\{a_k^{\text{VP}}, b_k^{\text{VP}}, c_k^{\text{VP}}\}$ are defined via the VP-EI integration scheme detailed in Appendix D.3.3. In particular, the following equalities hold:

$$ 
\begin{aligned}
\sqrt{1 - \lambda_{T-t_{k+1}, T-t_k}} &= \exp\left(\frac{1}{2}[\alpha(T - t_{k+1}) - \alpha(T - t_k)]\right) = \frac{1}{a_k^{\text{VP}}}, \\
\sigma^2 \lambda_{T-t_{k+1}, T-t_k} &= \sigma^2[1 - \exp(\alpha(T - t_{k+1}) - \alpha(T - t_k))] = \frac{c_k^{\text{VP}}}{(a_k^{\text{VP}})^2}.
\end{aligned}
$$

Therefore, we have $p_{k|k+1}^*(y_k|y_{k+1}) = N(y_{k+1}; a_k^{\text{VP}}y_k, c_k^{\text{VP}}I_d)$. In this case, the log-ratio $\log dp^\theta/dp^*$ can be computed as follows.

**Lemma 19.** Assume that the joint distributions $p^*$ and $p^\theta$ are respectively induced by the transition kernels defined in (37) and (38). Then, for any $y_{0:K} \in (\mathbb{R}^d)^{K+1}$, we have

$$ 
\begin{aligned}
\log \frac{dp^\theta}{dp^*}(y_{0:K}) &= \log Z + \log \frac{\pi^{\text{base}}(y_0)}{\gamma(y_K)} + \sum_{k=0}^{K-1} \log \frac{dp_{k+1|k}^\theta(y_{k+1}|y_k)}{dp_{k|k+1}^*(y_k|y_{k+1})} \\
&= \log Z + \log \frac{\pi^{\text{base}}(y_0)}{\gamma(y_K)} + \sum_{k=0}^{K-1} \frac{w_k^{\text{VP}}}{2} \|s_{T-t_k}^\theta(y_k)\|^2 + \sum_{k=0}^{K-1} \sqrt{w_k^{\text{VP}}} s_{T-t_k}^\theta(y_k)^\top \epsilon_k \\
\text{where } \epsilon_k &= \frac{1}{\sqrt{c_k^{\text{VP}}}} (y_{k+1} - a_k^{\text{VP}}y_k - b_k^{\text{VP}} s_{T-t_k}^\theta(y_k))
\end{aligned}
$$

*Proof.* This result lies on the following identity

$$ \frac{\mathrm{d}p_{k+1|k}^{\theta}(y_{k+1}|y_k)}{\mathrm{d}p_{k|k+1}^{*}(y_k|y_{k+1})} = -\frac{1}{2} \|\epsilon_k\|^2 + \frac{1}{2} \left\| \epsilon_k + \sqrt{w_k^{\mathrm{VP}}} s_{T-t_k}^{\theta}(y_k) \right\|^2 \quad \square $$

Based on Lemma 19, we simply deduce the expression of the discrete time versions of LV-based and KL-based DIS losses.

**Corollary 20 (Expression of LV-DIS loss).** When using the EI scheme to integrate SDE (21), the LV-based discrete time DIS loss (36) simplifies as

$$ \mathcal{L}_{\text{LV}}^{\text{DIS}}(\theta) = \text{Var} \left[ \sum_{k=0}^{K-1} w_k^{\text{VP}} s_{T-t_k}^{\theta}(Y_k)^{\top} \{s_{T-t_k}^{\hat{\theta}}(Y_k) - \frac{1}{2} s_{T-t_k}^{\theta}(Y_k)\} \right. \\ \left. + \sum_{k=0}^{K-1} \sqrt{w_k^{\text{VP}}} s_{T-t_k}^{\theta}(Y_k)^{\top} Z_k + \log \frac{\pi^{\text{base}}(Y_0)}{\gamma(Y_K)} \right], $$

where $\{Z_k\}_{k=0}^{K-1}$ are independently distributed according to $N(0, I_d)$, and $\{Y_k\}_{k=0}^K$ is recursively defined by $Y_0 \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$,

$$ Y_{k+1} = a_k^{\text{VP}}Y_k + b_k^{\text{VP}} s_{T-t_k}^{\hat{\theta}}(Y_k) + \sqrt{c_k^{\text{VP}}} Z_k, $$

with $\hat{\theta}$ being a detached version of $\theta$.

*Proof.* This result is an application from Lemma 19, where we have $\epsilon_k = Z_k + \sqrt{w_k^{\text{VP}}} \{s_{T-t_k}^{\hat{\theta}}(Y_k) - s_{T-t_k}^{\theta}(Y_k)\}$ for any $k \in \{0, \dots, K-1\}$. $\square$

**Corollary 21 (Expression of KL-DIS loss).** When using the EM scheme to integrate SDE (21), up to additional constants independent of $\theta$, the KL-based discrete time DIS loss (35) simplifies as

$$ \mathcal{L}_{\text{KL}}^{\text{DIS}}(\theta) = \mathbb{E} \left[ \sum_{k=0}^{K-1} \frac{w_k^{\text{VP}}}{2} \left\| s_{T-t_k}^{\theta}(Y_k^{\theta}) \right\|^2 + \log \frac{\pi^{\text{base}}(Y_0^{\theta})}{\gamma(Y_K^{\theta})} \right], $$
where {$Z_k$}$_{k=0}^{K-1}$ are independently distributed according to N(0, I$_d$), and {$Y_k^\theta$}$_{k=0}^K$ is recursively defined by $Y_0^\theta \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$,

$$Y_{k+1}^{\theta} = a_{k}^{\text{VP}} Y_{k}^{\theta} + b_{k}^{\text{VP}} s_{T-t_{k}}^{\theta} (Y_{k}^{\theta}) + \sqrt{c_{k}^{\text{VP}}} Z_{k}.$$

## D.5 CMCD SETTING

Here, we propose a discrete time alternative to the continuous time objective proposed by Vargas et al. (2024). Let $\theta \in \Theta$. Consider the variational diffusion processes $(Y_t^\theta)_{t \in [0, T]}$ defined by SDE (23), and $(X_t^\theta)_{t \in [0, T]}$, defined by (22).

Assume that we are given transition kernels $(\mathcal{A}, y_k) \mapsto p_{k+1|k}^\theta(\mathcal{A}|y_k)$ and $(\mathcal{A}, y_{k+1}) \mapsto q_{k|k+1}^\theta(\mathcal{A}|y_{k+1})$ respectively approximating $\mathbb{P}(Y_{t_{k+1}}^\theta \in \mathcal{A} | Y_{t_k}^\theta = y_k)$ and $\mathbb{Q}^\theta(X_{t_k}^\theta \in \mathcal{A} | X_{t_{k+1}}^\theta = y_{k+1})$. Then, for $K$ sufficiently large, the path measures $\mathbb{P}^\theta$ and $(\mathbb{Q}^\theta)^R$ may be approximated by the joint distributions $p_{0:K}^\theta \in \mathcal{P}^{(K+1)}$ and $q_{0:K}^\theta \in \mathcal{P}^{(K+1)}$ defined by $p_{0:K}^\theta = \pi^{\text{base}} \prod_{k=0}^{K-1} p_{k+1|k}^\theta$ and $q_{0:K}^\theta = \pi \prod_{k=0}^{K-1} q_{k|k+1}^\theta$.

Following the RDS methodology, we propose to approximate the KL-based continuous-time CMCD objective by

$$\mathcal{L}_{\text{KL}}^{\text{CMCD}}(\theta) = \mathbb{E} \left[ \log \frac{dp_{0:K}^{\theta}(Y_{0:K}^{\theta})}{dq_{0:K}^{\theta}} \right], \quad Y_{0:K}^{\theta} \sim p_{0:K}^{\theta}, \qquad (39)$$

and the LV-based continuous-time CMCD objective by

$$\mathcal{L}_{\text{LV}}^{\text{CMCD}}(\theta) = \text{Var} \left[ \log \frac{dp_{0:K}^{\theta}(Y_{0:K}^{\theta})}{dq_{0:K}^{\theta}} \right], \quad Y_{0:K}^{\theta} \sim p_{0:K}^{\theta}. \qquad (40)$$

For any $k \in \{0, \dots, K-1\}$, let us define $w_k^{\text{CMCD}} = \sigma^2 \delta_k / 4$.

To approximate SDEs (22) and (23) on time interval $[t_k, t_{k+1}]$, we may define the following transition kernels

$$q_{k|k+1}^{\theta}(\cdot|y_{k+1}) = N(y_{k+1} + \frac{\sigma^2 \delta_k}{2} (\nabla \log \pi_{t_{k+1}}(y_{k+1}) - h_{t_{k+1}}^{\theta}(y_{k+1})), \sigma^2 \delta_k I_d), \quad (41)$$

$$p_{k+1|k}^{\theta}(\cdot|y_k) = N(y_k + \frac{\sigma^2 \delta_k}{2} (\nabla \log \pi_{t_k}(y_k) + h_{t_k}^{\theta}(y_k)), \sigma^2 \delta_k I_d), \quad (42)$$

that can be obtained by applying the EM integration scheme, see Appendix B. In this case, the log-ratio $\log dp^\theta/dq^\theta$ can be computed as follows.

**Lemma 22.** Assume that the joint distributions $q^\theta$ and $p^\theta$ are respectively induced by the transition kernels defined in (41) and (42). Then, for any $y_{0:K} \in (\mathbb{R}^d)^{K+1}$, we have

$$\log \frac{dp^\theta}{dq^\theta}(y_{0:K}) = \log Z + \log \frac{\pi^{\text{base}}(y_0)}{\gamma(y_K)} + \sum_{k=0}^{K-1} \log \frac{dp_{k+1|k}^\theta(y_{k+1}|y_k)}{dq_{k|k+1}^\theta(y_k|y_{k+1})} \\ = \log Z + \log \frac{\pi^{\text{base}}(y_0)}{\gamma(y_K)} + \frac{1}{2} \sum_{k=0}^{K-1} w_k^{\text{CMCD}} \|u_k\|^2 + \sum_{k=0}^{K-1} \sqrt{w_k^{\text{CMCD}}} u_k^\top \epsilon_k,$$

where

$$u_k = \nabla \log \pi_{t_k}(y_k) + \nabla \log \pi_{t_{k+1}}(y_{k+1}) + h_{t_k}^{\theta}(y_k) - h_{t_{k+1}}^{\theta}(y_{k+1}), \\ \epsilon_k = \frac{1}{\sigma\sqrt{\delta_k}} \left( y_{k+1} - y_k - \frac{\sigma^2\delta_k}{2} \left\{ \nabla \log \pi_{t_k}(y_k) + h_{t_k}^{\theta}(y_k) \right\} \right).$$

**Proof.** Based on the notation of $u_k$ and $\epsilon_k$, we obtain the result by noting that

$$\log \frac{dp_{k+1|k}^{\theta}(y_{k+1}|y_k)}{dq_{k|k+1}^{\theta}(y_k|y_{k+1})} = -\frac{1}{2} \|\epsilon_k\|^2 + \frac{1}{2} \left\| -\epsilon_k - \sqrt{w_k^{\text{CMCD}}} u_k \right\|^2$$

□
Based on Lemma 22, we simply deduce the expression of the discrete time versions of LV-based and KL-based CMCD losses.

**Corollary 23 (Expression of LV-CMCD loss).** When using the EM scheme to integrate SDE (23), the LV-based discrete time CMCD loss (40) simplifies as

$$ \mathcal{L}_{LV}^{\text{CMCD}}(\theta) = \text{Var} \left[ \log \frac{\pi^{\text{base}}(Y_0)}{\gamma(Y_K)} + \sum_{k=0}^{K-1} w_k^{\text{CMCD}} \left\{ \frac{1}{2} \|u_k^{\theta}\|^2 + (h_{t_k}^{\hat{\theta}}(Y_k) - h_{t_k}^{\theta}(Y_k))^{\top} u_k^{\theta} \right\} \right. \\ \left. + \sum_{k=0}^{K-1} \sqrt{w_k^{\text{CMCD}}} (u_k^{\theta})^{\top} Z_k \right], $$

$$ \text{with } u_k^{\theta} = \nabla \log \pi_{t_k}(Y_k) + \nabla \log \pi_{t_{k+1}}(Y_{k+1}) + h_{t_k}^{\theta}(Y_k) - h_{t_{k+1}}^{\theta}(Y_{k+1}), $$

where $\{Z_k\}_{k=0}^{K-1}$ are independently distributed according to $N(0, I_d)$, and $\{Y_k\}_{k=0}^K$ is recursively defined by $Y_0 \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$,

$$ Y_{k+1} = Y_k + \frac{\sigma^2 \delta_k}{2} \{ \nabla \log \pi_{t_k}(Y_k) + h_{t_k}^{\hat{\theta}}(Y_k) \} + \sigma \sqrt{\delta_k} Z_k. $$

**Proof.** This result is an application from Lemma 22, where we have $\epsilon_k = Z_k + \sqrt{w_k^{\text{CMCD}}} \{h_{t_k}^{\hat{\theta}}(Y_k) - h_{t_k}^{\theta}(Y_k)\}$ for any $k \in \{0, \dots, K-1\}$. $\square$

**Corollary 24 (Expression of KL-CMCD loss).** When using the EM scheme to integrate SDE (23), up to additional constants independent of $\theta$, the KL-based discrete time CMCD loss (39) simplifies as

$$ \mathcal{L}_{KL}^{\text{CMCD}}(\theta) = \mathbb{E} \left[ \log \frac{\pi^{\text{base}}(Y_0^{\theta})}{\gamma(Y_K^{\theta})} + \frac{1}{2} \sum_{k=0}^{K-1} w_k^{\text{CMCD}} \|u_k^{\theta}\|^2 \right], $$

$$ \text{with } u_k^{\theta} = \nabla \log \pi_{t_k}(Y_k^{\theta}) + \nabla \log \pi_{t_{k+1}}(Y_{k+1}^{\theta}) + h_{t_k}^{\theta}(Y_k^{\theta}) - h_{t_{k+1}}^{\theta}(Y_{k+1}^{\theta}), $$

where $\{Z_k\}_{k=0}^{K-1}$ are independently distributed according to $N(0, I_d)$, and $\{Y_k^{\theta}\}_{k=0}^K$ is recursively defined by $Y_0^{\theta} \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$,

$$ Y_{k+1}^{\theta} = Y_k^{\theta} + \frac{\sigma^2 \delta_k}{2} \{ \nabla \log \pi_{t_k}(Y_k^{\theta}) + h_{t_k}^{\theta}(Y_k^{\theta}) \} + \sigma \sqrt{\delta_k} Z_k. $$

## E ANNEALED MCMC SAMPLERS

In this section, we present a line of Monte Carlo algorithms that aim at sampling from a sequence of distributions $\{p_k\}_{k=0}^K \subset \mathcal{P}(\mathbb{R}^d)$, with $K \ge 1$, defined such that $p_0$ is an easy-to-sample distribution and $p_k$ is harder and harder to sample as $k$ increases. In this section, we assume that each distribution $p_k$ is absolutely continuous with respect to the Lebesgue measure, with tractable unnormalized density and tractable score, where both can be evaluated pointwise.

**Reminders on the Metropolis-Adjusted Langevin Algorithm (MALA).** This sampling algorithm is the only one in this section to be designed so as to sample from a single distribution $\pi \in \mathcal{P}(\mathbb{R}^d)$ with density $\gamma: \mathbb{R}^d \to \mathbb{R}_+$ known up to a normalizing constant. It is an extension of the largely used Unadjusted Langevin Algorithm (ULA) (Roberts & Tweedie, 1996).

Given a number of steps $L \ge 1$, a step-size $\lambda > 0$ and an easy-to-sample distribution $p^{\text{init}} \in \mathcal{P}(\mathbb{R}^d)$, ULA builds a Markov chain $(X^\ell)_{\ell=0}^L$ defined by $X^0 \sim p^{\text{init}}$ and for any $\ell \in \{0, \dots, L-1\}$,

$$ X^{\ell+1} = X^\ell + \lambda \nabla \log \gamma(X^\ell) + \sqrt{2\lambda Z^\ell}, \quad Z^\ell \sim N(0, I_d), $$

which reduces to a time discretization of the Langevin diffusion that admits $\pi$ as invariant distribution. Building upon this recursion, MALA (Roberts & Tweedie, 1996) addresses the discretization bias by adding a Metropolis-Hastings acceptance/rejection step at each iteration, see Algorithm 5. In practice, one can easily adapt the step-size $\lambda$ by targeting a specific acceptance ratio. We follow this paradigm in our numerical experiments, see Appendix H for more details.
**Algorithm 5: Metropolis-Adjusted Langevin Algorithm (MALA)**

**Input:** Target distribution $\pi$ with unnormalized density $\gamma$, step size $\lambda > 0$, number of steps $L \ge 1$, initial distribution $p^{\text{init}}$

▷ Initialization

$$X^0 \sim p^{\text{init}}$$

for $\ell \in \{0, \dots, L-1\}$ do

▷ Propose a new sample

$$\tilde{X}^{\ell+1} = X^{\ell} + \lambda \nabla \log \gamma(X^{\ell}) + \sqrt{2\lambda Z^{\ell}}, \quad Z^{\ell} \sim N(0, I_d)$$

▷ Compute the acceptance ratio

$$\alpha = \min \left( 1, \frac{\gamma(\tilde{X}^{\ell+1})q(X^{\ell}|\tilde{X}^{\ell+1})}{\gamma(X^{\ell})q(\tilde{X}^{\ell+1}|X^{\ell})} \right) \text{ where } q(y|x) = N(y; x + \lambda \nabla \log \gamma(x), 2\lambda I_d)$$

▷ Accept or reject depending on $\alpha$

$$U \sim U([0, 1])$$

if $U \le \alpha$ then

$$\mid X^{\ell+1} = \tilde{X}^{\ell+1}$$

else

$$\lfloor X^{\ell+1} = X^{\ell}$$

**Output:** Sequence of samples $X^{0:L}$ such that $X^{0:L} \stackrel{\text{iid}}{\sim} \pi$ (approximately)

**Annealed Langevin MCMC (AL-MCMC).** Based on exact samples from $p_0$, the Annealed Langevin MCMC algorithm (Song & Ermon, 2019; 2020) consists in sequentially sampling from $p_{k+1}$ with MALA (Algorithm 5) initialized in the last samples obtained at step $k$. The core idea behind this procedure is that the warm initialization obtained by the easy-to-sample distributions (i.e., when $k$ is low) can help to overcome sampling issues, such as high energy barriers, in the more complex distributions (i.e., when $k$ is large). We refer to Algorithm 6 for the pseudo-code of this approach.

**Algorithm 6: Annealed Langevin MCMC (AL-MCMC) algorithm**

**Input:** Sequence of target distributions $\{p_k\}_{k=0}^K$, step-sizes $\{\lambda_k\}_{k=1}^K \in (0, \infty)^K$, number of MALA steps per level $L \ge 1$

▷ Initialization

$$X_0^{0:L} \stackrel{\text{iid}}{\sim} p_0$$

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

▷ Sample from $p_{k+1}$ by starting at samples obtained from $p_k$

$$X_{k+1}^{0:L} = \text{MALA}(p_{k+1}, \lambda_{k+1}, L, \delta_{X_k^L}), \text{ see Algorithm 5}$$

**Output:** Sequence of samples $X_{0:L}^{0:L}$ such that $X_k^{0:L} \stackrel{\text{iid}}{\sim} p_k$ (approximately)

**Sequential Monte Carlo (SMC) and extension.** The Sequential Monte Carlo (Neal, 2001; Del Moral et al., 2006) algorithm lies on the same paradigm as Algorithm 6 but proceeds in parallel over $N$ chains (which are referred as particles) while ensuring that the warm initialization is correct. To do so, SMC relies on an additional resampling step that may occur before certain MALA run, based on accumulated importance weights that we next define.

Let $N \ge 1$, $k \in \{0, \dots, K-1\}$, $(\tilde{w}_k^n)_{n=1}^N$ be unnormalized importance weights accumulated on each chain up to step $k$ $(\tilde{w}_k^n = 1/N$ if $k=0$), and $(x_k^n)_{n=1}^N$ be samples from $p_k$ obtained independently as the last samples of $N$ parallel MALA runs (or exactly obtained if $k=0$). Then, the SMC importance weights designed to sample from $p_{k+1}$, denoted by $(\tilde{w}_{k+1}^n)_{n=1}^N$, and their renormalized versions, denoted by $(w_{k+1}^n)_{n=1}^N$, are defined as

$$\tilde{w}_{k+1}^n = \frac{p_{k+1}(x_k^n)}{p_k(x_k^n)} \tilde{w}_k^n, \quad w_{k+1}^n = \frac{\tilde{w}_{k+1}^n}{\sum_{j=1}^N \tilde{w}_{k+1}^j}.$$

The SMC resampling step then consists in obtaining $N$ new samples $(y_{k+1}^n)_{n=1}^N$ by random selection among the original samples $(x_k^n)_{n=1}^N$ via a multinomial distribution (with replacement) defined by the normalized weights $(w_{k+1}^n)_{n=1}^N$. These novel samples will then be used as starting points for
running MALA on the target $p_{k+1}$. In practice, this resampling step only occurs adaptively based on the value of the accumulated Effective Sample Size (ESS) defined at step $k$ by

$$ESS_k = \frac{\left(\sum_{n=1}^{N} \tilde{w}_k^n\right)^2}{\sum_{n=1}^{N} (\tilde{w}_k^n)^2}.$$

We detail the whole SMC procedure in Algorithm 7.

**Algorithm 7: Sequential Monte Carlo (SMC) algorithm with adaptive resampling**

**Input:** Sequence of target distributions {$p_k$}$_{k=0}^K$, step-sizes {$\lambda_k$}$_{k=1}^K \in (0, \infty)^K$, number of MCMC steps per level $L \ge 1$, number of particles $N \ge 1$, ESS resampling threshold $\alpha \in (0, 1]$

▷ Initialization

$$X_0^{1:N} \stackrel{\text{iid}}{\sim} p_0, \quad \tilde{W}_0^{1:N} = 1/N$$

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

▷ 1. Compute the SMC importance weights (in parallel, for $n = 1, \dots, N$)

$$\tilde{W}_{k+1}^n = (p_{k+1}(X_k^n)/p_k(X_k^n)) \tilde{W}_{k+1}^n$$

▷ 2. Compute the accumulated ESS (in parallel, for $n = 1, \dots, N$)

$$ESS_{k+1} = \left(\sum_{n=1}^{N} \tilde{W}_k^n\right)^2 / \left(\sum_{n=1}^{N} (\tilde{W}_k^n)^2\right)$$

▷ 3. Resample $X_k^{1:N}$ if the ESS too low

if $ESS_{k+1} < \alpha N$ then

▷ Normalize the importance weights (in parallel, for $n = 1, \dots, N$)

$$W_{k+1}^n = \tilde{W}_{k+1}^n / \sum_{j=1}^N \tilde{W}_{k+1}^j$$

▷ Resample the particles

$$I^{1:N} \sim \mathcal{M}(\tilde{W}_{k+1}^1, \dots, \tilde{W}_{k+1}^N)$$

$$Y_{k+1}^{1:N} = X_k^{I^{1:N}}, \quad \tilde{W}_k^{1:N} = 1/N$$

else

$$\lfloor Y_{k+1}^n = X_k^n$$

▷ 4. Sample from $p_{k+1}$ by starting from sample $Y_{k+1}^n$ (in parallel, for $n = 1, \dots, N$)

$$\tilde{X}_{k+1}^{1:L} = \text{MALA}(p_{k+1}, \lambda_{k+1}, L, \delta_{Y_{k+1}^n}), \text{ see Algorithm 5}$$

$$X_{k+1}^n = \tilde{X}_{k+1}^L$$

**Output:** Sequence of samples $X_{0:K}^{1:N}$ such that $X_k^{1:N} \stackrel{\text{iid}}{\sim} p_k$ (approximately)

Recently, Phillips et al. (2024) proposed an extension of the SMC algorithm, *Particle Denoising Diffusion Sampler* (PDDS), where they consider the specific case where the intermediate distributions {$p_k$}$_{k=0}^K$ are defined as the approximate marginals of a denoising diffusion process. Below, we explain the fundamentals of PDDS.

Let $\pi \in \mathcal{P}(\mathbb{R}^d)$ be a hard-to-sample distribution. With the same notations as in Section 2.1, consider the noising process $(X^*)_{t \in [0, T]}$ induced by the following SDE

$$dX_t^* = f(t)X_t^*dt + \sqrt{\beta(t)}dW_t, \quad X_0^* \sim \pi. \qquad (43)$$

Then, assuming that $X_T^* \sim \pi^{\text{base}}$, the time-reversal of SDE (43) is given by

$$dY_t^* = -f(T-t)Y_t^*dt + \beta(T-t)\nabla \log p_{T-t}^*(Y_t^*)dt + \sqrt{\beta(t)}dB_t, \quad Y_0^* \sim \pi^{\text{base}}, \qquad (44)$$

where $p_t^*$ is the density of the marginal distribution associated to $X_t^*$, which is intractable in general. Assume that we are given a sequence of intermediate unnormalized densities $(p_t)_{t \in [0, T]}$ that approximate the target densities $(p_t^*)_{t \in [0, T]}$.

In particular, the time-reversed SDE (44) can be approximated by

$$dY_t = -f(T-t)Y_t dt + \beta(T-t)\nabla \log p_{T-t}(Y_t)dt + \sqrt{\beta(T-t)}dB_t, \quad Y_0 \sim \pi^{\text{base}}. \qquad (45)$$

Consider a time discretization {$t_k$}$_{k=0}^K$ such that $t_0 = 0$ and $t_K = T$, and the sequence of target distributions {$p_k$}$_{k=0}^K$ defined by $p_k = p_{T-t_k}$, such that $p_k$ is expected to be harder and harder to sample from as $k$ increases.
Let $k \in \{0, \dots, K-1\}$, define $\delta_k = t_{k+1} - t_k$. Since SDE (43) is linear, it can be exactly integrated on time interval $[T - t_{k+1}, T - t_k]$, see Lemma 2. We denote by $p_{k|k+1}^*$ the corresponding transition density. On the other hand, SDE (45) can only be approximately integrated on time interval $[t_k, t_{k+1}]$ (due to the score term in the drift). For instance, this can be done using the Euler-Maruyama transition kernel<fnref n="4" />

$$p_{k+1|k}(\cdot|y_k) = N(y_k - f(T - t_k)\delta_k y_k + \beta(T - t_k)\delta_k \nabla \log p_{T-t_k}(y_k), \beta(T - t_k)\delta_k I_d).$$

To sample sequentially from $\{p_k\}_{k=0}^K$, Phillips et al. (2024) suggest to apply the SMC procedure detailed in Algorithm 7, where the k-th resampling step (1-3) is now defined as follows.

0. Sample $\hat{X}_{k+1}^n \sim p_{k+1|k}(\cdot|X_k^n)$ (in parallel, for $n = 1, \dots, N$);

1. Compute the importance weights (in parallel, for $n = 1, \dots, N$)

$$\tilde{W}_{k+1}^n = \frac{p_{k+1}(\hat{X}_{k+1}^n) p_{k|k+1}^*(X_k^n | \hat{X}_{k+1}^n)}{p_k(X_k^n) p_{k+1|k}(\hat{X}_{k+1}^n | X_k^n)} \tilde{W}_k^n;$$

2. Apply the resampling procedure (based on the ESS value) to the particles $\hat{X}_{k+1}^{1:N}$ with importance weights given by $\{W_{k+1}^n\}_{n=1}^N$.

Hence, this novel resampling step amounts to reweight the samples obtained from the approximate time-reversed SDE (45) by using the exact noising scheme as proposal.

**Replica Exchange (RE).** The Replica Exchange algorithm (Swendsen & Wang, 1986) aims at sampling from the annealed distributions $\{p_k\}_{k=0}^K$ with $K + 1$ parallel Markov chains, i.e., the k-th chain targets $p_k$. For the sake of pedagogy, we assume here that $K$ is even.

At each level, RE sampling is done via a local MCMC sampler, such as MALA, which is expected to have poor performance for large $k$. To alleviate this issue, every L MALA steps, RE randomly performs a swapping between Markov chains. More specifically, each consecutive pair $(k, k+1)$ of Markov chains (i.e., either $(k, k+1) \in \{(0, 1), \dots, (K-2, K-1)\}$ or $(k, k+1) \in \{(1, 2), \dots, (K-1, K)\}$), with respective current states $X_k$ and $X_{k+1}$, is swapped with probability

$$p_{k,k+1} = \min \left( 1, \frac{p_k(X_{k+1}) p_{k+1}(X_k)}{p_k(X_k) p_{k+1}(X_{k+1})} \right).$$

We provide the pseudo-code of this algorithm in Algorithm 9.

**Algorithm 8: Swapping step of Replica Exchange**

**Input:** Sequence of target distributions $\{p_k\}_{k=0}^K$ with even $K$, current samples $X_{1:K}$, swapping state $s \in \{0, 1\}$

**if** $s = 0$ **then**

Indexes = $\{0, 2, \dots, K-2\}$

**else**

Indexes = $\{1, 3, \dots, K-1\}$

▷ Compute the swapping probability between $k$ and $k+1$ (in parallel, for $k \in Indexes$)

$$p_{k,k+1} = \min \left( 1, \frac{p_k(X_{k+1}) p_{k+1}(X_k)}{p_k(X_k) p_{k+1}(X_{k+1})} \right)$$

▷ Swap depending on $p_{k,k+1}$ (in parallel, for $k \in Indexes$)

$U \sim U([0, 1])$

**if** $U < p_{k,k+1}$ **then**

$$\lfloor (X_k, X_{k+1}) = (X_{k+1}, X_k) \rfloor$$

**Output:** Randomly swapped sequence $X_{0:K}$

<footnote n="4">⁴In the original paper, the authors consider the VP noising scheme and compute $p_{k+1|k}$ as the transition kernel obtained by Exponential Integration of SDE (45) on $[t_k, t_{k+1}]$.</footnote>
**Algorithm 9: Replica Exchange (RE) algorithm**

**Input:** Sequence of target distributions {$p_k$}$_{k=0}^K$ with even $K$, step-sizes {$\lambda_k$}$_{k=1}^K \in (0, \infty)^K$, total number of MCMC steps $M \ge 1$, swap frequency $S \ge 1$

▷ Initialization

$X_{0:K}^0$ iid $\sim p_0$

$s = 0$

for $m = 0, \dots, M-1$ do

if $(m+1) \bmod S == 0$ then

▷ Swap Markov chains for indexes determined by $s$

$X_{0:K}^{m+1} = \text{swapRE}(X_{0:K}^m, s)$, see Algorithm 8

$s = (s+1) \bmod 2$

else

▷ Sample locally (in parallel, for $k = 0, \dots, K$)

$X_k^{m+1} = \text{MALA}(p_k, \lambda_k, 1, \delta_{X_k^m})$, see Algorithm 5

**Output:** Sequence of samples $X_{0:K}^{0:M}$ such that $X_k^{0:M}$ iid $\sim p_k$ (approximately)

**Known limitations of annealed MCMC samplers.** In general, annealed MCMC samplers usually target a sequence of distributions {$p_k$}$_{k=0}^K$, which forms a geometric bridge between an easy-to-sample distribution $\pi^{\text{base}}$ and a hard-to-sample distribution $\pi$, i.e., $p_k \propto (\pi^{\text{base}})^{(K-k)/K} \pi^{k/K}$. However, this scheme suffer from mode-switching for multi-modal distribution $\pi$ (Phillips et al., 2024), i.e., the mode weights of the intermediate distributions may vary a lot as $k$ increases. Therefore, transitions between annealing levels is hard in practice, as samples can easily get stuck in a high-probability modes. Moreover, due to the use of reweighting schemes, SMC and RE are highly sensitive to the overlap between consecutive distributions, which requires to take $K$ large in complex multi-modal settings, hence increasing the complexity of the methods.

# F ENERGY-BASED MODELS

This section is dedicated to describing Energy-Based Models.

## F.1 SINGLE-LEVEL EBMS

Let $p \in \mathcal{P}(\mathbb{R}^d)$. We assume having access to samples from $p$ but not to its density. Energy-based Models (EBM) perform a density estimation task by modeling the density of $p$ with $p^\varphi(x) = \exp(-E^\varphi(x))/Z^\varphi$, where $E^\varphi: \mathbb{R}^d \to \mathbb{R}$ is a neural network and $Z^\varphi = \int_{\mathbb{R}^d} \exp(-E^\varphi(x))dx$ is an unknown normalizing constant. This model is trained by maximizing the log-likelihood of the model $\mathcal{L}(\varphi) = \mathbb{E}[\log p^\varphi(X)]$ where $X \sim p$, whose gradients can be computed as

$$ \nabla_{\varphi}\mathcal{L}(\varphi) = \mathbb{E}[\nabla_{\varphi}E^{\varphi}(X^{-})] - \mathbb{E}[\nabla_{\varphi}E^{\varphi}(X^{+})], \quad X^{-} \sim p^{\varphi}, \quad X^{+} \sim p. \quad (46) $$

Given samples from $p$ (referred to as positive samples) and samples from $p^\varphi$ (referred to as negative samples), one can compute a Monte-Carlo estimator of $\nabla_{\varphi}\mathcal{L}(\varphi)$ using (46). In practice, the negative samples are usually obtained by running a MCMC sampler on $p^\varphi$ (which is possible because it is known up to a normalizing constant). However, obtaining truthful negative samples is the main challenge of the EBM training procedure as $p^\varphi$ is expected to be as multi-modal as $p$.

## F.2 MULTI-LEVEL EBMS

Let $p \in \mathcal{P}^{(K+1)}$, with $K \ge 1$. For $k = 1, \dots, K$, we assume having access to samples from $p_k$ but not to its density. Regarding $p_0$, we assume it to be known entirely. Our goal is to approximate the density of $p_k$ at every level $k \in \{1, \dots, K\}$. A naive approach would consist in defining single-level EBMS for each $k$. Assuming that the distributions $p_k$ are not completely unrelated, doing this would be however very inefficient. This section presents alternatives to this approach, which are all implemented in our code.
## F.2.1 PRIOR WORKS

**Diffusion Recovery Likelihood (DRL).** Recently, Gao et al. (2021); Zhu et al. (2024) proposed the DRL framework, based on the extra assumption that for any $k \in \{0, \dots, K-1\}$, the conditional distribution $p_{k|k+1}$ is known and given by $p_{k|k+1}(\cdot | x_{k+1}) = N(\alpha_{k+1}x_{k+1}, \sigma_{k+1}^2 I_d)$ for some $\alpha_{k+1} \in \mathbb{R}$ and $\sigma_{k+1} > 0$. This is typically the case when the joint distribution $p$ is defined via a discrete time approximation of a linear SDE, see Appendix C. In the following, we denote $Y_k = \alpha_k X_k$ where $X_k \sim p_k$. In DRL, the authors suggest to use a multi-level EBM $p_k^\varphi(\cdot) = \exp(-E^\varphi(k, \cdot))/Z_k^\varphi$, where $E^\varphi: \{1, \dots, K\} \times \mathbb{R}^d \to \mathbb{R}$ is a neural network and the normalizing constant $Z_k^\varphi = \int_{\mathbb{R}^d} \exp(-E^\varphi(k, x))dx$ is intractable. In this case, the density of $p_k$ is approximated, up to a multiplicative constant, by $\exp(-E^\varphi(k, \cdot))$. Using Bayes rule, they define for any $k \in \{0, \dots, K-1\}$ the following conditional EBM

$$p_{k+1|k}^{\varphi}(y_{k+1} | x_k) \propto p_{k|k+1}(x_k | y_{k+1}/\alpha_{k+1})p_{k+1}^{\varphi}(y_{k+1}) \\ = (Z_k^{\varphi}(k+1, x_k))^{-1} \exp \left( -\frac{1}{2\sigma_{k+1}^2} \|x_k - y_{k+1}\|^2 - E^\varphi(k+1, y_{k+1}) \right),$$

where $Z_k^\varphi(k+1, x_k) = \int_{\mathbb{R}^d} \exp(-(2\sigma_{k+1}^2)^{-1} \|x_k - y\|^2 - E^\varphi(k+1, y)) dy$. Note that if $\sigma_{k+1}^2$ is small enough, for instance when $K$ is large, then the quadratic term dominates in the expression of $\log p_{k+1|k}^\varphi(\cdot | x_k)$, which means that this conditional distribution localizes on $x_k$. Therefore, by applying a first order Taylor expansion of the energy term at $x_k$, the following Gaussian approximation holds

$$p_{k+1|k}^{\varphi}(y_{k+1} | x_k) \approx N(x_k - \sigma_{k+1}^2 \nabla E^\varphi(k+1, x_k), \sigma_{k+1}^2 I_d).$$

Hence, if $\sigma_{k+1}^2$ is chosen small enough, this conditional distribution will be very easy to sample from. Based on this observation, the authors suggest to use the following maximum likelihood objective

$$\mathcal{L}^{\text{DRL}}(\varphi) = \sum_{k=0}^{K-1} \mathcal{L}_k^{\text{DRL}}(\varphi), \quad \text{with} \quad \mathcal{L}_k^{\text{DRL}}(\varphi) = \mathbb{E} \left[ \log p_{k+1|k}^{\varphi}(\alpha_{k+1} X_{k+1} | X_k) \right], \quad (47)$$

where $X_{k+1} \sim p_{k+1}$ and $X_k \sim p_{k|k+1}(\cdot | X_{k+1})$. In particular, the gradient of the single-level DRL loss defined in (47) can be expressed as

$$\nabla_{\varphi} \mathcal{L}_{k}^{\text{DRL}}(\varphi) = \mathbb{E}[\nabla_{\varphi} E^{\varphi}(k, Y_{k}^{-})] - \mathbb{E}[\nabla_{\varphi} E^{\varphi}(k, \alpha_{k} X_{k}^{+})],$$

where $X_k \sim p_k$, $Y_{k+1}^{-} \sim p_{k+1|k}^{\varphi}(\cdot | X_k)$ and $X_{k+1}^{+} \sim p_{k+1}$. Up to a multiplicative constant, the density of $p_k$ can then be approximated by $x_k \mapsto p_k^{\varphi}(\alpha_k x_k)$.

**Algorithm 10: Gibbs-within-Langevin MCMC transition kernel**

**Input:** Previous state $(k, x) \in \{0, \dots, K\} \times \mathbb{R}^d$, step-size $\lambda > 0$

▷ Sample $\tilde{x}$ given $k$ with standard MCMC

$\tilde{x} = \text{MALA}(p^\varphi(\cdot | k), \lambda)$ where $p^\varphi(\cdot | k) \propto p_k^\varphi(\cdot)$

▷ Sample $\tilde{k}$ given $\tilde{x}$ with a multinomial distribution

$\tilde{k} \sim \mathcal{M}(\exp(-E^\varphi(0, \tilde{x})), \dots, \exp(-E^\varphi(K, \tilde{x})))$

**Output:** Next state $(\tilde{k}, \tilde{x})$

**Diffusion Assisted EBM (DA-EBM).** On the other hand, (Zhang et al., 2023) suggest to solve the problem by jointly modeling indexes and states. Here, the multi-level EBM is expressed as $p_k^\varphi(x) = \exp(-E^\varphi(k, x))/Z_k^\varphi$, where $E^\varphi: \{0, \dots, K\} \times \mathbb{R}^d \to \mathbb{R}$ is a neural network and the normalizing constant $Z_k^\varphi = \sum_{k=0}^K \int_{\mathbb{R}^d} \exp(-E^\varphi(k, x))dx$, which can be learn by maximizing the following maximum likelihood objective

$$\mathcal{L}^{\text{DAEBM}}(\varphi) = \mathbb{E}[\log p_k^{\varphi}(X_k)], \quad k \sim \text{U}(\{0, \dots, K\}), \quad X_k \sim p_k. \quad (48)$$

The gradient of objective (48) can be written as

$$\nabla_{\varphi} \mathcal{L}^{\text{DAEBM}}(\varphi) = \mathbb{E}[\nabla_{\varphi} E^{\varphi}(k^{-}, X^{-})] - \mathbb{E}[\nabla_{\varphi} E^{\varphi}(k^{+}, X_{k}^{+})],$$

where $(k^{-}, X^{-}) \sim p^{\varphi}, k^{+} \sim \text{U}(\{0, \dots, K\})$ and $X_{k}^{+} \sim p_{k+}$. In the same fashion as Kim & Ye (2023), the authors suggest to perform the negative sampling from $p^{\varphi}$ by doing a Gibbs-within-Langevin procedure, whose transition kernel is summarized in Algorithm 10.
## F.2.2 OUR APPROACH: USING AN ANNEALED MCMC SAMPLER AS BACKBONE

In this paper, we also suggest to learn a multi-level EBM approximating $p$, defined by $p_k^\varphi(x) = \exp(-E^\varphi(k, x)) / Z_k^\varphi$, where $E^\varphi: \{0, \dots, K\} \times \mathbb{R}^d \to \mathbb{R}$ is a neural network and the normalizing constant $Z_k^\varphi = \int_{\mathbb{R}^d} \exp(-E^\varphi(k, x)) dx$. We propose to learn $\varphi$ by maximizing a simple yet novel maximum likelihood joint objective

$$ \mathcal{L}_{\text{ours}}(\varphi) = \sum_{k=0}^{K} \mathcal{L}_{k}^{\text{ours}}(\varphi), \quad \mathcal{L}_{k}^{\text{ours}}(\varphi) = \mathbb{E}[\log p_{k}^{\varphi}(X_{k})], \quad X_{k} \sim p_{k}, \qquad (49) $$

with single-level gradient given by

$$ \nabla_{\varphi} \mathcal{L}_{k}^{\text{ours}}(\varphi) = \mathbb{E}[\nabla_{\varphi} E^{\varphi}(k, X_{k}^{-})] - \mathbb{E}[\nabla_{\varphi} E^{\varphi}(k, X_{k}^{+})], \quad X_{k}^{-} \sim p_{k}^{\varphi}, \quad X_{k}^{+} \sim p_{k}. $$

Unlike previous multi-level EBM algorithms, the negative sampling phase can be simply done in our framework by leveraging annealed MCMC algorithms presented in Appendix E on the sequence of EBM densities $\{p^\varphi(k, \cdot)\}_{k=0}^K$. This allows us to kill two birds in one stone as (i) we get negative samples for each single-level EBM and (ii) we overcome the negative sampling limitations of each EBM by leveraging the correlation between the consecutive levels. In practice, we use Replica Exchange as default annealed MCMC sampler because of its massively parallel capabilities. We summarize our training procedure in Algorithm 11. For ease of reading, we presented the continuous time analog of this approach in the main part of the paper, see Section 3.3.

### Algorithm 11: Multi-level EBM training using annealed MCMC samplers

**Input:** Number of training iterations $N \ge 1$, initial parameter $\varphi_0$, batch size $B \ge 1$, number of MCMC steps $L \ge 1$, predefined annealed MCMC sampler

**for** $n = 0, \dots, N-1$ **do**

▷ Get $B$ samples from $p_k$ (in parallel, for $k = 0, \dots, K$)

$\{X_{k,b}^{+}\}_{b=1}^{B}$ iid

▷ Get $B$ samples from each $p_k^{\varphi_n}$ using the annealed MCMC sampler

$\{X_{k,b}^{-}\}_{b=1,k=0}^{B,K} = \text{annealed MCMC}(\{p_k^{\varphi_n}\}_{k=0}^K, B)$

▷ Apply a stochastic gradient descent on $\varphi_n$

Compute the MC estimator $\hat{g}_n$ of the gradient of the loss $\mathcal{L}_{\text{ours}}(\varphi_n)$ defined in (49)

$$ \hat{g}_n = B^{-1} \sum_{k=0}^{K} \left( \sum_{b=1}^{B} \nabla_{\varphi} E^{\varphi_n} (k, X_{k,b}^{-}) - \sum_{b=1}^{B} \nabla_{\varphi} E^{\varphi_n} (k, X_{k,b}^{+}) \right) $$

Update $\varphi_n$ to $\varphi_{n+1}$ with Adam optimizer based on $\hat{g}_n$

**Output:** Optimized parameter $\varphi_N$

### Multi-level EBM parameterization as GMM tilting.
To reduce the computational footprint of multi-level EBMs in the context of LRDS, we suggest to define our EBM model as a tilting of a learned GMM. Let $\bar{\varphi} = \{w_j, m_j, \Sigma_j\}_{j=1}^J$ be the parameters learned by GMM-LRDS (Algorithm 3) to approximate $\hat{\pi}^{\text{ref}}$ and denote by $p_t^{\bar{\varphi}}$ the t-th marginal density of the resulting reference process. Building on this, we define $p_t^\varphi$ as a tilting of $p_t^{\bar{\varphi}}$ by

$$ \log p_t^\varphi(x) = \log p_t^{\bar{\varphi}}(x) - \text{NN}^\varphi(t, x), $$

where $\text{NN}^\varphi: [0, T] \times \mathbb{R}^d \to \mathbb{R}$ is a neural network with parameter $\varphi$. In this scenario, we design the initial $\varphi_0$ to ensure that $\text{NN}^\varphi_0 = 0$. In practice, in the case where the GMM shares very few similarities with $\hat{\pi}^{\text{ref}}$ close to $t = 0$, we rather consider

$$ \log p_t^\varphi(x) = \mathbf{1}_{t>t_{\text{lim}}}(t) \times \log p_t^{\bar{\varphi}}(x) - \text{NN}^\varphi(t, x), $$

where $t_{\text{lim}} \in [0, T]$. For most of our numerical experiments, we found that $t^{\text{lim}} = 0.2$ brought satisfying results.
# G ADDITIONAL METRICS FOR VARIATIONAL DIFFUSION-BASED METHODS

When specifically using annealed VI methods, one may consider additional metrics to evaluate the performance of the corresponding samplers. For completeness purpose, we provide in this section the expressions of the sampling metrics presented by Blessing et al. (2024), for all discrete time variational diffusion-based methods presented in this paper: RDS (including PIS and DDS), DIS and CMCD. All of these metrics are implemented in our code.

Given a fixed level of time discretization $K \ge 1$, these performance criteria require the evaluation of a variational importance weight function $w: (\mathbb{R}^d)^{K+1} \to \mathbb{R}_+$, assessing the quality of the variational approximation, that is specific to each variational setting.

* For RDS (including PIS and DDS), we define

$$w(y_{0:K}) = \frac{\mathcal{Z}^{\text{ref}}}{\mathcal{Z}} \frac{\mathrm{d}p_{0:K}^{\theta}(y_{0:K})}{\mathrm{d}p_{0:K}^{\star}(y_{0:K})} = \frac{\mathrm{d}p_{0:K}^{\theta}(y_{0:K})\gamma^{\text{ref}}(y_{K})}{\mathrm{d}p_{0:K}^{\text{ref}}(y_{0:K})\gamma(y_{K})},$$

where $p^{\theta} \in \mathcal{P}^{(K+1)}$ and $p^{\star} \in \mathcal{P}^{(K+1)}$ are respectively defined as discrete time approximations of the path measures $\mathbb{P}^{\theta}$ and $(\mathbb{P}^{\text{ref}})^R$; see Appendix D.3 for more details.

* For DIS, we define

$$w(y_{0:K}) = \frac{1}{\mathcal{Z}} \frac{\mathrm{d}p_{0:K}^{\theta}(y_{0:K})}{\mathrm{d}p_{0:K}^{\star}(y_{0:K})},$$

where $p^{\theta} \in \mathcal{P}^{(K+1)}$ and $p^{\star} \in \mathcal{P}^{(K+1)}$ are respectively defined as discrete time approximations of the path measures $\mathbb{P}^{\theta}$, induced by SDE (21), and $(\mathbb{P}^{\star})^R$; see Appendix D.4 for more details.

* For CMCD, we define

$$w(y_{0:K}) = \frac{1}{\mathcal{Z}} \frac{\mathrm{d}p_{0:K}^{\theta}(y_{0:K})}{\mathrm{d}q_{0:K}^{\theta}(y_{0:K})}$$

where $p^{\theta} \in \mathcal{P}^{(K+1)}$ and $q^{\theta} \in \mathcal{P}^{(K+1)}$ are respectively defined as discrete time approximations of the path measures $\mathbb{P}^{\theta}$, induced by SDE (23), and $(\mathbb{Q}^{\theta})^R$, where $\mathbb{Q}^{\theta}$ is induced by SDE (22); see Appendix D.5 for more details.

Note that, for each variational setting, $w$ does not depend on the normalization constant of the target distribution: in particular, $\mathcal{Z}$ simplifies in the DIS and CMCD expressions. This also stands for $\mathcal{Z}^{\text{ref}}$ in the RDS setting.

## G.1 'REVERSE' PERFORMANCE CRITERIA

**Definition.** Standard variational metrics, referred to as 'reverse' performance criteria by Blessing et al. (2024), are defined as expectations with respect to the variational distribution.

(a) The Evidence Lower Bound (ELBO), expected to be maximized, defined by

$$\text{ELBO} = -\mathbb{E}[\log w(Y_{0:K}^{\theta})], \quad Y_{0:K}^{\theta} \sim p_{0:K}^{\theta}.$$

For any RDS setting with intractable $\mathcal{Z}^{\text{ref}}$, this criterion may assess the performance of the RDS sampler itself, but cannot be used for numerical comparison with different variational settings.

(b) The MC estimation of the normalizing constant $\hat{Z}_r$, expected to be close to $\mathcal{Z}$, defined by

$$\hat{Z}_r = \mathbb{E}[w(Y_{0:K}^{\theta})], \quad Y_{0:K}^{\theta} \sim p_{0:K}^{\theta}.$$

For any RDS setting with intractable $\mathcal{Z}^{\text{ref}}$, this criterion cannot be used.

(c) The normalized Explained Sum of Squares (nESS$_r$), expected to be close to 1, defined by

$$\text{nESS}_r = \frac{\mathbb{E}[w(Y_{0:K}^{\theta})]^2}{\mathbb{E}[w(Y_{0:K}^{\theta})^2]}, \quad Y_{0:K}^{\theta} \sim p_{0:K}^{\theta}.$$
**Computation.** Based on the results from Appendix D.2, we detail the simulation of $Y_{0:K}^{\theta}$ and the computation of $\log w(Y_{0:K}^{\theta})$ for each setting.

* **RDS setting:** $Y_0^{\theta} \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$,

$$Y_{k+1}^{\theta} = a_k Y_k^{\theta} + b_k \{s_{T-t_k}^{\text{ref}}(Y_k^{\theta}) + g_{T-t_k}^{\theta}(Y_k^{\theta})\} + \sqrt{c_k} Z_k, \quad Z_k \sim \text{N}(0, I_d),$$

and the importance weight function is given by

$$\log w(Y_{0:K}^{\theta}) = \sum_{k=0}^{K-1} \frac{w_k}{2} \|g_{T-t_k}^{\theta}(Y_k^{\theta})\|^2 + \sum_{k=0}^{K-1} \sqrt{w_k} g_{T-t_k}^{\theta}(Y_k^{\theta})^{\top} Z_k + \log \frac{\gamma^{\text{ref}}}{\gamma}(Y_K^{\theta}),$$

where $\{w_k, a_k, b_k, c_k\}_{k=0}^{K-1}$ depend on the noising scheme and the discretization setting.

* **DIS setting:** $Y_0^{\theta} \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$,

$$Y_{k+1}^{\theta} = a_k Y_k^{\theta} + b_k s_{T-t_k}^{\theta}(Y_k^{\theta}) + \sqrt{c_k} Z_k, \quad Z_k \sim \text{N}(0, I_d)$$

and the importance weight function is given by

$$\log w(Y_{0:K}^{\theta}) = \sum_{k=0}^{K-1} \frac{w_k}{2} \|s_{T-t_k}^{\theta}(Y_k^{\theta})\|^2 + \sum_{k=0}^{K-1} \sqrt{w_k} s_{T-t_k}^{\theta}(Y_k^{\theta})^{\top} Z_k + \log \frac{\pi^{\text{base}}(Y_0^{\theta})}{\gamma(Y_K^{\theta})},$$

where $\{w_k, a_k, b_k, c_k\}_{k=0}^{K-1}$ come from the EI discretization of the VP noising scheme.

* **CMCD setting:** $Y_0^{\theta} \sim \pi^{\text{base}}$, and for any $k \in \{0, \dots, K-1\}$,

$$Y_{k+1}^{\theta} = Y_k^{\theta} + \frac{\sigma^2 \delta_k}{2} \{\nabla \log \pi_{t_k}(Y_k^{\theta}) + h_{t_k}^{\theta}(Y_k^{\theta})\} + \sigma \sqrt{\delta_k} Z_k, \quad Z_k \sim \text{N}(0, I_d),$$

and the importance weight function is given by

$$\log w(Y_{0:K}^{\theta}) = \frac{1}{2} \sum_{k=0}^{K-1} w_k^{\text{CMCD}} \|u_k^{\theta}\|^2 + \sum_{k=0}^{K-1} \sqrt{w_k^{\text{CMCD}}} (u_k^{\theta})^{\top} Z_k + \log \frac{\pi^{\text{base}}(Y_0^{\theta})}{\gamma(Y_K^{\theta})},$$
with $u_k^{\theta} = \nabla \log \pi_{t_k}(Y_k^{\theta}) + \nabla \log \pi_{t_{k+1}}(Y_{k+1}^{\theta}) + h_{t_k}^{\theta}(Y_k^{\theta}) - h_{t_{k+1}}^{\theta}(Y_{k+1}^{\theta}).$$

## G.2 'FORWARD' PERFORMANCE CRITERIA

**Definition.** As observed by Blessing et al. (2024), the 'reverse' variational metrics presented above may not be able to quantify mode collapse. Hence, the authors propose novel variational metrics, referred to as 'forward' performance criteria, that are defined as expectations with respect to the ground truth distribution.

(a) The Evidence Upper Bound (EUBO), expected to be minimized, defined by

$$\text{EUBO} = \mathbb{E}[\log w(Y_{0:K})], \quad Y_{0:K} \sim p_{0:K}^{*}.$$

For any RDS setting with intractable $\mathcal{Z}^{\text{ref}}$, this criterion may assess the performance of the RDS sampler itself, but cannot be used for numerical comparison with different variational settings.

(b) The MC estimation of the normalizing constant $\hat{Z}_f$, expected to be close to $\mathcal{Z}$, defined by

$$\hat{Z}_f = (\mathbb{E}[w(Y_{0:K})^{-1}])^{-1}, \quad Y_{0:K} \sim p_{0:K}^{*}.$$

For any RDS setting with intractable $\mathcal{Z}^{\text{ref}}$, this criterion cannot be used.

(c) The normalized Explained Sum of Squares (nESS$_f$), expected to be close to 1, defined by

$$\text{nESS}_f = (\mathbb{E}[w(Y_{0:K})^{-1}]\mathbb{E}[w(Y_{0:K})])^{-1}, \quad Y_{0:K} \sim p_{0:K}^{*}.$$
**Computation.** Based on the results from Appendix D.2, we detail the simulation of $Y_{0:K}$ and the computation of $\log w(Y_{0:K})$ for each setting.

* **RDS setting with EM scheme:** $Y_K \sim \pi$, and for any $k \in \{0, \dots, K-1\}$,

$$Y_k = S_{k+1}(t_k)Y_{k+1} + S_{k+1}(t_k)\sigma_{k+1}(t_k)Z_k, \quad Z_k \sim N(0, I_d),$$

where $S_k(t) = \exp\left(\int_{T-t_k}^{T-t} f(u)du\right)$, $\sigma_k^2(t) = \int_{T-t_k}^{T-t} \frac{\beta(u)}{S_k(u)^2}du$, and the importance weight function is given by

$$\log w(Y_{0:K}) = -\sum_{k=0}^{K-1} w_k^{\text{EM}} g_{T-t_k}^{\theta}(Y_k)^{\top} \left\{ s_{T-t_k}^{\text{ref}}(Y_k) + \frac{1}{2} g_{T-t_k}^{\theta}(Y_k) \right\} \\ -\sum_{k=0}^{K-1} \sigma_{k+1}(t_k) g_{T-t_k}^{\theta}(Y_k)^{\top} Z_k + \sum_{k=0}^{K-1} \left( \frac{1}{S_{k+1}(t_k)} - 1 + f(T-t_k)\delta_k \right) g_{T-t_k}^{\theta}(Y_k)^{\top} Y_k \\ + \log \frac{\gamma^{\text{ref}}}{\gamma}(Y_K),$$

as a consequence of Lemma 10.

* **RDS setting with EI scheme (PBM):** $Y_K \sim \pi$, and for any $k \in \{0, \dots, K-1\}$,

$$Y_k = \frac{Y_{k+1}}{a_k^{\text{PBM}}} + \frac{\sqrt{c_k^{\text{PBM}}}}{a_k^{\text{PBM}}} Z_k, \quad Z_k \sim N(0, I_d),$$

and the importance weight function is given by

$$\log w(Y_{0:K}) = -\sum_{k=0}^{K-1} w_k^{\text{PBM}} g_{T-t_k}^{\theta}(Y_k)^{\top} \{s_{T-t_k}^{\text{ref}}(Y_k) + \frac{1}{2} g_{T-t_k}^{\theta}(Y_k)\} \\ -\sum_{k=0}^{K-1} \sqrt{w_k^{\text{PBM}} g_{T-t_k}^{\theta}(Y_k)^{\top} Z_k} + \log \frac{\gamma^{\text{ref}}}{\gamma}(Y_K),$$

as a consequence of Lemma 13.

* **RDS setting with EI scheme (VP):** $Y_K \sim \pi$, and for any $k \in \{0, \dots, K-1\}$,

$$Y_k = \frac{Y_{k+1}}{a_k^{\text{VP}}} + \frac{\sqrt{c_k^{\text{VP}}}}{a_k^{\text{VP}}} Z_k, \quad Z_k \sim N(0, I_d),$$

and the importance weight function is given by

$$\log w(Y_{0:K}) = -\sum_{k=0}^{K-1} w_k^{\text{VP}} g_{T-t_k}^{\theta}(Y_k)^{\top} \{s_{T-t_k}^{\text{ref}}(Y_k) + \frac{1}{2} g_{T-t_k}^{\theta}(Y_k)\} \\ -\sum_{k=0}^{K-1} \sqrt{w_k^{\text{VP}} g_{T-t_k}^{\theta}(Y_k)^{\top} Z_k} + \log \frac{\gamma^{\text{ref}}}{\gamma}(Y_K),$$

as a consequence of Lemma 16.

* **DIS setting:** $Y_K \sim \pi$, and for any $k \in \{0, \dots, K-1\}$,

$$Y_k = \frac{Y_{k+1}}{a_k^{\text{VP}}} + \frac{\sqrt{c_k^{\text{VP}}}}{a_k^{\text{VP}}} Z_k, \quad Z_k \sim N(0, I_d),$$

and the importance weight function is given by

$$\log w(Y_{0:K}) = -\frac{1}{2} \sum_{k=0}^{K-1} w_k^{\text{VP}} \|s_{T-t_k}^{\theta}(Y_k)\|^2 - \sum_{k=0}^{K-1} \sqrt{w_k^{\text{VP}} s_{T-t_k}^{\theta}(Y_k)^{\top} Z_k} + \log \frac{\pi^{\text{base}}(Y_0)}{\gamma(Y_K)},$$

as a consequence of Lemma 19.
• CMCD setting: $Y_K \sim \pi$, and for any $k \in \{0, \dots, K-1\}$,

$$Y_k = Y_{k+1} + \frac{\sigma^2 \delta_k}{2} \{ \nabla \log \pi_{t_{k+1}}(Y_{k+1}) - h_{t_{k+1}}^{\theta}(Y_{k+1}) \} + \sigma \sqrt{\delta_k} Z_k, \quad Z_k \sim N(0, I_d),$$

and the importance weight function is given by

$$\log w(Y_{0:K}) = -\frac{1}{2} \sum_{k=0}^{K-1} w_k^{\text{CMCD}} \|u_k^{\theta}\|^2 - \sum_{k=0}^{K-1} \sqrt{w_k^{\text{CMCD}}} (u_k^{\theta})^{\top} Z_k + \log \frac{\pi^{\text{base}}(Y_0)}{\gamma(Y_K)} \\ \text{with } u_k^{\theta} = \nabla \log \pi_{t_k}(Y_k) + \nabla \log \pi_{t_{k+1}}(Y_{k+1}) + h_{t_k}^{\theta}(Y_k) - h_{t_{k+1}}^{\theta}(Y_{k+1}),$$

as a consequence of Lemma 22.

# H IMPLEMENTATION DETAILS

## H.1 DETAILS ON TARGET DISTRIBUTIONS

**Gaussian mixture from Figure 2.** The distribution from Figure 2 is a 16-dimensional mixture between $N(\mathbf{m}_1, \Sigma_1)$ and $N(\mathbf{m}_2, \Sigma_2)$ with weights 0.7 and 0.3 respectively. Let $\mathbf{1}_d$ be the $d$-dimensional vector with all components equal to 1, then we have

$$\mathbf{m}_1 = -0.6 \times \mathbf{1}_{16}, \quad \mathbf{m}_2 = +0.6 \times \mathbf{1}_{16}.$$

Moreover, denote $R_{\theta}$ a rotation matrix with angle $\theta$ between the first and last axes, then

$$\Sigma_1 = R_{\pi/4} \text{diag}(10^{-2}, \dots, 10^{-2}, 10^{-1}) R_{\pi/4}^T, \quad \Sigma_2 = R_{\pi/6} \text{diag}(10^{-1}, 10^{-2}, \dots, 10^{-2}) R_{\pi/6}^T.$$

In practice, we estimate the weight of the first mode $w_1$ by computing a Monte Carlo estimator of

$$\int \mathbf{1}_{N(x; \mathbf{m}_1, \Sigma_1) > N(x; \mathbf{m}_2, \Sigma_2)}(x) \pi(x) dx. \qquad (50)$$

**Bi-modal Gaussian mixture.** This Gaussian mixture is used in Figure 1 and also in Section 5. It is a $d$-dimensional mixture between $N(-\mathbf{1}_d, \Sigma)$ and $N(\mathbf{1}_d, \Sigma)$ with weights 2/3 and 1/3 respectively. The different values of $\Sigma$ are recapped in Table 3. Just like the previous experiment, we also use the formula (50) to estimate the strongest mode weight with Monte Carlo approach.

**Multi-modal Gaussian mixtures.** Given a number of mixture components $L > 2$, we define a $d$-dimensional Gaussian mixture $x \in \mathbb{R}^d \mapsto \sum_{\ell=1}^L w_\ell N(x; \mathbf{m}_\ell, 0.5 I_d)$, where the mean locations $\{\mathbf{m}_\ell\}_{\ell=1}^L$ are independently distributed according to $U([-L, L]^d)$ and the mixture weights $\{w_\ell\}_{\ell=1}^L$ are strictly increasing, defined with a constant geometrical increment such that $w_L/w_1 = 3$. To assess the mode weight recovery in this multi-modal setting, we compute the total variation distance between the exact mode weight histogram and the mode weight histogram computed from Monte Carlo approximation.

**Rings distribution.** Let $T: \mathbb{R}^+ \times [0, 2\pi] \to \mathbb{R}^2$ be the transformation from polar to cartesian coordinates. Let $\theta \sim U([0, 2\pi])$ and $R \sim N(r, \sigma^2)$ with $r, \sigma > 0$, we say that $T(R, \theta)$ is distributed as a ring of radius $r$ and width $\sigma$. The Rings distribution is a mixture of rings with radiuses $r = 1, 3, 5$ and width $\sigma = 0.1$. The weights of each ring within the mixture are 1/9, 3/9 and 5/9. In this specific case, we use more MCMC chains and GMM components than the number of modes. The initial points of the MCMC samplers are fetched by drawing the same amount of uniformly distributed points on each ring.

**Checkerboard distribution.** We divide the square $[-4, 4]^2$ into 16 squares of size $2 \times 2$. This distribution is defined as a mixture of uniform distributions on those squares in an interleaving fashion. We assign a weight of 18.75% to each 4 squares on the left and 6.25% to each 4 squares on the right. The initial points of the MCMC samplers are taken as the middle of each square.
<table><caption>**Table 3: Covariance settings for the bi-modal Gaussian mixture.** We denote by logdiag(a, b) with $0 < a < b$ the diagonal matrix in $\mathbb{R}^{d \times d}$ whose diagonal is a log-linear interpolation between a and b. Moreover, $Q_d \in \mathbb{R}^{d \times d}$ is a random orthogonal matrix built by doing the QR decomposition of a random matrix distributed according to $U([0, 5]^{d \times d})$. The seed used to generate this last matrix is fixed at 42.</caption><thead><tr><th>Type</th><th>Difficulty</th><th>Cond. number</th><th>Covariance</th></tr></thead><tbody><tr><td>Diagonal</td><td>Isotropic</td><td>1</td><td>(5 × 10<sup>-2</sup>)<sup>2</sup> × I<sub>d</sub></td></tr><tr><td>Diagonal</td><td>Medium</td><td>10<sup>2</sup></td><td>(5 × 10<sup>-2</sup>)<sup>2</sup> × logdiag(10<sup>-2</sup>, 1)</td></tr><tr><td>Diagonal</td><td>Hard</td><td>10<sup>4</sup></td><td>(5 × 10<sup>-2</sup>)<sup>2</sup> × logdiag(10<sup>-4</sup>, 1)</td></tr><tr><td>Full</td><td>Medium</td><td>10<sup>2</sup></td><td>(5 × 10<sup>-2</sup>)<sup>2</sup> × Q<sub>d</sub> × logdiag(10<sup>-2</sup>, 1) × Q<sub>d</sub><sup>T</sup></td></tr><tr><td>Full</td><td>Hard</td><td>10<sup>4</sup></td><td>(5 × 10<sup>-2</sup>)<sup>2</sup> × Q<sub>d</sub> × logdiag(10<sup>-4</sup>, 1) × Q<sub>d</sub><sup>T</sup></td></tr></tbody></table>

**The $\phi^4$ field system.** The $\phi^4$ model is a simplified framework often used as a continuous version of the Ising model, aiding in the exploration of phase transitions within statistical mechanics. As per Gabrié et al. (2022), we focus on a discretized version of this model, set on a one-dimensional grid with size $d = 32$. Each configuration is represented as a $d$-dimensional vector $(\phi_i)_{i=1}^d$. To further constrain the system, we fix the field to zero at both ends by setting $\phi_0 = \phi_{d+1} = 0$.

The negative log-density of the distribution is then expressed as follows

$$ \ln \pi_h(\phi) = -\beta \left( \frac{ad}{2} \sum_{i=1}^{d+1} (\phi_i - \phi_{i-1})^2 + \frac{1}{4ad} \sum_{i=1}^{d} (1 - \phi_i^2)^2 + h\phi_i \right). \quad (51) $$

We selected parameter values that make the system bimodal, setting $a = 0.1$ and the inverse temperature $\beta = 20$, while adjusting the value of $h$. We define $w_+$ as the statistical frequency of configurations where $\phi_{d/2} > 0$, and $w_-$ as the frequency where $\phi_{d/2} < 0$. When $h = 0$, the system is symmetric under the transformation $\phi \to -\phi$, so we anticipate $w_+ = w_-$. However, for $h > 0$, the negative mode becomes more prevalent.

For large values of $d$, the relative likelihood of the modes can be approximated using Laplace expansions at the 0th and 2nd orders. Letting $\phi_+^h$ and $\phi_-^h$ represent the local maxima of (51), these approximations provide the following results respectively

$$ \frac{w_{-}}{w_{+}} \approx \frac{\pi_{h}(\phi_{-}^{h})}{\pi_{h}(\phi_{+}^{h})}, \quad \frac{w_{-}}{w_{+}} \approx \frac{\pi_{h}(\phi_{-}^{h}) \times |\det H_{h}(\phi_{-}^{h})|^{-1/2}}{\pi_{h}(\phi_{+}^{h}) \times |\det H_{h}(\phi_{+}^{h})|^{-1/2}}, \quad (52) $$

where $H_h$ is the Hessian of the function $\phi \to \ln \pi_h(\phi)$. We compute the modes $\phi_+^h$ and $\phi_-^h$ using a gradient descent on the distribution's potential. Since we do not have access to ground truth samples in practice, we compare the Laplace approximations defined in (52), considered as our ground truth, to a Monte Carlo estimation of $w_- / w_+$.

**Bayesian logistic regression models.** Finally, we evaluate the performance of a Bayesian logistic model, defined for any pair $(x, y) \in \mathbb{R}^{\text{dim}} \times \{0, 1\}$ by the likelihood $p(y|x; w, b) = \text{Bernoulli}(y; \sigma(w^T x + b))$, where $w \in \mathbb{R}^{\text{dim}}$ is a weight vector, $b \in \mathbb{R}$ is an intercept and $\sigma$ is the sigmoid function. Given a dataset $\mathcal{D} = \{(x_j, y_j)\}_{j=1}^M \subset \mathbb{R}^{\text{dim}} \times \{0, 1\}$ of size $M$, we aim to sample from the posterior distribution $p(w, b|\mathcal{D}) \propto p(w, b) \prod_{j=1}^M p(y_j|x_j; w, b)$ where $p(w, b) = p(w)p(b)$ is a fixed prior distribution. Following Blessing et al. (2024), we consider four real-world settings of binary classification problem: **Ionosphere** (dim = 35, $M = 351$), **Sonar** (dim = 61, $M = 208$), **German Credit** (dim = 25, $M = 1000$), **Breast Cancer** (dim = 31, $M = 569$). Each of these datasets is randomly split into a training subset $\mathcal{D}_{\text{train}}$ of size $0.8M$ and a test subset $\mathcal{D}_{\text{test}}$ of size $0.2M$. In this setting, we define the target distribution $\pi = p(w, b|\mathcal{D}_{\text{train}})$ and evaluate the sampling quality by computing the (unnormalized) average predictive log-likelihood $\log p(w, b|\mathcal{D}_{\text{test}})$, which is expected to be maximized.
<table><caption>Table 4: Complexity of each algorithm during the overall sampling procedure. We track the number of target evaluations and neural network evaluations. We consider that computing $\pi$ is as expensive as computing $\nabla \log \pi$. We assume an infinite parallel computational capabilities. Note that the cost of the training procedure can be amortized to sample multiple times.</caption><thead><tr><th>Method</th><th>Number of target evaluations</th><th>Number of neural network evaluations</th></tr></thead><tbody><tr><td>SMC</td><td>O(KM)</td><td>0</td></tr><tr><td>RE</td><td>O(M)</td><td>0</td></tr><tr><td>PIS/DDS/DIS/CMCD</td><td>O(NK)</td><td>O(NK)</td></tr><tr><td>PDDS</td><td>O(A(N + KM))</td><td>O(A(N + KM))</td></tr><tr><td>iDEM</td><td>O(AN)</td><td>O(ANK)</td></tr><tr><td>LRDS</td><td>O(N)</td><td>O(NK)</td></tr></tbody></table>

In each case, the prior is carefully chosen as a Gaussian distribution with the following parameters:

* Ionosphere: $p(w) = N(0, 5.25 I_{34})$ and $p(b) = N(4.25, 0.25^2)$,

* Sonar: $p(w) = N(0, 4.5 I_{60})$ and $p(b) = N(-2.5, 0.5^2)$,

* German Credit: $p(w) = N(0, 1.25 I_{24})$ and $p(b) = N(3.25, 0.5^5)$,

* Breast Cancer: $p(w) = N(0, 3.75 I_{30})$ and $p(b) = N(31, 2^2)$.

We draw the reader's attention to the fact that these posterior distributions don't exhibit explicit multi-modality features. Therefore, for EBM-LRDS, the tilting EBM parameterization is rather based on a Gaussian approximation of the target distribution than a GMM approximation. The initial point of the MCMC samplers are sampled from the prior distribution.

## H.2 COMPUTATIONAL COMPARISON OF THE PRESENTED METHODS

This section aims at comparing the computational budget of each algorithm. In particular, we track the number of neural network evaluations and the number of target evaluations. For each algorithm, we take the following notations.

* For SMC: $K$ (number of annealing steps), $M$ (number of MCMC steps per level)

* For RE: $M$ (total number of MCMC steps)

* For RDS (PIS/DDS/LRDS), DIS and CMCD: $N$ (number of training steps), $K$ (number of discretization steps)

* For PDDS: $N$ (number of training steps), $M$ (number of MCMC steps per level), $K$ (number of discretization steps), $A$ (number of adaptation steps)

* For iDEM: $N$ (number of training steps), $K$ (number of discretization steps), $A$ (number of adaptation steps)

We detail the complexities of each algorithm in Table 4. Note that, in LRDS, the complexity of MCMC sampling (to compute $\hat{\pi}^{\text{ref}}$) or the complexity of learning the reference process were ignored as they are negligible compared to the training and sampling budgets; see Appendix H.3 for details.

## H.3 HYPER-PARAMETERS OF EACH SAMPLING ALGORITHM

In this section, we detail every hyper-parameter for each algorithm involved in the computations of the results in Section 5. The computationally-related parameters were chosen to ensure a comparable execution clock-time between the different algorithms on the same hardware.

**Construction of the MCMC dataset for $\hat{\pi}^{\text{ref}}$.** We build the dataset of reference samples by running MCMC samplers initialized in the modes of the target distribution. We run 4 chains per mode. We use the Metropolis-Adjusted Langevin algorithm (MALA), see Algorithm 5, as default sampler, except for the Checkerboard distribution where we use the Random Walk Metropolis Hastings (RWMH) algorithm as the score is constant. In both cases, we adapt the step size automatically for 8192 warmup steps to aim at a 70% acceptance rate. The datasets are 60000 samples long.
**Annealed MCMC methods (SMC, RE).** For both algorithms, we incorporate prior knowledge into the base distribution by taking $\pi^{\text{base}} = N(\mathbf{m}, \Sigma)$ where $\mathbf{m}$ and $\Sigma$ are computed using maximum likelihood on samples from $\hat{\pi}^{\text{ref}}$. In each case, the annealed density path is a linearly-spaced geometric interpolation, where $\pi_k \propto (\pi^{\text{base}})^{(K-k)/K} \pi^{k/K}$. For SMC (see Algorithm 7), we use $K=128$ annealing levels with $N=1024$ particles and $L=64$ MCMC steps per level (with a 4096 steps warmup) and set the ESS resampling threshold as $\alpha=1$. For RE (see Algorithm 9), we use $K=128$ annealing levels with 256 independent chains per level and a swapping frequency every $S=8$ steps for a total of $M=32$ steps (with a 8192 steps warmup). Both algorithms use MALA as MCMC backbone with a step-size automatically tuned to achieve 70% acceptance rate.

**Variational diffusion-based methods (PIS, DDS, DIS, CMCD, LRDS).** For PIS and DDS we use the implementation provided by Richter et al. (2023), while we use the implementation provided by Blessing et al. (2024) for DIS and CMCD. For PIS, DDS and DIS, we set the hyper-parameter $\sigma$ as advised by Appendix I.2 by computing the mean and the variance of the samples from $\hat{\pi}^{\text{ref}}$ to add prior knowledge. For CMCD, we sample a version of $\pi$ which is shifted by $\mathbf{m}$, defined as the mean computed on samples from $\hat{\pi}^{\text{ref}}$, and we set $\sigma$ as the square root of the average variance estimated on those samples. For PIS, DIS and DDS, we use $K=100$ linearly spaced time discretization steps while CMCD learns its own time discretization grid of size $K=128$ with a final time set at $T=10$.

For all of them, we perform 4096 optimization steps with a batch of size 2048. The neural network at stake is a Fourier MLP, as in Zhang & Chen (2022), with 4 layers of width 64. In the case of PIS, DDS and DIS, we use a target-informed parameterization by adding $NN(t) \times \nabla \log \pi(x)$ (where NN is a time-dependent scalar neural network) to the Fourier MLP, as suggested by the respective authors. As recommended by Vargas et al. (2024), we do not consider this extra-parameterization in CMCD, since the drift of the generative process is already informed by $\pi$. We highlight that, by default, we do not use this target-informed parameterization in LRDS, since the reference process is specifically designed to be it-self target-informed, hence avoiding useless evaluations of the target score. Additionally, as recommended by Zhang & Chen (2022), we design the LRDS guidance network such that $g^{\theta_0} = 0$. This ensures that the very first sampling phase is driven solely using the reference process.

For all diffusion-based algorithms and target distributions, we ensured that the noising schedules were chosen to ensure that the target distribution gets successfully noised at time $T$. For LRDS, we use an exponential integration of the respective time-reversed SDE as default transition kernel; see Appendix D.3.2 and Appendix D.3.3.

**Adaptive diffusion-based methods (iDEM, PDDS).** For iDEM, we use the implementation provided by the original paper (Akhound-Sadegh et al., 2024). In order to provide prior information leveraging $\hat{\pi}^{\text{ref}}$, (i) we standardize the target distribution by using the empirical mean and variance of $\hat{\pi}^{\text{ref}}$ (ii) we preload the buffer with samples from $\hat{\pi}^{\text{ref}}$. Following the design choices of the original paper, we use a Variance Exploding noising scheme with a geometric variance schedule $\sigma_t$. Because we standardized the target distribution, we decided to take $\sigma_T = 5$ for each target distribution to ensure that the distribution is properly noised. We used the same hyper-parameters as in their Gaussian mixture experiment with 400 epochs instead of 1000 to ensure computational fairness. For PDDS, we also use the implementation provided by the original paper (Phillips et al., 2024). In order to provide prior information leveraging $\hat{\pi}^{\text{ref}}$, we standardized the target distribution by leveraging the empirical mean and variance of $\hat{\pi}^{\text{ref}}$ thus bypassing the original VI step of their algorithm. We also significantly increased the default computational budget by using 128 discretization steps, 2048 particles and training for 100000 steps with a batch size of 512. The rest of the hyper-parameters are the default ones.

**LRDS reference fitting details.** For GMM-LRDS, the EM algorithm is taken from Pedregosa et al. (2011). We use a diagonal covariance for Rings distribution, Checkerboard distribution, diagonal Gaussian mixtures and Bayesian posterior distributions, and use a full covariance otherwise. In the full covariance case, we regularize the covariance to ensure their positivity. For EBM-RDS, we fit multi-level EBMs with Algorithm 11 using the Replica Exchange annealed sampler (see Algorithm 9) as backbone. This sampler has the advantage of being massively parallel and thus suited for GPU computations. We perform 200 epochs with a batch size of 32 for each noise level. To increase gradient accuracy and enhance efficiency, we accumulate the gradients over 10 steps where
we keep the same negative samples but update the positive ones. We perform 32 MCMC steps to sample a batch of negative samples with a swap happening every 8 steps. We keep 16 MCMC steps out of 32 to compute the expectations. To compensate the short length of our MCMC chains, we leverage a persistent buffer to kickstart the chains at each noise level. Lastly, we utilise the GMM tilting EBM parameterization detailed at the end of Appendix F where the neural network parameterized as $x \mapsto NN(t, x)^T x$ as suggested by Salimans & Ho (2021). We leverage the GMM tilting initialization of the EBM to perform exact MCMC sampling at the very first gradient step. Moreover, we preconditioned the network as advised by Karras et al. (2022) by leveraging the Gaussian approximation of $\hat{\pi}^{\text{ref}}$.

# I ADDITIONAL RESULTS

## I.1 LINK BETWEEN THE TARGET AND THE REFERENCE PROCESSES

Let $\pi^{\text{ref}} \in \mathcal{P}(\mathbb{R}^d)$ be an arbitrary distribution. We recall that under mild assumptions on $\pi^{\text{ref}}$, see e.g., Cattiaux et al. (2023), the time-reversal of $\mathbb{P}^{\text{ref}}$ is associated to the SDE

$$dY_t^{\text{ref}} = -f(T-t)Y_t^{\text{ref}}dt + \beta(T-t)s_{T-t}^{\text{ref}}(Y_t^{\text{ref}})dt + \sqrt{\beta(T-t)}dB_t, \quad Y_0^{\text{ref}} \sim \mathbb{P}_T^{\text{ref}}. \quad (53)$$

**Link with Schrödinger Bridge.** Since $\mathbb{P}^*$ and $\mathbb{P}^{\text{ref}}$ are based on the same noising diffusion process, these path measures are linked by the relation $\mathbb{P}^* = \pi \otimes \mathbb{P}_0^{\text{ref}}$. Hence, using the KL chain rule (Léonard, 2014), it is clear that $\mathbb{P}^*$ is solution to the following minimization problem over path measure space

$$\arg \min \{\text{KL}(\mathbb{P} | \mathbb{P}^{\text{ref}}) : \mathbb{P} \in \mathcal{P}(C_T), P_0 = \pi\},$$

which is often referred to as a half-Schrödinger bridge problem in stochastic control literature.

**Link with Doob's h-transform.** This relation on path measure space can also be written as $\mathbb{P}^* = \frac{d\pi}{d\pi^{\text{ref}}} \cdot \mathbb{P}^{\text{ref}}$, see Appendix D.1. Therefore, solely based on the SDE (53) describing $(\mathbb{P}^{\text{ref}})^R$, $(\mathbb{P}^*)^R$ may be expressed as a Doob's h-transform of $(\mathbb{P}^{\text{ref}})^R$ via the SDE

$$dY_t = -f(T-t)Y_t dt + \beta(T-t)\{s_{T-t}^{\text{ref}}(Y_t) + h_{T-t}(Y_t)\}dt + \sqrt{\beta(T-t)}dB_t, Y_0 \sim \mathbb{P}_T^*,$$

where the Doob's control function $h: [0, T] \times \mathbb{R}^d \to \mathbb{R}^d$ is defined, for any $y_t \in \mathbb{R}^d$ and any $t \in [0, T]$, by

$$h_{T-t}(y_t) = \nabla \log \mathbb{E}_{(\mathbb{P}^{\text{ref}})_{T|t}^R} \left[ \frac{d\pi}{d\pi^{\text{ref}}}(Y_T) \mid Y_t = y_t \right] = \nabla \log \mathbb{E}_{(\mathbb{P}^{\text{ref}})_{0|T-t}^R} \left[ \frac{d\pi}{d\pi^{\text{ref}}}(X_0) \mid X_t = y_t \right].$$

Additionally, we have for any $x_t \in \mathbb{R}^d$,

$$\begin{aligned} \frac{p_t^*}{p_t^{\text{ref}}}(x_t) &= \frac{1}{p_t^{\text{ref}}(x_t)} \int_{\mathbb{R}^d} p_{t|0}^*(x_t|x_0) d\pi(x_0) \\ &= \int_{\mathbb{R}^d} \frac{d\pi}{d\pi^{\text{ref}}}(x_0) \frac{p_{t|0}^{\text{ref}}(x_t|x_0) d\pi^{\text{ref}}(x_0)}{p_t^{\text{ref}}(x_t)} \\ &= \mathbb{E}_{(\mathbb{P}^{\text{ref}})_{0|t}} \left[ \frac{d\pi}{d\pi^{\text{ref}}}(X_0) \mid X_t = x_t \right]. \end{aligned}$$

By combining previous computations, we obtain $h_t = \nabla \log p_t^*/p_t^{\text{ref}} = g_t$ for any $t \in [0, T]$.

## I.2 OPTIMAL SETTING OF ISOTROPIC GAUSSIAN REFERENCE DISTRIBUTION

The goal of this section is to explain how to set the hyperparameter $\sigma \in (0, \infty)$ in PIS and DDS settings so as to 'optimize' the reference distribution $\pi^{\text{ref}} = N(0, \sigma^2 I_d)$, when targeting a multi-modal distribution. The same reasoning can also be applied in the DIS setting to set the base distribution $\pi^{\text{base}}$, chosen as an isotropic Gaussian distribution.

Let $\pi \in \mathcal{P}(\mathbb{R}^d)$ be our target distribution. Assume that we are given a diagonal Gaussian approximation $\hat{\pi} = N(\mathbf{m}, \Sigma)$ of $\pi$, with mean $\mathbf{m} \in \mathbb{R}^d$ and covariance $\Sigma = \text{diag}(\Gamma^2) \in \mathbb{R}^{d \times d}$, where
$\Gamma \in (0, \infty)^d$. We propose to set $\pi^{\text{ref}}$ as 'close' as possible to $\tilde{\pi}$ in PIS and DDS by solving the following variational problem

$$ \arg \min \{\text{KL}(\tilde{\pi} | N(0, \sigma^2 I_d)) : \sigma \in (0, \infty)\}. $$

Let $\sigma > 0$. We have

$$ \begin{aligned} \text{KL}(\tilde{\pi} | N(0, \sigma^2 I_d)) &= \frac{1}{2} \left[ \log \frac{|\sigma^2 I_d|}{|\Sigma|} - d + \sigma^{-2} \mathbf{m}^T \mathbf{m} + \text{Tr}(\sigma^{-2} \Sigma) \right] \\ &= \frac{1}{2} \left[ d \log(\sigma^2) - \sum_{j=1}^{d} \log \Gamma_j^2 - d + \sigma^{-2} \mathbf{m}^T \mathbf{m} + \sigma^{-2} \sum_{j=1}^{d} \Gamma_j^2 \right]. \end{aligned} $$

Moreover, the derivative with respect to $\sigma^2$ is given by

$$ \frac{\partial}{\partial \sigma^2} \text{KL}(\tilde{\pi} | N(0, \sigma^2 I_d)) = \frac{1}{2} \left[ d\sigma^{-2} - \sigma^{-4} \left( \mathbf{m}^T \mathbf{m} + \sum_{j=1}^{d} \Gamma_j^2 \right) \right]. $$

In particular, the optimal solution of this variational problem is

$$ \sigma_{\pi}^{2} = d^{-1} \left( \mathbf{m}^{T} \mathbf{m} + \sum_{j=1}^{d} \Gamma_{j}^{2} \right), $$

which motivates us to set $\pi^{\text{ref}} = N(0, \sigma_{\pi}^2 I_d)$.

In practice, when we do not have samples from the target distribution, we rather rely on a diagonal Gaussian approximation of the empirical distribution $\hat{\pi}^{\text{ref}}$ (obtained with local MCMC samplers) and set $\pi^{\text{ref}} = N(0, \sigma_{\hat{\pi}^{\text{ref}}}^2 I_d)$.

**Application to Gaussian mixtures.** In the particular case where $\pi$ is a Gaussian mixture, we have an analytical formula for $\tilde{\pi}$.

Let $\pi = \sum_{j=1}^{J} w_j \mathcal{N}(\mathbf{m}_j, \Sigma_j)$ be our target distribution with $w_{1:J} \in \Delta_J$, where the j-th component has mean $\mathbf{m}_j \in \mathbb{R}^d$ and covariance $\Sigma_j = \text{diag}(\lambda_j) \in \mathbb{R}^{d \times d}$ with $\lambda_j \in (0, \infty)^d$. Consider $X \sim \pi$ and $Z \sim \mathcal{M}(w_1, \dots, w_J)$. We first have

$$ \mathbb{E}[X] = \mathbb{E}[\mathbb{E}[X | Z]] = \sum_{j=1}^{J} w_j \mathbf{m}_j . $$

We define the *diagonal variance* of $\pi$ by $\text{diagVar}(\pi) = (\text{Var}[X]_1, \dots, \text{Var}[X]_d) \in (0, \infty)^d$ where $\text{Var}[X]_i = \mathbb{E}[(X_i - \mathbb{E}[X]_i)^2]$ for any $i \in \{1, \dots, d\}$. By the law of total variance, we have

$$ \begin{aligned} \text{Var}[X]_i &= \mathbb{E}[\text{Var}[X | Z]_i] + \text{Var}[\mathbb{E}[X | Z]]_i \\ &= \sum_{j=1}^{J} w_j (\lambda_j^2)_i + \sum_{j=1}^{J} w_j ((\mathbf{m}_j)_i - \mathbb{E}[X]_i)^2 . \end{aligned} $$

By characterisation of Gaussian distributions, $N(\mathbb{E}[X], \text{diagVar}(\pi))$ is the closest diagonal Gaussian distribution to $\pi$ in the Kullback-Leibler sense. Therefore, in this particular setting, $\tilde{\pi}$ may be defined by $\mathbf{m} = \mathbb{E}[X]$ and $\Gamma^2 = \text{diagVar}(\pi)$. We use this approach when computing $\sigma_{\pi}^2$ in the numerical experiment commented in Figure 1.

## I.3 FAILURE OF LOCAL MCMC SAMPLERS ON MULTIMODAL DISTRIBUTIONS

As we show in Figure 7, *local MCMC samplers such as the Random Walk Metropolis Hastings (RWMH) algorithm (Metropolis et al., 1953), the Metropolis-adjusted Langevin Algorithm (MALA) (Roberts & Tweedie, 1996), the Hamiltonian Monte Carlo (HMC) algorithm (Duane et al., 1987; Brooks et al., 2011) or the No-U Turn Sampler (NUTS) (Hoffman & Gelman, 2014) tend to produce Markov chains that get trapped in modes. The resulting samples are thus not representative of the global landscape.*
<figure><img src="image_9.png" /><figcaption>Figure 7: Samples from different MCMC samplers when sampling from a bi-modal Gaussian mixture in 2 dimensions - There are 4 different chains (in different colors). The MCMC samplers ran for 4096 warmup steps before producing those 1024 samples. We only display 1 sample every 4 steps.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Increasingly expressive GMM approximations of the rings distribution - (Top) Approximation using diagonal covariances (Bottom) Approximation using full covariances</figcaption></figure>

## I.4 LIMITATIONS OF GAUSSIAN MIXTURE MODELS

Figure 8 shows the progressive GMM approximation of the Rings distribution by leveraging an increasing number of components. This figure shows that, in this case, getting a good approximation of the distribution requires significantly more components than the number of modes (here, 3). This highlights the limitation of the approximation power of Gaussian Mixture models.

## I.5 FURTHER EXPERIMENTAL RESULTS

In this section, we display detailed results for the experimental settings presented in Section 5, as well as additional results that assess the robustness and the superiority of LRDS compared to its competitors for a large diversity of settings. When ground truth samples are available (i.e., in all settings except $\phi^4$ and Bayesian regression), we display statistical estimations of the integral probability metrics defined in Appendix B, namely, the regularized Wasserstein distance with regularization $\varepsilon = 10^{-3}$ ($W_{2,\varepsilon}$), the Maximum Mean Discrepancy (MMD) and the sliced Kolmogorov-Smirnov distance (KS).

### I.5.1 HIGH-DIMENSIONAL BI-MODAL GAUSSIAN MIXTURES

Below, we display the results of all considered samplers, including LRDS, when targeting the bi-modal Gaussian mixtures defined in Appendix H.1 with increasing dimension $d \in \{8, 16, 32, 64\}$. To assess the sampling quality at a global level, we compute for each sampler the absolute error 'Mode Err.' when estimating the strongest mode weight (in %, comprised between 0%, the best, and 66.7%, the worst) via Monte Carlo approximation. All metrics are computed based on 8192 samples. In particular,