Learned Reference-based Diffusion Sampler for multi-modal distributions
Authors: Maxence Noble
Louis Grenioux
Marylou Gabrié
Alain Oliviero Durmus
Venue: ICLR 2025 Poster (2025)
Keywords: sampling from multi-modal densities
diffusion models
variational inference
stochastic optimal control
multi-level energy-based models
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We propose LRDS, a novel practical diffusion-based sampler for multi-modal distributions that leverages prior knowledge on the locations of the modes.

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.

Published as a conference paper at ICLR 2025
LEARNED REFERENCE-BASER DIFFUSION SAMPLING FOR MULTI-MODAL DISTRIBUTIONS
Maxence Noble", Louis Grenioux', Marylou Gabrié & Alain Oliviero Durmus
CMAP, CNRS
École polytechnique, Institut Polytechnique de Paris
91 120 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 ac- exact samples and relying solely evaluations of unnormalized densities. The resulting samplers approximate the time-reversal of a noising diffusion pro- cess, bridging the target distribution an easy-to-sample base distribution. In practice, the performance of these methods heavily depends on key hyperparam- eters 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 meth- ods. Building on existing approaches, we introduce Learned reference-basele Dif- fusion Sampler (LRDS), a methodology specifically designed to leverage prior knowledge 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 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.
INTRODUCTION
We consider the problem of sampling from a probability density known up to a normalizing constant. More precisely, we consider a target distribution æ E (RR)) with probability density given by I (x)/ where y Rd -> JR,, can be pointwise evaluated and the normalizing constant Z Jiz )x)) de is intractable. This problem appears in wide variety of applications such as Bayesian statistics (Liu & Liu| 200T Kroese et al. 20IT), statistical mechanics (Krauth) 2006) or molecular dynamics (Stoltz al. 2010). In particular,
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 notoriously hard problem, yet, maybe surprisingly, even the location of the modes is known, sampling л remains a very challenging problem (Noé et al.) 2019 Pompe al.) 2020 Grenioux et al. 2023). In this work, we aim address this specific issue and will assume that we have access to the location of the modes as prior information on π. However, do not assume have access a priori ground truth samples from π.
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)
guaranteed to be efficient for high-dimensional target distributions that are log-concave or satisfy functional incqualities (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 il- lustration. To alleviate this issue, annealing strategies have been proposed. They consist in bridging an easy-to- pusy-to--ampll sample base distribution hase to the target distribution π via intermediate distributions
Both authors contributed equally.
Corresponding authors: (maxence noble e-burri lot, louis, greni grenioux)3pp apol yytchhi wtechnique.ed. que

---

Published as a conference paper at ICLR 2025
This sequence should form smooth path in the distribution space, so that sampling from Tx+1 starting from samples from T5 is relatively simple. Typically, up to normalizing constants again, the densities of these annealed distributions are chosen as geometric interpolation between π base and y, and are sampled either sequentially (Neal 2001 Del Moral al., 2000) 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 al. 2009 Tawn et al. 2020 Syed et al., 2022 Hence, recovering the proportions of the different modes of IT is hard in practice.
{TT) k= with K 2 1, To -hhsse and TKK
Alternatively, the bridging distributions { 1k 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 at the foundation of diffusion-baase -haood 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 al. 2024). Following this approach, several diffusion- -aased samplers relying exclusively
on Monte Carlo methods have recently shown promises for sampling from multi-modas 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 popular alternative to MCMC methods. Given class of parameterized (variational) distributions that are easy to sample from, they aim to find the optimal parameters that minimize a fixed discrep- ancy metric with respect to π.. Standard VI settings involve mean-field approximations (Wainwright & Jordan 2008), mixture models (Arenz et al. 2023), and more recently normalizing flows (Rezende al. 2014]] Papamakarios et al. (021). that typically optimize reverse Kullback-Letis Leibler (KL) ob- jective. 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 ex- tended target distribution with marginals corresponding to the { This approach has been explored following explicit bridging paths, i.e., { }
admit tractable unnormalized densities
(Wu et al. 2020; Arbel et al. [2021] Geffner & Domke [202]] Matthews et al.) 2022] Doucet et al.) 2022] Gettner & Domkee 2203)] Vargas et al.) (024). but the resulting samplers may suffer from explained above. On the other hand, recent class of annealed VI samplers
mode switching,
using implicit biffusion-basee paths has emerged (Tzen & Raginsky [2019]
Holdijk et al., 2023
Pavon 2022 Zhang & Chen 2022, Berner et al.) 2223 Vargas et al. /2023 b Zhang & Chen 2022
Berner et al. 2023 Vargas et al.] 2023b [Akhound-Saderr
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- 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 ap- proaches: (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, more computa- tionally intensive scheme, takes advantage of Energy- -asee Models for harder sampling problems.
We show that GMM-LRDS and EBM.-LDS outperform previous diffusion-based samplers on challenging multi-modal settings.
Notation. For any measurable space (X, X)). we denote by P (X) the space of probability mea- sures defined on (X, X²). Unless specified, if X is a topological space, &' is defined as the Borel a-field of X. For any T > 0. we denote by C, C( T), RRd) the space of continuous functions from [0, T) to Rd endowed with the uniform topology; hence (CC,) is the set of continuous-time stochastic processes (or path measures) defined on J]. For any P (CT). we denote by pp its time-reversal, defined such that if (Xx)reto.r] ~ IP. then (XT-i)tejo.r) ~ 1. For any t E [0, 77, we denote by P, the marginal distribution of P at time The density of the Gaussian distribution with mean m E IRd and covariance Σ E ddd is denoted by x + N(x; m, Σ) and A refers to the Σ/= my = 1) where J > 1.
J-dimensional simplex, ie, A = ((w;)/=1 E [0,
2

---

Published as a conference paper at ICLR 2025
Finally, for ease of understanding, we will denote noising diffusion processes by (Xt) driven by Brownian motion (We))t(o.T) while (Yt))e(o.7) will refer to denoising diffusion processes, driven by Brownian motion (B0) In this paper, è 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) con-
tinuous 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, 7] by the linear SDE
dX, = f(t))xd BB(t)dW, Xo 1,~m π
(1)
where the horizon T > is fixed, (W.) [0, 2 + (0, xc) and / [0, + R. We denote by Px the path measure associated to the SDE
standard d-dimensional Brownian motion, ,
[ With appropriate choices of f and 3, this SDE admits explicit transition kernels and Xr is approximately, or exactly, distributed according to an easy-to-sample base distribution huse
In the case where 1 = hhuue,, 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 (209));
Zhang & Chen 2022 : Vargas et al. (2023b), where n base 6g. Otherwise, when we only have Fý. R
π hase we define this setting as an 'ergodic' noising scheme. A well-known example is the Variance Preserving (VP) noising process for which hase N(O, a² Ia) with a > 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 arbitrary noising scheme. We refer to Appendix more details.
Denote by p. the density of P w.r.t. the Lebesgue measure. Under mild assumptions on म, f and B, see e.g.. Cattiauu et al. (023). the time-reversal of the noising process PA i.e., the distribution of (Xx-<)<(oor)) and denoted by (P*) is associated to the SDE
(2) where (Be) is another standard d-dimensional Brownian motion. By definition of the time- reversal, holds that Yr и. Therefore, if we were able simulate this diffusion process, we would obtain approximate samples from π. We therefore refer to (Y))eejo,r)) as the target process. However, the scores (V log p. [e((T involved in the drift function of ( 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 201T) since samples from æ are not available 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 reference process as described next.
dY = - t) + 8(T - t)V log + /3(T- (jdB: Yo F7.
Variational reference- approach. This approach relies on reference process that is so- lution SDE ( with initial condition Xo ~ no where nd E ((Rd) chosen such that the marginal scores of the associated path measure pref are tractable. Similarly т,, we assume that the probability density of ee is known up a multiplicative constant and is given by x 7 (x)) z where - Rdd + R, and gre = J& re (x))r is not necessarily tractable. Finally, for any t [0. denote the marginal scores a V log pfs, with of being the density of Ippe w.r.. the Lebesgue measure, Based on the tractable reference scores, the SDE describing the target pro- cess (2) can be rewritten as
dY = - -iVddt B(T () (VI) + (Y.) + BTT Y P7.
(3)
where 9: V lo /pppppp"" is now the only intractable term. This formulation the time-reversed SDE has already been largely used in the diffusion model literature, especially for conditional gen- erative models; see e.g., Dhariwal & Nichol (202)) In this specific context, Si is called a guidance Doob-h transform control, see Ap-
term. From probabilistic perspective, 9r be interpreted pendix[I.I] Exploiting the reference-bart based formulation given by (3) Richter et al. (2023) propose to estimate the guidance terms (J)) rather than the target scores log )ee[o.T]-
3

---

Published as a conference paper at ICLR 2025
Table I: Connection to prior works. Here, E 10,00 1] x Idd + R refers to a neural network.
Method
Noising
(()
2
(x)
Reference parameters
PIS Zhang & Chen 2022
PBM(a)
o(x;0,"PT1s)
2 o//0 (T-(Ttt)
(tuned) (tuned)
DDS Vargas al. 2023a
VP(z)
N(r;0,0Fla)
x/0²
GMM-LRDS (Section
3.2
arbitrary
N,, 2n)
anulytical
{uy, my Ey1 (leurned) ip (leamne)
EBM-LRDS (Section
arbitrury
expl EP 0,))) unknown
(,,
To this end, the authors consider the class of variational path measures (IP") JAEO C (CT) where for any 0 E θ IPth is associated the diffusion process
dr - d) dt + 3(T - t) ())"" (X,")d -B Y.o
bise
(4)
where g" [0. 2 Rd is typically neural network. In particular, choosing n- Gaussian distribution recovers exactly the previous VI methods Puth Integral Sampler (PIS) (Zhang & Chen
2022) and Denoising Diffusion Sampler (DDS) (Vargas et al., 2023a), see Table [ Furthermore,
Richter et al. (2023) propose to minimize the following LV-based objective
(iv( Var ()(d))//((()))))))))) (.r)))
Y
~ pp
(5)
where 0 is a detached version of 0, meaning that gradients with respect to will not see ô. 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, [021] has the benefits improve mode exploration, to avoid costly computations of the loss gradients and to have zero variance at the optimal solution (Richter et al.j [003 Section 2.3). Finally, the LV loss (5) can be further explicited, detailed in the following proposition.
/ base and there exists 0 E such that of = 9t. Then, the loss defined (5) achieves optimal solution at H* and, setting @ = log(> /7), simplifies as
Proposition Assume that P1
Lix Var S/ B(T - t) 19 (x.))
de ["
VB(T 1797- (dB.("))BBBB dB, -((r))
(6)
2
This result is an adaptation of (Richter et al. pendix D.I for completeness. We also emphasize that the assumption made Proposition[]] I only needed for ergodic noising schemes. In practice, Richter et al. (2023) only consider two specific ver- sions of this continuous-time variational loss, where π is determined by the PIS and DDS settings. In the following, we will refer to these extensions of PIS and DDS LV-PIS and LV-DDS. In con- trast, we keep a general perspective, by letting -- completely arbitrary, and describe our proposition of practical implementation of a Reference-basee Diffusion Sampler (RDS) in this context.
2023 Lemma 3.1), which proof is restated in Ap-
Riccter
2.2 REFERENCE BASED DIFFUSION SAMPLING IN PRACTICE
Discrete-time setting. In practice, the parameterized process (Y, (ee cannot be simulated ex- actly and SDE can only be numerically solved with small step-size using, for example, the Euler-Maruyamp, (EM) or Exponential Integration (El) (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 (1x k such that to = 0, lx = T and K 2 1. Given an arbitrary reference distribution ndd, we propose to approximate the continuous-time objective defined in (6 by the following discrete time RDS objective
CRDS(0) (0) Var
K
wh (Y&) { of (Yt) Σ
1 (
K
Vwref- (Yx) Ze o(Yr)
k=0
where (-) are independently distributed according to N(O, {M k is recursively defined by Yo ~ π base and for any l: E (O...., .., K
1-1
-1),
Yx+1 axYi + bu buloff (Yx) + o (Yx) + Volls.
(8)
4

---

Published as a conference paper at ICLR 2025
with 0 being detached version of f and {wy, ag. bR, Ck k= being tractable coefficients that de- pend on the choice of discretization and noising schemes. In particular, if K is sufficiently large, {Ex) has approximately the same distribution as (Y: k==- We refer to Appendix D.2 for 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" exhibited by Richter et al. (2023) are only designed for LV-PIS and LV-DDS, and may be seen particular instances of
K-1
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 time objective ( for some N 2 1. Like other variational diffusion-based objectives, the RDS loss is simulation-based meaning that, the whole process { defined in (8) must be simulated each SGD step. The complete training procedure of RDS is summarized in Algorithmß] stated in Appendix [A] After training, the final parameter 0 is fixed, and approximate samples from ग are obtained by simulating {
from g", using (8) and taking YK-
LEARNED RDS FOR MULTI-MODAL DISTRIBUTIONS
In this section, we highlight the crucial role played by the reference distribution म- in RDS when targeting multi-modal distributions. We first show in Section 3.T the limitations of previous vari- ational reference-basee methods and give intuition on how प. should be chosen. Based on this observation, we introduce Learned Reference-based Diffiusion Sampler (LRDS), a novel diffusion- based sampler defined as a practical instance of RDS, where ाh is learned from prior knowledge on π,, namely the location of their modes. We present two approaches based on Gaussian Mixture Models (GMM-LRDS), see Section 3.2, and Energy Energy-Based Models (EBM-LRDS), see Section 3.3
3.1 ON THE IMPORTANCE OF THE REFERENCE DISTRIBUTION
Choosing a Gaussian distribution for गd ensures that the entire reference process is Gaussian. Therefore, the scores (xf*) [=e[o.T] can be analytically computed. We will refer to this case Gaus- sian 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 []. i.e., x- parameter is fixed before running the optimization over 0 and a priori disconnected from the target distribution. To assess the sensitivity of the algorithms with respect to the hyperparameter 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 D left). We observe that the performance is highly dependent on the value of o, with an optimal accuracy achieved when the isotropic variance is directly related to statisties of the target distribution; see Appendix I.2 for details. However, it is hard to estimate this quantity without access to ground truth samples.
where the hyperparameter is a E (0, oo). This
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 re as a Gaussian mixture with the same components as the target with their relative weights as hyperparameter (see Figure [ right). We find that this design of π re, 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 to adopt the following paradigm for RDS: given a multi-modal target distribution ग, we aim to design N" close to π, 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 learn well-suited reference distribution.
AAthough the authors do not provide any expression of discrete-ting version of 6 in their paper, discrete-time loss is implemented for LV-PIS and LV-DDS in their official codebase available at https: //gi wwtp:thww.o..ctuubbbe.eeeeee sde sampler
5

---

Published as a conference paper at ICLR 2025
Mode weight estimation error
0.4
LV-PIS
LV-DDS
PBM
0.3
VP
0.2
0.1
0.0
10
10 10 10
10'
0
10el
Figure Illustration of the decisive role of the reference distribution. Here, we target 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. (Leftt): Results for LV-PIS and LV-DDS when varying the value of their hyperparameter (which directly determines श as shown in Table D]. The green dotted line represents the optimal variance for Gaussian approximation m Appendix 1.2 for related computations. (Right) Results for RDS in PBM and VP settings when setting ग as a Gaussian mixture with the same modes as x, but w is replaced by W Details on the design of this experiment given in Appendix H.T
LRDS pipeline. Our algorithm LRDS proceeds in three main steps:
(a) Obtain reference samples. Given a standard MCMC sampler which targets æ, such as the Metropolis-Adjusee Langevin Algoiihm (MALA) (Roberts & Tweedie 1996), we simulate multi- ple Markov chains that are initialized in the target mode locations. We refer to the obtained samples as the reference samples and denote by and the corresponding empirical distribution. As recalled in Appendix I.3] in presence high energy barriers between the modes, æ fails at estimating the global energy landscape. Moreover, /re does not admit density w.r.t. the Lebesgue measure.
(b) Define the reference process. In this stage, we aim (i) learning re continuous approxi- K-1 mation of भt and (ii) computing the marginal scores times (T- -== of the induced reference process preí i.e., the noising process defined by SDE (T). where Xo T- In practice, set ge either as a Gaussian Mixture Model, which relies on light parameterization, or as an Eneegy-BBases 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, describe two practical implementations of stage (b) LRDS using different parameterizations of the reference distribution.
3.2 GAUSSIAN MIXTURE MODEL LRDS
We first propose set re as a Gaussian mixture to integrate multi-modality in RDS. In particular, we have z and any reference density is parameterized as ref(() Σラ= (-; N.mmm,,3.. E,) where J > 1 is the number of reference modes, {wj (wj)//=1 A, are the reference weights, m, Rd and 2. ddd are respectively the mean and the covariance of the j-th reference mode. In this set- ting, each marginal of [pee also a Gaussian mixture whose parameters
be simply deduced from 7 making their score fully tractable, see Appendix B Therefore, the only difficulty at stage (b) lies in the estimation of the reference parameters ;p wy m, E,>} for some vance. If > 1, we refer to this version of the algorithm as Gaussian Mixture Model LRDS (GMM- LRDS). this case, we learn &p by running the Expectation Minimization (EM) algorithm (Dempster et al. 1977) with samples from /re If J 1, LRDS reduces to Gaussian LRDS (G-LRDS), and we learn parameters (m), ) 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 summa- rize the whole GMM-LRDS methodology Algorithm 3 provided in Appendix A and illustrate in Figure 2 the behaviours of GMM-LRDS and G-LRDS for bi-modal target distribution. Our results show that choosing a Gaussian reference distribution leads to mode collapse, while defining re as in GMM-LDRS enables RDS recover the relative weights accurately.
fixed in ad-

---

Published as a conference paper at ICLR 2025
Target
Gaussian reference
Gaussian Mixture reference
Mode weight 70.0%
Mode weight 49.8%
Mode weight 50.0%
Reference samples
G-LRDS samples
GMM-LRDS samples
Mode weight 50.0%
Mode weight 0.0%
Mode weight 68.1%
Figure 2: Comparison between G-LRDS and GMM-LRDS. Here, the target distribution is the same 16-dimensional Gaussian mixture as in Figure (top left), see Appendix H. for more details. For illustration purpose. projections along the first two coordinates of "Mode weight' refers to the effective weight of the left mode. Reference samples (bottom left) obtained by rumning MALA sampler initialized in both target modes: each color (orange/purple) depicts one MALA Markov chain. In particular, none of them mixes hetween the modes. Running G-LRDS with 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).
used. In each cell, the value
3.3 ENERGY-BASHR 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 correct approximation of qrel, regardless of the number of reference modes J. see Appendix 1.4 for illustrations. To provide more flexibility, we propose a second version of LRDS where the reference distribution is parameterized by Energy-Based Model (EBM), i.e., red( (I) = exp(-E*(x) where E Rd + R is a neural network with parameters and the normalizing constant g J&x exp( EV (x))dd 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 ne by p insist on its parametric nature.
At first sight, it seems natural to learn the density nd and the scores (st 16 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 Score Matching (SM) objective with samples
from ard (Hyvärinen & Dayan 20055 Vincent 2011). On the other hand, one could learn p by maximizing the standard Maximum Likelihood (ML. objective LLL, : o + [logp" (X)), where
X ~ tto,, with corresponding gradient given by & +
EV, (X+)) +),, where
X- ~ ~p and X+ oto Nevertheless, computing gradients of objective L requires to sample from pf (negative sampling), which a well known hurdle EBM training since p is expected to be as multi-modal as se and thus π.
Therefore, we rather suggest to learn a path of parametric distributions (pî )ee(o.r] as the marginal
distributions of )e defined the noising process induced by SDE ①I. starting at X ~ frel, To do so, we set (pi )eefo,r, as a multi-level EBM, i.e., PÝ (x) = (t, where EY [o, 7] x R2 + TR is now a time-dependent neural network with parameters yp and the normalizing constants 2 f exp( F ())x are still intractable. In the RDS framework, exp(-E* (0,,)) x)) and sx) = -V (t, To learn this multi-level EBM, we propose to maximize the integrated MI. objective oP - 1ed Ellog pÝ (X In this case, we can leverage the correlations between the single-level EBMs to alleviate their in- dividual issue of negative sampling, a strategy that has already been investigated in several works
this amounts consider 7 (z)
(ao et al. 2021 Zhu et al. 2024 Zhang et al. 2023).

---

Published as a conference paper at ICLR 2025
Here, since (Xx")) defines path of increas- ingly simpler distributions, annealed MCMC sam- plers can be conveniently used to sample from the multi-level EBM densities. This allows us to kill two birds in one stone: (i) obtain nega- tive samples for each single-level EBM, which is needed to compute the gradient of the ML objec- tive, and (ii) we overcome the individual sampling issues at every level thanks annealing. This method is completely detailed in Algorithm [ of Appendix F together with previous literature on multi-level EBMs. Morever, Appendix E provides details on annealed MCMC samplers.
Target
GMM
EBM
Exact samples
GMM-LRDS samples
EBM-LRDS samples
This version of LRDS called Thergy-Based Model LRDS (EBM-LRDS) is summarized in Algo- rithm 4 of Appendix A In Figure we illustrate the superiority of EBM-LRDS over GMM-LRDS for target distribution that exhibits complex ge- ometry. While the GMM fails at capturing the lo- cal energy landscape of R. which results in poor performance of GMM-LRDS, EBM-LRDS cap- tures well the target distribution since the reference EBM recovers the complex geometry of the target.
Figure 3: Comparison between GMM-LRDS and EBM-LRDS in a multi-modal setting. Here, we target the 2-dimensional Rings distri- bution, which has 3 unbalanced modes repre- sented by the rings. (Leftt): Target density 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).
density (top)
4 RELATED WORKS
In the past years, numerous parametric methods have been proposed to approximate a diffusion process (EO)) ~ TP induced by an SDE with intractable drift, such ( bridging easy-to- sample distribution hhose to a target distribution æ. Based on deep learning techniques, they suggest to learn a path measure Ip", induced by a neural SDE that admits a neural network (parameterized by 0 E θ be optimized) as drift function. Below, we review two main classes of such algorithms.
Variational diffusion-based methods. To optimize 6, several approaches have adopted varia- tional formulation on path measure space, which consists in minimizing a divergence-based loss with samples from ps We adopt the same perspective in the present work. For instance, recent line of works has defined IP as the time-reversal of a noising diffusion process, such as PIS (Zhang & Chen 2022). DDS (Vargas et al. 2023a) and Time-reversed Diffission Sampler (DIS) (Berner al., 2023 In these methods, the neural network used in Fd 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 costly procedure, since each training step requires a full simula- tion of ipo 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 target-informed parameterization. On the other hand, Vargas et al. (2024) propose a different perspective by defining as a controlled version of an annealed Langevin diffusion, that is expected to follow a prescribed path of tractable marginal den- sities. Here, both IP and I are induced by neural SDEs. The authors present two versions of their algorithm. Controlled Monte Carlo Diffusiio (CMCD), using either reverse KL 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 Appendix D.2
erregenne,,Riccher
Adaptive diffusion-bassee approaches. To alleviate the computational difficulties of divergence- based losses, recent works have proposed to learn Pd with an adaptive procedure (Phillips et al.) 2024 Akhound-Sadegh et al.] 2024). These methods iterate two steps which consist of (a) ob- taining approximate target samples by sampling from pd and (b) optimizing / using those samples via learning techniques usually restricted to generative modeling. In particular, AAkhoundSSadehh et al. (2024) present Iterated Denoising Energy Matching (iDEM), where I corresponds to the time- reversal of Variance- Explodingg noising diffusion process (Song et al. 2021). In their setting, stage (a) is conducted by running the SDE induced by pd and stage (b) relies on a novel energy-matching loss which directly depends on æ. On the other hand, the Particle Denoising Diffusion Sampler
8

---

Published as a conference paper at ICLR 2025
(PDDS) (Phillips et al. 2024) (a) introduces SMC-based scheme when sampling from IP" 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 pe, 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 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 diffusiion-bsse methods. implemented with the LV loss - LV-PIS (Zhang
& Chen 2022). LV-DDS (Vargas al. 2023a), LV-DIS (Berner et al. 2023) and LV-CMCD (Vargas et a 2024); (c) adaptive diffusion- based approaches - -DEMM { Aakhound-Sadeee et al. (2224 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 mode Fl The details of each target distribution can be found in Appendix H.I
General experimental setting. To ensure fair comparison with previous approaches, we make sure that all competing methods are informed as LRDS of prior knowledge on the target distri- bution. More specifically: we set π hhaas as a Gaussian approximation of भ in SMC, RE or CMCD; we choose a based on Gaussian isotropic approximation of Ą for PIS, DDS or DIS (see Sec- tion 3.1: we standardize the target distributions using the empirical mean and variances of *nd for iDEM and PDDS. Additionally, we pre-fill iDEM's training buffer with samples from şre 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 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 π is expensive to compute. We refer to Appendix H.3 for complete details of the implementation of these methods.
base drift in CMCD;
High-dimensional Gaussian mixtures. We first consider a synthetic but challenging setting, where ı bi-modal Gaussian mixture whose modes are N(- E,), with weight w1 = 2/3, and N(1,, Σ2), with weight W2- = 1/33, where 1 is the d-dimensional vector with all components equal to 1, and for i E (1,2]. D; E _dxdd is a diagonal positive matrix with conditioning number equal to 100. For each method, the ground truth weight wi is estimated by Monte Carlo estimator 10. We report the estimation error (w1 = for increasing values of d in Table 2 Note that mode collapse occurs when Φú E {0, 1}. We observe that GMM-LRDS outperforms competing methods in all the considered dimensions. In Appendix I.3 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 ot from statistical mechanics. Next, we sample from the ID ot model, previously studied (Gabrié et al. 2022 Crenioux et al. 2024).. At the chosen temperature, the distribution has two well distinet 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-/w,, //,, and compare the results with Laplace approximation (0-th and 2-nd orders), see Appendix H. 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 Complete results showing the failure of the competing methods are provided in Appendix [.5.4] Given the satisfying results obtained with the lightweight GMM-LRDS sampler, we did not run EBM-LRDS for this target distribution.
In Appendix G we give expressions of additional variational metrics proposed by Blessing al. 2024 to quantify mode collapse the specific case of RDS.

---

Published as a conference paper at ICLR 2025
Table 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.
3.0
Laplace approx.
Laplace approx. (2nd order)
2.5
GMM-LRDS
1.5
Algorithm
d= 16
32
d 64
SMC
11.4%.9.19
15.8%
15.2%=7.3%%
1.0
RE
16.5%++.39
15.9% 14%
17.0%=11.4%
0.000
0.001
0.002
0.003
IN-PIS
6.0%
33.2% 1%
33.0%=20.12
LV-DDS
11.8%%49.3%
31.5%
13.1%.20.1%
LV-DIS
16.4%%t0.5%
16.5% 4%
16.8%=20.6%
Figure 4: Estimation of the relative
LV-CMCD
32.2%=15.49
50.1%%+..88
16.3%+10.699
weight of oi modes with increasing h, averaged over 16 sampling runs.
IDE
33.3%.10 0%
66. 7%.7%
11.7%620.49
PDDS
0.8%%
66.7%c+..0
N/A
GMM-LRDS
1.7%%%.6%
2.7%:2..8%
4.1 122.6%%
GMM LRDS
EBM LRDS
SMC
RE
LV-PIS
LV-DDS
LV-DIS
LV-CMCD IDEM
Target
O
O
O
O
O
O
Figure 5: Samples obtained for Rings distribution. Reasonable results could not be obtained with PDDS due to numerical issues.
LV-PIS
LV-DDS
iDEM
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, nents for GMM-LRDS and leverage the Replica Exchange algo- rithm 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 ex- act for the Checkerboard distribution, see Figure while being able to correctly estimate the ground truth mode weights.
GMM LRDS
EBM LRDS
consider J 64 compo-
Target
Figure 6: Samples obtained for Checkerboard distribu- tion. All sampling methods except LRDS fail to provide accurate samples.
Bayesian logistic regression models. Lastly, we provide in Appendix [I3.3] 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 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 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-worke sampling tasks on non-euclidean spaces which may benefit from the flexible definition of the reference process as EBM. An interesting test case of this kind that shall be considered is the sampling the Boltzmann distribution of proteins in interna coordinates.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGEMENTS
We would like to thank Lorenz Richter, Julius Berner, Denis Blessing, Junhua Chen, Francisco Var- gas 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 learn- ing: Foundations, generative models and non-conver optimisation when work on this paper was undertaken. Part the work of AD 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 Ageney. 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 (AD011))
REFERENCES
Tara Akhound-Sadegie Jarrid Rector-Brooks, Joey Bose, Sarthak Mittal, Pablo Lemos, Cheng-Har Liu, Marcin Sendera, Siamak Ravanbakhsh, Gauthier Gidel, Yoshua Bengio, Nikolay Malkin, and Alexander Tong. Iterated denoising energy matching for sampling from boltzmann densi-
ties. In Proceedings of the 41st International Conference Machine Learning, volume 235 of Proceedings of Machine Learning Research, pp. 760-786. PMLR, 21-27 Jul 2024. URL
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 Yc, Michael Volpp, and Gerhard Neumann. A unified per- spective natural gradient variational inference with gaussian mixture models. Transactions on Machine Learning Research, 2023. ISSN 2835-8856. URL
[orumpid=ttBjsx4E3
Julius Berner, Lorenz Richter, and Karen Ullrich. An optimal control perspective on diffusion- based generative modeling. Transactions on Machine Learning Research, 2023. URL https:
Denis Blessing, Xiaogang Jia, Johannes Esslinger, Francisco Vargas, and Gerhard Neumann. Be- yond ELBOs: A large-scale evaluation of variational methods for sampling. In Proceed- ings of the 4lst International Conference on Machine Learning. volume 235 of Proceedings of Machine Learming Research, 4205-4229 PMLR, 21-27 Jul 2024. URL https://
rroceddingg..ll.pprrss//// http://ss//////llssssnnn.......
Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Li Meng. Handbook of Markov Chain Monte Carlo. Chapman and Hal//CRC, May 2011. ISBN 9780429138508 doi: 10. 1201/b10905 URI.
ht.tp://ax.do..org/ wtt.:/....r..1...0/11o.o...
Patrick Cattiaux, Giovanni Conforti, Ivan Gentil, and Christian Léonard. Time reversal of diffusion processes under finite entropy condition. Annales de I' l'Institu Henri Poincuré, Probabilités et Statistiquess, 59(4): 1844 1881, 2023.
doi: 10. (0.1214/22AAHHP1332 URLhttps: /ttp///ooo////////////o/oooooo/.
10.1214/22 A111112010
Marco Cuturi.
Sinkhorn distances: Lightspeed computation of optimal transport. 26,
2013. URL
File/ F=D2A222077722221332222233333332
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.00 230.:0499.0 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 Arnad Doucet. Target score matching. arXiv preprint 2xXvv:242.008677 2024.

---

Published as a conference paper at ICLR 2025
Pierre Del Moral, Amaud Doucet, and Ajay Jasra. Sequential Monte Carlo samplers. Journal of the Royal Statistical Society Series B: Statistical Methodology, 68(3):41 4-66, 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 syn- thesis. In Advances in Neural Information Processing Systems, 2021. URL https://
Arnaud Doucet, Will Grathwohl, Alexander G Matthews, and Heiko Strathmann. Score-based diffu- sion meets annealed importance sampling. Advances in Neural Information Processing Systems, 15:21482-21494, 2022.
Simon Duane, A.D. Kennedy, Brian J. Pendleton, and Duncan Roweth. Hybrid monte carlo. Physics Letters B, 195(2):216-2 22,222,, 1987. ISSN 0370-2693. doi: https://doi.on 10. 00.1001600377--66331779999999999
URL http::::/www.sceeee. enwedicinnceeecco rect..om/ /ccc///ceenee//////ccccc
articcle/p1100370263333333
Alain Durmus and Éric Moulines. Quantitative bounds of convergence for geometrically ergodic
Markov chain in the Wasserstein distance with application 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- Adaptive Monte Carlo augmented with normalizing flows. Proceedings of the National Academy of Sciences, 119(10):e200942111, 2022.
Marylou Gabrié, Grant M. Rotskoff, and Eric Vanden-Eijnden. Adaptive Monte Carlo aug- mented with normalizing flows. Proceedings of the National Academy of Sciences, 119(10): e2109420119, 2022. doi: 10. 10. htt33/paa..21042000 https///ttttt//////oo/orrrrr//o/ wttp://p.ppp..g/ooo... abs/ 2os//10.173/// http:///aaa...1222.02
URL https:://////
Ruiqi Gao, Yang Song, Ben Poole, Ying Nian Wu, and Diederik P Kingma. Learning energy-based models by diffusion recovery likelihood. In International Conference Learning Representa-
tions, 2021.
Tomas Geffner and Justin Domke. Meme variational inference via uncorrected hamiltonian anneal- ing. Advances in Neural Information Processing Systems, 34:639-651 2021.
Tomas Geffner and Justin Domke. Langevin diffusion variational inference. In International Con- ference 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.
r.press/
Louis Grenioux, Maxence Noble, Marylou Gabrié, and Alain Oliviero Durmus. Stochastic lo- calization via iterative posterior sampling. In Proceedings of the 41st International Confer- ence Machine Learning, volume 235 of Proceedings of Machine Learning Research, pp. 16337-16376. PMLR, 21-27 Jul 2024. URL htRlhttps:/// ptttps://proceedis ngs.ml r.press/v235.. grenioux24a.htmi
Arthur Gretton, Karsten M Borgwardt, Malte 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.
12

---

Published as a conference paper at ICLR 2025
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-166) 2014. http:::/:m...ooo/pp//ppo/
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.
//openreview.mettee formm d= Tfw.T JAh
Xunpeng Huang, Hanze Dong, Yifan HAO, Yian Ma, and Tong Zhang. Reverse Diffusion Monte Carlo. In The Twelfth International Conference Learning Representations, 2024a. URL
Xunpeng Huang, Difan Zou, Hanze Dong, Yi-An Ma. and Tong Zhang. Faster sampling with- out isoperimetry via biffusion-based monte carlo. In Proceedings of Thirty Seventh Conference Learning Theory, volume 247 of Proceedings of Machine Learming Research, pp. 2438- 2493. PMI.R, 30 Jun-03 Jul 2024b. URL http:::://pproeeeengggeeee r.press/v2477 huang24a.html
Aapo Hyvärinen and Peter Dayan, Estimation of non-normalized statistical models by score match- ing. Journal of Machine Learning Research, 6(4), 2005.
Ghassen Jerfel, Serena Wang, Clara Wong-Fannjiang, Katherine A. Heller, Yian Ma, and Michael Jordan. Variational refinement for importance sampling using the forward Kullback-Leible di- vergence. In Proceedings of the -hirethSeventh Conference on Uncertainty in Artificial Intelli- http:://coooeeeingg....o
gence, pp. 1819-1822 1829. PMLR. December 2021.
v161/jerfel2la.hee.
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. Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 16955-16977. PML.999777.......M PMLR, 23-29 Jul 2023. URL
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 133-0-470-1777344
Christian Léonard. Some properties of path measures. Séminaire de Probabilités XLVI, pp. 207-230, 2014.
Jun 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. International Conference on Machine Learning, pp- 15299, PMLR, 2022.
Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Ed- ward Teller. Equation of State Calculations by Fast Computing Machines. The Journal of Chem- ical Physics, 1166):188711092 06 1953, ISSN 0021-9606. doi: 10.1063/1.1699114 URL
htpp://66...cg//00006/........
Laurence Illing Midgley, Vincent Stimper, Gregor N. C. Simm, Bernhard Schölkopf. and José Miguel Hernández-Lobato. Flow annealed importance sampling bootstrap. The Eleventh International Conference Learning Representations, 2023. URL https://
Radford M Neal. Annealed importance sampling, Statistics and computing, 11:125-139 1399, 2001.
Radford M Neal. Meme using hamiltonian dynamics. arXiv preprint arXiv: /20i /206. 1901, 2012.
13

---

Published as a conference paper at ICLR 2025
Alexander Quinn Nichol and Prafulla Dhariwal. Improved denoising diffusion probabilistic models. In International Conference on Machine Learning, pp. 3162-81711 PMLR, 2021.
Frank Noé, Simon Olsson, Jonas Köhler, and Hao Wu. Boltzmann generators: Sampling equilib- rium states of many-body systems with deep learning. Science, 365(6457)eaaw) 147, Septem- ber 2019. ISSN 0036-8075, 1095-9203.
doi: 10.1126/science.n aaweencee.awwllln... 1147
URL https:
Nikolas Nüsken and Lorenz Richter. Solving nigh-dimensionne hamilton- jacobi-bellman 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 Lak- shminarayanan. 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-leamn Machine learning in python. Journal of machine learning research, 12(0ct))2825-2833 2011.
Gabriel Peyré, Marco Cuturi, et al. Computational optimal transport: With applications to data science. Foundations and Trends@ Machine Learning, 1(5-6))355-6077, 2019.
Angus Phillips, Hai-Dang Dau, Michael John Hutchinson, Valentin De Bortoli, George Deli- giannidis. and Arnad Doucet. Particle denoising diffusion sampler. In Proceedings of the 41st International Conference on Machine Learning. volume 235 of Proceedings of Ma- chine Learning Research, pp. 40724. PMIR, 21-27 Jul 2024. URI. https://
proceedi ngs.m r.press/ v235/ 00223/pl0 phi hs24a.html
Emilia Pompe, Chris Holmes, and Krzysztof Eatuszyáski. A framework for adaptive MCMC targeting multimodal distributions. The Annals of Statistics, 48(5):2930- doi: 10.121419-A0S19166 URI.
2952, Octoher 2020. ISSN 0090-5364, 2168-8966.
48/
imodall
distrrlbuttions/00112224.. htt.://2111111111991.1....... Statistics.
20s1916.ful1 Publisher: Institute of Mathematical
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: wtt.:/wwwww.eeee/oooiieee--eeee-
D J Rezende, S Mohamed, and D Wierstra. Stochastic backpropagation and ap- proximate inference in deep generative models. Proceedings of The 31st Inter- national Conference in Machine Learning, Beijing China, 32:1278-, 2014. URL http: htl:://jn http://.iiiooooaa////////33377-- ASA-
http/////reeeee11/e 4.htm $5Cnpapers3:
//publ wwloo//aaioo///iooooooooo
arXiv:
1401.4082v3 ISBN: 9781634393973.
Lorenz Richter, Julius Berner, and Guan-Homg Liu. Improved sampling via learned diffusions. In
ICML Workshop on New Frontiers Learning, Control, and Dynamical Systems, 2023. URL
wewevveewwrrr ew. net. forum? beulgyDDiee0
Gareth O. Roberts and Richard I.. Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 2(4):341-363, 1996. ISSN 13507265. URL http:
/ www.:/stoo..rrggttaaeee/3331111.
Tim Salimans and Jonathan Ho. Should EBMs model the energy or the score?, 2021. URL https:
/
RNb
14

---

Published as a conference paper at ICLR 2025
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In International conference on machine learn- ing. pp. 2256-2265, PMLR, 2015.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distri- bution. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. URL https: Eile/ $11-990222777701377713acc5520122
Yang Song and Stefano Ermon. techniquu
techniques for training score-based generative models. In Advances in Neural Information Processsing Systems, volume 33, pp. 12438 12448. Curran Asso- ciates, Inc., 2020. URL http:pss//pcccee.....s........./
waaaapperrrieeeeppeeeeeeeeeeeeee files/paper/ 631 |00.0000000111111111100000000330
2020/file/
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://
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. Phys- ical Review Letters, 57(21):6077-2004 November 1986. doi: 10. 10.1100/hhyRRee.Lett...000
URL http:///hhaaarrr////00//////ssss American Physical Society.
Publisher:
Saifuddin Syed, Alexandre Bouchard-Côté, George Deligiannidis, and Arnaud Doucet. Non-
reversible parallel tempering: a scalable highly parallel mcme 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 temper- ing. 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 Conferencee on Learning Theory, pp. 3084-3114 14. PMLR, 2019.
Francisco Vargas, Will Sussman Grathwohl, and Arnaud Doucet. Denoising diffusion samplers.
In The Eleventh International Conference Learning Representations, 2023a. URL https:
7 htreenrveiieennett/roumi........ TAbuu.
Francisco Vargas, Andrius Ovsianas, David Fermandes, Mark Girolami, Neil D Lawrence, and Niko- las Nüsken. Bayesian learning via neural Schröddinge-FFöllmm 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 Learn- ing Representations, 2024. URL U http:::::///penneieeee..eee
Pascal Vincent. A connection between score matching and denoising autoencoders. Neural compu- tation, 23(7):1661-1674, 2011.
Martin J. Wainwright and Michael Jordan. Graphical Models, Exponential Families, and Variational Inference. Foundations and Trends@ in Machine Learning, 1:1-305, 2008. doi: 10. 10111/2200000000..
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 Informa- tion Processing Systems, 0999333-994,, 2020.
Qinsheng Zhang and Yongxin Chen. Path Integral Sampler: a stochastic control approach for sam- pling. In The Tenth International Conference on Learning Representations, 2022.
15

---

Published as a conference paper at ICLR 2025
Xinwei Zhang, Zhiqiang Tan, and Zhijian Ou. Persistently trained, diffusion-assisted energy- based models. Stat, 12(1)==625. 2023. doi: https//ddbiorg//..100//ttttuuu URL https:
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: http:////pppppeeeeeeeeeeeo/e/oaa
16

---

Published as a conference paper at ICLR 2025
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 computa- tions. Appendix C we describe the noising diffusion schemes considered in this work, namely the Pinned Brownian Motion and the Variance-Presse Preserving setting, and propose novel schemes. For RDS and related variational diffusion-baset methods, we dispense details on their continuous-time and corresponding discrete-time variational objectives in Appendix [ We respectively provide an overview of existing annealed MCMC methods and multi-level Energy-Based models in Appendix F 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 implementa- tion of RDS and other competing methods given in Appendix
Finally, present additional
theoretical and experimental results in Appendix
Our codebase is available at hhttps://giit httt::s////tttb..coomm h2w.hhoo442ee0 sampler Irds
Notation. Let E (CT) and E [0, 2J. For any x, E IRd we denote by pa E (CC)) the path measure associated to the stochastic process (Xt))E(o.T) P conditioned X, Ir. Furthermore, for any E (RR),
Pir CT stands for the path measure Jg= Pi d))()) For any K 1, denote the of joint distributions P ((R k+1) by g(K+1), Finally, we adopt the following notation to design sample batches: for any L 2 L 2 0 and K > K 2 0, we denote by L K:R the
sequence samples {X R.L -.KE
A PSEUDO-CODES OF RDS-BASED ALGORITHMS
We respectively give sampling procedures and training procedures of general version of RDS in Algorithm and Algorithm 2 Relying on this, we derive the complete training schemes for GMM- LRDS (Algorithm and EBM-LRDS (Algorithm
Algorithm 1: Refereecce-basd Diffusion Sampling (RDS): sampling stage
Input: Time discretization ] of [0, T). reference scores - neural network
K-1 9º [0, x IR Rd, variational coefficients (ax, by, C )= defined Appendix D.2
b Initialization
Yu T K- (=x)===
(Zx
N(O, Ta) O ..... K
==
for
K- 1 do
> Compute the l--th diffusion step
Ya+1 Ya hn ( (Vx) oZ- (Y+) ))))+ + VoET
Output: Discrete time process Ya.k approximating (Y;
Algorithm 2: sefference-based Diffusion Sampling (RDS) training stage
Input: Time discretization (4x j of J0, T]. target density 7, reference density 7 and scores
( 1kzo neural network d [o, T] Rd - Rd initialized such that 9 = 0, number training iterations N, batch size B, variational coefficients (wk, ag. by, Cx k= defined Appendix [D22
for 0....,N do
> Simulate B trajectories of the process (Y7") )eejur 1:B i.d. SamplingRDS(g"1 see Algorithm 1
b Apply stochastic gradient descent step on 0, Compute MC estimator (RDS(O.)
loss LRDS (0..) defined (7) with samples L:
Compute the gradient V & (f.)) and update 9., u with Adam optimizer
Output: Learned guidance 9
17

---

Published as a conference paper at ICLR 2025
Algorithm 3: Gaussian Mixture Model LRDS (GMM-LRDSS) training stage
Input:
Time discretization {tx K E== of [0, TI, target density 7, location of the target modes )= number of Markov chains per target mode M 2 1, size of Markow chains Noot 2 1, effective size of Markov chains Nef E [1, No numher of reference modes J
P. (a) Obtain reference samples
For each i E {1,.... 7), build via MALA M Markov chains (X O:Niot |m= of size Niot starting Keep the last New samples of each Markov chain define and = (X & (b) Define the reference process
I.M
Noot -Nar--N 1m1mm1
if J = then
Fit on nT Gaussian model with parameterized density 7th (Maximum Likelihood estimation)
else
Fit on nT J-component Gaussian Mixture Model with parameterized density 7 (EM algorithm)
Set get = 7 and compute { It k analytically, see Appendix
K-1
P (c) Run the RDS variational optimization
Run 2
Output: Trained GMM-L.RDS sampler
Algorithm 4: Energy-Based Model LRDS (EBM-LRDS) training stage
Input:
Time discretization {lx ) of [o, T). target density n. location of the target modes (xx 1/=
number of Markov chains per target mode M 2 1, size of Markov chains Noo: 1, effective size Markov chains Neæ E [1, Nioe]
> (a) Obtain reference samples
For each E [4). build via MALA Markov chains [X }M
size Niot starting at x, J.M
Keep the last Net samples each Markov chain to define aTd 0e (X Mor--Natt::oo I=l,m=1
(b) Define the reference process
Based on nn fit multi-level EBM E using Algorithm
Set (x) = exp(- E0, x)) and s7-1,, (z) = V (T 2)) for any E {0...., K 1}
b (c) Run the RDS variational optimization
Run Algorithm
Output: Trained EBM-RDS sampler
B PRELIMINARIES
Linear SDE integration. We first dispense useful lemma to compute exact integration in SDEs with linear drift.
AO (Integrabilityy conditions on f and (). Coefficients f: [0, + R and 3 [0,2] + (0, oo) are such that (a) f is integrable on (0, T) and (b) B is integrable on (0, T).
Lemma 2. Let T > 0 and b E Rdd Consider the SDE defined on T] by dy, = f(t + b)dt + / BB where coefficients f and B verify AO Then, for any pair of time-steps (s. t) such that T 2 >>> 2 0, the conditional distribution of Ye given Y, Ya E Rd, denoted by pus (llw),
verifies
(1))) = (exp(S, f(u)du)y. + (exp(I] f(u)du) - 1) b, 1tt B(u) exp(2 J f(r)dr)dula).
Proof. Assume 0 Define the function Ç : t E [0,T] + exp(- Jd f())du) and consider the stochastic process (Z4)eeoo.r defined by Z, = (()),, for any L T). By ito's formula, we have dz2 = ((t)((t)bdt + <(t) (())dB = ((t)bbt + (Bt))BB Therefore, for any time-steps
(s. 1) such that 2 > 8 2 0, we have
<(t)); <(s) (C(s) - [(t)){ + f" (
and then
Y, exp(S. f())du)Y, + (exp(S.' f(u)du) - 11)++ 6+ 1 BB(u)exp(J. ((r)drd
which gives the result using îto's Ito's isometry and the fact that Y, is independent from (B, - B,) B.)re(s.T)
18

---

Published as a conference paper at ICLR 2025
General SDE integration. Consider general SDE of the form dy, (X) + dx + /B(t)dB, defined on [0, 7), where f. / satisfy A and g [0, x TRd + TR is black-box function. Let (s,t) be a pair of time-steps such that T 2 <> conditional distribution of Y given Y, = Ma Rd. 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
O. We aim to compute the
(a) Euler-Maruyama (EM) scheme, which consists of exactly integrating on time interval [s, 4] the SDE defined by dY,, =
BB(s))BBB with constant drift and constant volatility.
(b) Exponential Integration (EZ) 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, = [()) + (o))))) BB(jdB,, with linear drift, where the intractable term is frozen at ys (Va).
Computation of scores for Gaussian mixtures. Below, we provide a useful lemma to compute the marginal scores of noising diffusion process applied to general Gaussian mixture.
Lemma 3. Let π E P (Rd)) Consider the noising process (Xe))e(o.T] defined by SDE (2) where Xu ~ a. Then, for any t (), the conditional distribution of X given Xo Io E Rd is defined by the Gaussian kernel
ago (2)) =
where S(t) exp(So f (u)du) and o² (t) = So ((u)/S(u)ddu.
In particular, a is Gaussian distribution i.e., = N(m, E) with mean E R and covariance ) E Rdxd, then any E [0, T1, the marginal distribution of X2 is Gaussian, defined by Pt = N(S(t)m,
(tt)
(10)
On the other hand, if π is a Gaussian mixture ie.. π = E₁ E,=,,,((m Nym, Ey) where 2 is the number of components, {wz}/=1 E As are the component weights, m, E Ed and D3 E gaxd are respectively the mean and the covariance of the j-th component, for any [0, the marginal distribution of X, is a Gaussian mixture, defined by
Pt
wy ((5((()))))(())5)) + (())))))))))))))
(11)
Note that differentiating provides closed-form solutions of the score of the stochastic process.
Proof. The proof of 9 directly follows from Lemma Then, 10) and [1 are obtained by apply- ing the affine mapping defined by the transition kernel 9
to the target distribution n.
In particular, II demonstrates that any noising diffusion scheme applied to Gaussian mixture does not show any mode switching.
Integral probability Metrics. Finally, dispense the definition of standard probability metries, that will be used to assess the performance of the considered samplers, when having access to ground truth samples. Consider two distributions y 30 (RR) and 9 (IR") that aim to compare.
We recall that the 2- Wasserstein distance between and v/ is given by
W2(1,, inf ([ || xn
20l122 dm (No, x1) (Rd x Re), no 14, m u31/2,
where m, denotes the i-th marginal of n for E {0, 1}. In practice, we rather turn to its entropy regularized version (Peyré et al. 2019)
Wa.- (1, v/) inf [e< || =o ( - H()) : (Rd Red). 14, = v/
where 8 > is a regularization haper-parameter and HT) Sx* x log T x1 ) (z0, (1)
S
refers to as the entropy of When having access to samples from f and v, compute statistical
19

---

Published as a conference paper at ICLR 2025
estimation of (2,,(1,)) 2) via Sinkhorn algorithm (Cuturi 2013) with 10-3, based on the GPU- oriented implementation from https:: / http:///i/n///aaaaaaaaaaaeeo////
sinkhorn
We also consider the Maximum Mean Discrepancy (MMD) (Gretton et al. 2012). which quantifies the dissimilarity between j and by comparing their mean embeddings in reproducing kernel Hilbert space H4, with kernel k, as followed
MMD(H, 2)
sup JEHx-Il/flx $1
([f(X) - f(O)) X~ M-MYYNV
In particular, (M,D(,,)) = 0 if and only if μ almost surely. When having access to samples (xr )== from ju and {yr] from v, we compute an unbiased statistical estimation of MMD(μ, v) given by
2 m]n-1] E Eisicign 2,) + m[m-1] E Eisicssm k(y,
mm j=1 k( By)
where k(x,8) = exp || yll (2a)) is the squared exponential kernel whose positive bandwidth is fixed as the median the squared pairwise distances computed on the joint set {2,) U {ys 3 (Gretton et al. 2012).
Finally, we recall that the Kolmogorov-Smiroov (KS) distance between μ and v is defined by KS u) sup (x) Fo(x))
ZERd
where Fn (respectively F,)) denotes the cumulative distribution function of j (respectively v/). When having access to samples from μ and v/, we implement statistical estimation of the sliced KS distance (Grenioux et al., 2023) using 128 random projections.
C DETAILS ON NOISING DIFFUSION PROCESSES
In this section, consider noising diffusion process (X,)) induced by SDE for an arbitrary target distribution (RR)) and T > and denote by P E (CC) the corresponding path measure. Below, we provide detailed results on this diffusion process. depending on the choice of f and 3: we first provide details the general setting in Appendix C.I then, we fully describe the Pinned Brownian Motion and the Mariance-Preserviie 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.I GENERAL NOISING SCHEME
Here, we consider the most general form of SDE (
Lemma 4. Let E 9 (Rd). Assume that and 3 both verify AD Then, for any pair of time-steps
(s. t) such that T 2 t> s 2 0, the conditional distribution of X, given X, = Fg E RRd. is defined
by the Gaussian kernel
(1..) N (exp(log S(t) log S()) ((s))
(12)
where S(t) exp(J] f (u)du) and o² (t) 1n B(u)/S(w)2du.
Proof. The proof of (12) directly follows from Lemma 2
C.2 PINNED BROWNIAN MOTION
For any t [0, T), define att) S a())dt and consider the where f(t) a(T)-o(0) Then, SDE can be rewritten
3(1)
8())x, dt VB(J)W, X m. dx, aT)) a(t)
20

---

Published as a conference paper at ICLR 2025
This noising scheme, known 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 b-transform on the uncontrolled SDE dX, = V (t,, associated to path measure Q. to O at time Interestingly, can be shown that the resulting path measure (P*) is solution to the following stochastic optimal control problem
arg min Q) PP P(=
which is often referred to as a half Schrödinger Bridge problem.
Lemma 5. Let E P (R.). Assume that ft) a(i) and B verifies AD Then, for any pair
B(t)
of time-steps (s, 1) such that T 2 > 8 2 0, the conditional distribution of X given X, = Is E TRd defined by the Gaussian kernel
(1((-(r))) - N ((T) - a(t)
(T)) (s(s) Tsa (a(T) ((())(a(()) -0(s)] -a))
(13)
a(T) 1 a(s))
Since px)( Su (010 (xo), it results that IP.. do-
Proof. The proof of 13) directly follows from Lemma
Following Lemma 5 the PBM is an 'exact' noising scheme and we have huse = 5o in this setting. Moreover, under mild assumptions
the time-reversed SDE writes
dY, = B( ( aT - Vlogpi- di VB(T / Yo = 0.
Y (T))
On the choice of the 3-schedule. Previous works have considered constant schedule 3(() o² (and therefore, axtt) = o²t), where a > 0 can be arbritrarily chosen, see e.g.. Zhang Chen 2022); Richter et al. (2023). We also follow this setting in practice.
C.3 VARIANCE- VARIANCE-PREESERVING PRESERVING DIFFUSION
For any & E TJ, denote (t) = 1-d (())dt and consider the case where f(t) = -B())220² with a and being such that fo (ss)ds >> 1. Then, SDE ( can rewritten as
(t)) dt + B(t)dW, Xo x.
dX, 20²
This noising scheme, known as the Variance- Variance-Preserrvnng (VP) scheme (Song et al. 2021) and previ- ously considered by Vargas al. (2023a), is -eattin-Uhhennbck diffusion process and is largely used score-based generative models.
Lemma 6. Let (RR)) and O. Assume that f(t) 22t and verifies such that Jo B(s)ds >> 1. Then, for any pair of time-steps (s, 1) such that T 2 >>> s 2 0, the conditional
distribution of X, given X, = Rdd is defined by the Gaussian kernel
(-1x(-1)))) (VI->e, Ao, 02 >e,eedd)
(14)
where Autt = 1 exp(a(s) a())) Since (x) Sga (r10(c/2o)d it results that Pj N(O, 12)-
Rd
Proof. The proof of (4) directly follows from Lemma along with fact that Ao RV 1.
Following Lemma 6 the VP scheme is an 'ergodic' noising scheme, converging exponentially fast to the Gaussian distribution N(0, o² La); therefore, we have base N(0, o² La) in this setting. Moreover, under mild assumptions on n, the time-reversed SDE writes as
3(T t) (Y, + 2 log (X))d 3 Yo F
dy 2
after simple linear time reparameterization.
21

---

Published as a conference paper at ICLR 2025
On the choice of the 3-schedule. Previous works have considered linear schedule ((t) Smin t/T) + Smax (t (T) where Smin = 0.1, 8max E {10, 20} and T = 1, see e.g.. Song et al. (2021) Richter et al. (2023); Reu et al. (2024) cosine parameterization (Nichol & Dhan- wal 2021] Vargas et al., 2023a), which has been proved to perform hetter in generative modeling. did not observe any significant difference between these set- tings, Hence, we fix the linear 3-schedule to be the default setting for our numerics (except DDS, originally implemented with the cosine schedule), and let a be arbitrarily chosen.
In sampling experiments,
C.4 ADDITIONAL EXAMPLES OF NOISING DIFFUSION PROCESSES
Below, we present original noising schemes, but did not consider them in our experiments.
Pinned Orustein-Uhlenbect Motion. We first propose to consider the case where f(t) =
2002 c A(u) du). Then, SDE
( can be rewritten as
B(1) dX,, 2 coth(S,' 2 du) + /B(t)dW,, Xo NT.
p(w)
Analogously to the PBM, one can show that this noising scheme can be obtained by applying Doob's h-transform the Armstein-Uhlenber process induced by the SDE dX, - X,dt + /3(()W, to hit 0 at time T. We refer to this noising diffusion scheme as the Pinned Ornstein- Uhlenbeck Motion (POUM).
Lemma 7. Let n E (RR). Assume that f(1) = - a(t) coth( du) and verifies A Then, for any pair of sime-steps (s,t) such that T 2 >>> 8 2 0, the conditional distribution of X given X, Ix E Bdd is defined by the Gaussian kernel (-10)) given by
N (( siinhss"" AQm du)]
Iz. Te, sinh( f-" ((w) du)² coot(( f-"
2 d() coth( 2f dud du) Le
sinhhs." zime du)
(15)
Since (z) fge (r) it results that IPG. = Sø-
Proof. The proof of (15) directly follows from Lemma]
Following Lemma the POUM scheme is an 'exact' noising scheme and we have hase Sy in this
setting. Moreover, under mild assumptions on æ, the time-reversed SDE (2) writes as
dY: 8T -)) (coth( 2
8(m) 2 Ye + 2V log (X))dt + / /(T - t)dB,, Yo = 0.
(0, 20) be a square integrable function on (0, Denote atl) S {())du. Inspired by Dai et al. (2023), we propose to consider [(t) =
Gaussianized Pinned Brownian Motion. Let [0, T
((T) and (()) a(T)-a(() Then, SDE can be rewritten as
<(A)
20(0)
((1)X8 dx, a(T) a(t) di + (T))) - a(t)
2C(1)
dW,, Xo TO
Lemma 8. Let π E (Rd)) Assume that f(t) -(()- and 3(t) ((T) where
c(t)
20(0)
Ç: [0, (0, oc) is square integrable function (0, T). Then, for any pair of time-steps (s, t)
such that T 2 t> 2 0, the conditional distribution of Xx given X, = Is E lad is defined by the Gaussian kernel
(1)(-1xa)) = N (a(())a(())
a(t)) ( (TT)) (x)) Mas (1 (1-10(2)-a - -(s)))
la
(16)
Since (z( = far (r)o((1/oa)dd((( it results that IP = N(O,La))
Proof. The proof of (6)) directly follows from Lemma
22

---

Published as a conference paper at ICLR 2025
Following Lemma 7 the noising scheme given above is 'exact' with base distribution given by hase (o, Ia). this setting being quite unusual in diffusion model community. Based on this, refer to this noising diffusion process the Gaussianized Pinned Brownian Motion (GPBM). Moreover, under mild assumptions on π, the time-reverse SDE 2 writes as
-t -t) dy, a(T) aa 1)
2C(T - 1) V (TT)) -) dBr Yo dBd
{Y, + 2V (X)) dt
N(O,
D VARIATIONAL OBJECTIVES FOR DIFFUSION-BASED METHODS
In this section, first provide details on the continuous-time framework of variational diffusion- based methods (including RDS) in Appendix D.I Then, we describe in Appendix [.2] the RDS discrete-time formulation, related the variational loss LRDS given (7). and compare this frame- work to discrete-time formulations of previous variational methods.
In the rest of the section, we consider target distribution т (RR) and recall that P denotes the path measure associated to the noising SDE II initialized at Xo ~ π.
D.1 CONTINUOUS TIME FORMULATION
D.1.1 RDS SETTING
Let th E (R.d). We recall that jprel denotes the path measure associated to the noising SDE starting at Xo xeef First note that IP* @ pref Therefore have the following relation F di .pef i.e., Pa is absolutely continuous with respect pref, with Radon-Nikodym derivative given dr
by (dP*/dr )dpe)))))))pppp (X (dz/dx (Xo). for any diffusion process (X,))ej0.7] with forward time direction.
Below, we first provide the proof of Proposition
Proof of Proposition ! Assume that 1.7. = 1Pgel = base and there exists 0* such that gf = gt- Under those two assumptions, it is clear that (Y," )<EO..T] i.e. O achieves optimal solution. We now detail how to obtain (6)
~ (F*) therefore, we have LLv (f*) = 0,
Consider the detached diffusion process (Y:²) (E(0,1 ~ pd. Using the fact that P log-ratio in the objective (5) can be computed as
dw dip
peet the
dPd log d(Ip)) R (.2z)) log
(Yd
dIpo
(M10.2) + loo 2Td (x-9) -log grel
log Z
(17)
(Frer) nt
Moreover, by applying Girsanov's theorem the path measures (Fret) and po with the assumption prof = n basc we have
d(prrt) (0.m) (
dpp
T B(T
(1) 2
at {" VB(T (V)) dB, (18)
log
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 ( by combining (17) and [8)
In contrast, previous reference-based works from Zhang & Chen 2022 Vargas al. (2023a) relied on the reverse KL continuous time objective
dip" d(Ipt)
C (0) (IP (1*)")= = E log
(
Yoo ~ pd
(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 presented next.
23

---

Published as a conference paper at ICLR 2025
Proposition 9. Assume there exists f θ such that of. 9t. Then, the loss defined in
19
achieves optimal solution at 8* and, setting & = log(? /n), it simplifies as
Cx) = 5 B - 2
(x,)l dt + e(x?))
(20)
where the equality holds, up to constants independent of 0.
Proof. For any 0 E O we have by the KL chain rule
k (0) KL(TT hase IPT) Ey
KL(P%% (PP))))
where the first term does not depend on 0. Then, it is clear that CkL is minimized when 0 0". The result of (20) directly follows from and [8): in particular, the additional constants are log g and log 2.
one may recognize a standard formulation of stochastic control- affine problem with terminal cost where the uncontrolled process is given by pret.
Based on the simplification 20
D.1.2 DIS SETTING
In their paper, Berner et al. (203)) seek to approximate the time-reversal of IP* for the VP noising
scheme, see Appendix C.3 To do so, they consider the class of variational path measures (P" Jece C
g (CC) where for any B E θ, p is associated to the SDE
dy - 1) yo + 3 -Y)
BBB(T l)dB, Yo N(0, (21)
where the 8-schedule is detailed in Appendix C.3 and s aims approximating the marginal scores P i.e., s (t.,) Vlogp; The resulting variational
[0. x Rd is a neural network that
method, DIS, then consists in minimizing the KL-based objective 0 + KL(P" (P+)*), later
extended by Richter al. (2023) to the LV-based setting. We refer to (Richter et al.
2023 Section
3.2.) for detailed expressions of these continuous-timg losses.
D.1.3 CMCD SETTING
Alternatively, Vargas al. 2024 propose to approximate diffusion process with marginals pre- scribed by a curve of distributions (TT). such that To = π base and Ty = æ, with unnnormalized densities that assumed to be tractable. To do so, they consider the class of pairs of variational path measures ((P9, (2) ()))
ece C (CT)) where, for any E O, Qo is associated to the SDE
dx
(V log - - (x adW,,
X T
(22)
2
and pd is associated the SDE
x²
dy = (V
+ (Y,")) adB, Yo ~ hase
(23)
2
with h" [ x Rd + Rd being a neural network. The resulting variational method, CMCD, consists in minimizing the KL-based objective 0 - KL(F" | () 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 (4x 1k such that In 1x T and K > In the rest of
the section, we denote &g = f 1x-
In the only of PBM, consider to to avoid numerical integration error (typically 10-1).
24

---

Published as a conference paper at ICLR 2025
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.I). EI integration in the PBM case (Appendix (.3.22 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.4land CMCD in Appendix DD.5 All of the presented settings implemented in our code.
D.3 GENERAL RDS SETTING
Let E @ and let nt E (Rd) be an arbitrary reference distribution. Consider the variational diffusion process (Y2 )ejo. defined by SDE (4. We recall that under mild assumptions on n". see e.g., .es..Cattiaux al. (2023), the time-reversal of pee is associated to the SDE
dy f(T t)y,' dt B(T t)s (Ytet) dt
ttdB, yo Pr
(24)
Assume that we given transition kernels (A, yk) PK+ (Aly)e) and (A, yk) + Pk+1 (Alyk) respectively approximating P(Ya., E AA Y = Wx) and P(Ytet (241 E A Y = uk). Then, for K sufficiently large, the path measures 1oo and (Loof) Rt may respectively be approximated by the joint distributions Pg:K- E (K(K+)) and PUKk E ((K+1), defined by P6:K = π hake [ Pk+IjK and ref = π. hase II ref Pk+z Using the relation P = dr Ireef (P*) may also be approximated by the joint distribution Pū.k da Po:k E (KK+))
pff
Using this correspondence between path measures and joint distributions, we propose to approxi- mate LIy defined in (5) by the discrete time LV-hased objective
dy dp
CRDS (0) Var log
(Y6.k)
Yo
Pu.k
(25)
In similar manner, Cku defined in 19 may be approximated by the discrete time KL-based objec- tive
dpb-k
Xx = E
log dpú.k
(Yolk)
Yo.k ~ P...
(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 dpo.k (Yd.k) = log dp (Yo.k) + log dp
2nd (Yx + +log Z - log z
(27)
dp::K
where 0 E (0,0). In the following sections, building upon (27, we derive the exact expression of the variational loss (5 for each combination of noising and discretization schemes. We emphasize that these derivations are particular applicable to PIS and DDS settings by taking re as given in Table
D.3.1 EM-BASED RDS
In this section, we consider an arbitrary noising scheme that fits within the framework detailed in Appendix C.I For any l (0, ..... K 1}, define the coefficients
wk B(T tx)
FM
FM
S(T Ix J8k
JMM B(T 1x )ok
ok
3(T
=x
To approximate SDEs (4 and 4 time interval [ta, t.+ we may define the following transition kernels
Pk+1jk (14k) N(a yu + BEM ster s (yk). IMI2)
(28)
Ph+1ije (2w) ) = N(a + 6 ( (yx) + gf - (yx)}, ck La)
0
(29)
25

---

Published as a conference paper at ICLR 2025
that be obtained by applying the EM integration scheme, Appendix B In this case, the log- rog-ratio log dp" /dpt can be computed as follows.
Lemma 10. Assume that the joint distributions pe and pdd are respectively induced by the transition
kernels defined in 28 and 299. Then, for any yo.k E (RR))++, we have
dp log dp (vo:k)
K-1
EM
K-1
Σ
k
19 -t
(yx))1
+
Σ
V /f gi (vb)
k=0
k=0
where EK
yk+1
Ok
FM
yk
by 1t (yk) + - (vx13)
(30)
Proof. In this proof, for ease of reading, we will respectively denote dx b% EMM and u EM by ak. bks Gj. and Wk. Let yo.k E (Id))++ We have
log dp" dpref (Vo.k)
K-
dp + | (3k- 18%k)
=
log
dy (2/=-1 l16)
||xx11
1
lle + (bk/ (ol//e))rruu (1)))
2
K
[1
9T (11x) ck: Tek.
20k
k=
where eg is defined by 30
The result of Lemma 10 directly follows by noting that Wis = /2/Ck.
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
K Σ
(RDs( = Var
EM w
-1 91-4 (Yx)
107-1 (Va) (
- -u (Y+)
FM w -r- (Yk) Zx log (Vk)
k=0
where {7x K-1 are independently distributed according to N(O, Ia), and (Yu) is recursively -]}
defined by Yo ~ hlase, and for any E (0..... ...... K
Y
FM FMY2 dg Ya
bj
ref (Yx) + of (Yw) + VOBM Zz ZM (1T
(31)
with being a detached version 0
k=0 Proof. Consider the discrete time process {Ya 1 defined by recursion 31
With the same nota-
tion as in Lemma 10 we have for any l E {O....,
Wo ( (Yk) -
Z +
recalling that w EM (bEM)2 /C Using 27
finally obtain the expression of the RDS loss.
For completeness, we also provide a similar result for the KL-based objective given in 26 Corollary 12 (Expression of reverse KI. loss EM setting). When using the EM scheme to integrate the SDE 4 up to additional constants independent of 0, the discrete time KI. loss simplifies
K
EM k
ng (X 7
LKI (0) = E
Σ
(Y2)11 + log g
2
26

---

Published as a conference paper at ICLR 2025
where {Zk] K-1 are independently distributed according N(0, Ia), and (YP) Tk=0 is recursively k== defined by Yo Thuue, for any k E {0, ................................ K 1-1} 1}
FM sef (YE) + u1 (YE) + +V V& Zx R+1 = Ukk Y
FM
1
As observed by Vargas et al. 2023a), there may be significant bias raising from the use of the Euler-Manue Mareymaa discretization when using diffusion-bast baased sampling methods. This motivates us to consider El 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 define the coefficients
We recall the notation a(t) = So {())du. For any & E {o...., K 1),
(a(T) - a(T (a)) tx) a(T - tx+1))
PBM
w
a(T) a(T +x+1)
a(T) aTT [k11)
PBM
0g
a(T) aT - )
PBM aT - (x) - a(T (x+1),
PBM (a(T) - a(T - (4+1)))((TT - (b) a(T (2+1))
a(T) - a(T 1k)
To approximate SDEs kernels
24 and 44 on time interval (4x+ 12+1). we may define the following transition
ref Pk 1/1 ( (kk = N(a + -. -T- (yk) 1, PRMI.).
(32)
ph 112 (11Bk) = N(a yk + 1PBM ( (yIc) + -. (W:)). Ck
PBMT2)
(33)
that can be obtained by applying the EI scheme, see Appendix B In this case, the log ratio log dpp /dpetl can be computed as follows.
Lemma 13. Assume that the joint distributions pe and pdd are respectively induced by the transition kernels defined in and (3). Then for any yo.k E K have
log dp" dp (Jo:k)
= Σ
K-1
w
PBM
119f
K
Σ
v wffuu gf (uk)) TC
2
k=0
K=0
where 6 VeRRM
yk+t1
yk
1
WPa ( (8/ (1E) + - (uk))))
Proof. The proof is exactly the same as Lemma 10 noting that we still have w.BM (6PUM)3
/ PBM
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 Corollary II
Corollary 14 (Expression RDS loss El setting PBM scheme). When using the El scheme to integrate SDE implemented with PBM, the RDS loss simplifies
1 (Yr)
CRDS(B) = Var
Σ k=
PBM w
---
gt to (Yr) ( 20 i (Yr)
K
PBM
-1- (Yx.) Zx + log 2 (Yk)
W
?
k=0
27

---

Published as a conference paper at ICLR 2025
T== R-1 are independently distributed according Ia), and {Yx) is recursively defined by Yo ~ husee, and for any k E .......K--
where {Z
PBM Yk + bg LPBM (SRM(ST- 11T-1e (Ys) + 9 (Yx) + V& Zk Y++1 = dg
with being a detached version 0. 0.
For completeness, we also provide similar result for the KL-based objective given in (26)
Corollary 15 (Expression reverse KI. loss - El setting PBM noising). When using the El scheme to integrate SDE 4 implemented with PBM, up to additional constants independent of 0 the discrete time KL loss simplifies as
(0) =
K
W
PBM
9 (Y2))1 + log (X2) 7
7
2
K=
where {Zx k== are independently distributed according to N(0, Ia), and ( 1 Tk=o is recursively defined by Yo ~ hase and for any k E {o...... ............ -13
ye
PBMY0 0kk
DPRM (5 18 (Y,) - (Y2)) + v o
CBB Zk
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 a<t) S ((u)du. For any k E fo...., K 1},
define A = exp ((al)
4)) a(T - (4+1)) 1 and the coefficients
2 (V1 + d 1)2
aa(T (x) -a(T tx+1)
w
2 tanh
4
/1 - d b 20 202 da
VP
Ok
c d²
To approximate SDEs kernels
24 and (4) on time interval [t,, f1. we may define the following transition
re Pk+1
1 ( (yk))
N(a VP yk + DYP stef 87-14 (yk); ck XPTa)
(34)
( N(a VP + by. ( (vk) + of (yk)}, PTT)
that can be obtained by applying the El scheme, Appendix B In this case, the log ratio log dpp /dy can be computed as follows.
Lemma 16. Assume that the joint distributions oe and pnd are respectively induced by the transition kernels defined in and 34). Then, for any yo.k E (IR))++, have
log ap dp (yo.k)
K-1
Σ
w
VP
19 -42 (2k)l1
+
K-1
Σ km0
V
-4- (Tx) c
where EN
yk++1
0g y/s
VP
byP rd
- (2x) + -4- (261))
Proof. The proof is exactly the same as Lemmal
10 noting that we still have wXP = //
By combining Lemma 16 with (2)). we deduce the expression of the RDS loss in this setting. The proof of the following is the same as in Corollary 11
28

---

Published as a conference paper at ICLR 2025
Corollary 17 (Expression of RDS loss El setting VP scheme). integrate SDE implemented with VP, the RDS loss simplifies as
When using the El scheme to
1 ( (86)
LRDS(0) Var
VP Wk 9T
(Yx)
2 2 (Yx))
K
}nd
+
W -T (Y&) T Zx + log (Yx)
VP
k=
H where {Zu LK-1 are independently distributed according to N(0, Ia), and {Yx is recursively
defined by Yo ~ buse, and for uny k E (0.....K
-
Ys41= = a Y + by 11 (Yk) o - (Y+) ZxTZz.
with 0 being a detached version of 0.
For completeness, we also provide a similar result for the KL.-based objective given in (26)
Corollary 18 (Expression of reverse KI. loss El setting VP noising). When using the scheme to integrate SDE (4) implemented with VP, up to additional constants independent of 0, the discrete time KL. loss simplifies as
VP &
(V2)) + log - ((RR) 7
(k)(0) = E
2
where {
K-1 are independently distributed according to N(O, Iz), and (VE) is recursively k=0
defined by Yo hure, and for any k E ,O.....K.
K. -11. 1}.
rd y = a Y Dg 18T-4
(YL) + 1 (Y + / Zx
D.4 DIS SETTING
Here, we propose discrete time alternative to the continuous time objective proposed by Berner ct al. (202)) and Richter et al. (2023). Let A e. Consider the variational diffusion process (Y, kejo,ry defined by SDE 21
(Alys) and (A,yk+1) P (Alu&+1) respectively approximating P(Yo, E E Y yyk) and P(X,, EA Xte. 1). Then, for K sufficiently large, the path measures ipo and (I** may be approximated by the joint distributions Ph:K E (K++)) and Po-k (K++1) defined by PO. = hase II Ph and
Assume that are given transition kernels (A,yx)
Po.k = π I -1 PRJK+1
Following the RDS methodology, propose to approximate the KL-based continuous-ting time DIS objective by
dp dp
C (0) IE log
(vo:))
Yoxx Pu-k
(35)
and the LV-based continuous-time DIS objective by
drl (x6x)) log dpp.
CDIS (0) Var
Yox ~ po.K-
(36)
To approximate SDES 21 and (I on time interval [tx, [--1] in the case of the VP noising scheme, we may define the following transition kernels
Pije+i (-1xx+1) = N(V
2 - La)
(37)
Ph+t1u(-1Bn)) = N(ax yu + 6hP + sT (yk),
0
VP
(38)
29

---

Published as a conference paper at ICLR 2025
where Ag.1 is defined in Lemma@with T and == = T- txp and {ax b )) are defined via the VP-EI integration scheme detailed in Appendix D.3.3] In particular, the following equalities hold:
= a² XT -T-1III-- = exp (a(T' (x+1) ([ff-41)]]
\T -4-11 -1-11 = exp
- th+1) - a(T -a(T -twyl (a
VP
1
VP (a
Therefore, we have Pick+1 (yx (k+1))
N YP yks yku cx La). In this case, the log-ratio
log dpp /dp' can be computed as follows.
Lemma 19. Assume that the joint distributions p' and pd are respectively induced by the transition kernels defined in 437 and (38). Then, for any yo.k E (RR)K+]
k+1,,,,, we have
log dp" dp-
huve (H0)
K
dp (3k+119n)))
(VO.K) = llog Z + log
log
(UK)
dp (3/ 13N++1)
k=0
π buse (yo)
K-
VP
K-1
= log Z + log where Eg
k
T- (1+))1++
Σ
V wk - (3k) Tex
vp
7(bK)
k=0
k=
(y&+1
VP ag yk
bk Sq
(yx))
Proof. This result lies on the following identity
dyp 1/r 114-11992
||Ex|1
1
cx + v V
YP s (9h)
1 (11+)4+1)
2
Based on Lemma 19 KL-based DIS losses.
we simply deduce the expression of the discrete time versions of LV-based and
Corollary 20 (Expression of LV-DIS loss). When using the El scheme to integrate SDE 21 the
LV-based discrete time DIS loss (36) simplifies as
K
CDV = Var
wk ST
vp
(Yk)
{st ( (Yk)
2 8T (Ya)}
1
k=
V /W - (Yx) Zu log hase (Yo)
7(Kx)
where
K- {== are independently distributed according Ia), and (Yr) is recursively defined by Yo baze and for any E (O...., 1-1},
Ys+1 0% YPYn 6% s s7 (Yv) VeT
with i being detached version of 0.
Proof. This result is application from Lemma 19 where have ey Zs + VWX (PP(GT- (Yx)
- (Yk)) for any k (O...., K-
Corollary 21 (Expression of KL-DIS loss). When using the EM scheme to integrate SDE (1 up to additional constants independent of 0, the KL-based discrete time DIS loss 35 simplifies as
K
UP k
T base (X,1)
CRiS (0)
E
9
(X
log
-1a
>(Y
30

---

Published as a conference paper at ICLR 2025
R==0 where {Zk] K-1 are independently distributed according N(0, Ia), and (YE) K
Tk=o is recursively
defined by Yo Thure, and for any k E {0. .....
1 1}.
y as Y + b/ 20 -u (VP) + VP 7x
D.5
CMCD SETTING
Here, propose discrete time alternative to the continuous time objective proposed by Vargas ct al. 2024). Let O. Consider the variational diffusion processes (Y?" (ET0)o)) defined by SDE
23 and (x )eE defined by (22).
Assume that we are given transition kernels (A, yk) + P1+11 (Alus) and (A, (k+1))) + &ekt11 (Alyk+1) respectively approximating P(X22 E A Yt = yk) and Q²(X;, E X 1)+1). Then, for K sufficiently large, the path measures IPº and (Q+) may be approximated by the joint distributions pb.K E (g(K+1) and g6.k @(K+1) defined by Po.k = n I P and K-1
base
G..K - I 0g k+11'
Following the RDS methodology, we propose to approximate the KL-based continuous-time CMCD objective
dp
CCCC E
log doh:
(Yo.k)
Yo.k ph:
(39)
and the LVV-ased continuous-time CMCD objective by
log dpb.
CCMC (0) =
LV
Var
(Yo.K) dg
Yu.k ~ Pu:
(40)
For any & E {O...., K 1}. let us define w
CMCD o² 8
To approximate SDEs (22) and (23 on time interval [t,, 11 we may define the following transi- tion kernels
026x (V log Teasi (yx+1) hi (2x+ )), 8z La),
(1 = N(y++1 +
(41)
2
0²
P1+1)1(-111) = N(yk + log Tea (yk) + hi, (yn)), 2 Ia)
(42)
that can be obtained by applying the EM integration scheme, Appendix B In this case, the log-ratio log dp" /dq² can be computed as follows.
Lemma 22. Assume that the joint distributions d and pd respectively induced by the transition kernels defined in 441) and (42). Then, for any Yo:k E ( we have
log dp" dg (Wo.k) = log Z + log
base (y1)
K
dp (1k++ l9k)
log
72K)
dgf |++1 (y/ 13k++)
k=
base (yo)
K-1
K-1
log Z + log
w
CMCD
||uell
w
CMCD un ck
km
where
Uk
V log Tra (2))
log Tezta (1++1) + h1g (yx) hl (H++1)
V log Tec (1n) + 12, (uk))
Uh++1
55
2
Proof. Based on the notation of us and ex. we obtain the result by noting that
dpk+11k (No++ 1 (yx)
log do (yn lon+1)
||cx|1
+ 2
1
c6
CMCD
ui
31

---

Published as a conference paper at ICLR 2025
Based on Lemma 22 KL-based CMCD losses.
simply deduce the expression of the discrete time versions of LV-based and
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
Tuuee (Ya)
K-1
LCMOD (0) = Var LV
log
W
CMCD
|
(1122 (Vk) - -he, (Vx)) Tuf
((K))
k=0
K-1
+
uuk
(ux) (
with ul log nr, (Y&) V log Tex (Yk+1) + (Ys) Mie.. (Ye41),
where 12x -1 are independently distributed according to N(0, Ia), and {Yu)kso is recursively defined by Yo ~ blase and for any E (O.....K.. -11,.
0²d2 Yk+1 = Yv +
log Te (Va) + 11 (Yx)) + +oVo
2
Proof. This result is application from Lemma 22 where have ca 7x +
u (1d (Yr)
h9, (Yk)) for any & E {O....,K K 1}.
Corollary 24 (Expression of KL-CMCD loss). When using the EM scheme to integrate SDE 23 up to additional constants independent of 0, the KL-based discrete time CMCD loss simplifies as
π hase (YO))
K-1
KL CCMCD(0) = E log
CMCD is
|
k=U
with usg V log Tee (Y V log Thess (Y + hi (YF) h2 (Ye+1).
where {Za K-1 are independently distributed according to N(O, Ia), and (Y Ik-o is recursively I== defined by Yo hhvee for any k E (O....,.. K-1}.
σ20x (Vlog Te, (Y + 11, (Y 1h,
y x==1= = Yp +
o
2
E ANNEALED MCMC SAMPLERS
In this section, present line of Monte Carlo algorithms that aim sampling from sequence of & distributions {px} ) k= C P (RR),, with K 2 1, defined such that Po is an easy-to sample distribution and pæ is harder and harder sample as k increases. In this section, we assume that each distribution Pa 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-Adjuste Langevin Algorithm (MALA). This sampling algo- rithm is the only one in this section to be designed so as to sample from single distribution E (Red) with density 7 Rd R2 known up to normalizing constant. It is an extension of the largely used Unadjusted Langevin Algorithm (ULA) (Roberts Tweedie 1996).
Given number of steps L > 1, step-size À > O and an easy-to-sample distribution pinit E 90 (RR), ULA builds a Markov chain (X ) defined by X ~ pimit and for any & E (O...., L 1).
X = x AV + V2AZ' z' N(0,
which reduces to a time discretization of the Langevin diffusion that admits π as invariant distribu- tion. Building upon this recursion, MALA (Roberts & Tweedie 1996) addresses the discretization bias by adding Metropolis Hastings acceptanceenejeetion step at each iteration, Algorithm 5 In practice, one can easily adapt the step-size À by targeting a specific acceptance ratio. We follow this paradigm in numerical experiments, Appendix for more details.
32

---

Published as a conference paper at ICLR 2025
Algorithm 5: Metropolis-Aduuste Langevin Algorithm (MALA)
Input: Target distribution ı with unnormalized density 7 step size À 0, number of steps & 2 1, initial
distribution y
P Initialization
X
for
E {0, L- 1 do
> Propose new sample
x = X + XV (x*) + V21Z', z' N(,,
> Compute the acceptance ratio
(1, >(X jox"
min
where q(v x)) = N(y; AV (x:), La)
> Accept or reject depending on a
U(00, 1J)
U S a then
x X
clse
L x X
Output: Soquence of samples X0 such that X lid (approximately)
Annealed Langevin MCMC (AL-MCMC). Based on exact samples from Po. the Annealed Langevin MCMC algorithm (Song & Ermon 2019 2020 consists in sequentially sampling from Pk-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 casy-to-sample distributions (i.e., when is low) can help to overcome sampling issues, such as high energy barriers, in the more complex distributions (iee,, when le is large). We refer Algorithm 6 for the pseudo-code of this approach.
Algorithm 6: Annealed Langevin MCMC (AL-MCMC) algorithm
Input: Sequence of target distributions {px K steps per level L 2
Jm step-sizes [Me] E (0, o0)x number of MALA
Initialization
X ~
for k O.....K. K- do
b Sample from Pk++1 by starting at samples obtained from Pk
K+1 MALA(PL41, A L,xxp) xxx) Algorithm
Output: Sequenee of samples X such that x ud Pv (approximately)
Sequential Monte Carlo (SMC) and extension. The Sequential Monte Carlo Neal 2001
Del Moral et al. (006) algorithm lies the same paradigm Algorithm proceeds in parallel over N chains (which are referred as particles) while ensuring that the warm initialization is correct. additional resampling step that may occur before certain MALA run,
To do so, SMC relies based on accumulated importance weights that we next define.
Let N 2 1, k E t0...., K 1}. (wi) V be unnormalized importance weights accumulated on each chain up to step h (wy 1/N if k 0), and (xx)) be samples from Ps obtained independently as the last samples of N parallel MALA runs (or exactly obtained if = 0). Then, the SMC importance weights designed to sample from Pk+1, denoted by (wi (wi++))m and their renormalized versions, denoted by (wF+)) are defined as
Pk+1(x2) (x2)
W K++
W
Wh Wh+1
(x
Σ W
The SMC resampling step then consists in obtaining N new samples (yE+)) D by random selection among the original samples (xg) via multinomial distribution (with replacement) defined by the normalized weights (wi+1)m 1 These novel samples will then be used starting points for
33

---

Published as a conference paper at ICLR 2025
running MALA the target Pk+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
.V (Em= m= Wk
ESS₁
Σ rumm
(wh)
We detail the whole SMC procedure in Algorithm
Algorithm 7: Sequential Monte Carlo (SMC) algorithm with adaptive resampling
Input: Sequence of target distributions {px 1K step-sizes 18 E (0, oo) number of MCMC steps per level I, 2 1, number of particles N 2 1, ESS resampling threshold a E (0.
Initialization iid
Po. W 1/N XD
for k
O O.....K K - do
> I. Compute the SMC importance weights (in parallel, for ,,...,N))
(xx))p((X()) W > 2. Compute the accumulated ESS (in parallel, for n = I,...,N))
W =
ESS (En=) WF) (EZ)))))))
>> Resample X if the ESS too low
ESS < aN then > Normalize the importance weights (in parallel, for I,...; N)
WT+1 / ΣØ= W k41
> Resample the particles
I I:N ~ M(W M(W c+1...., W
M x
WEN 1/N
else
Yx+1 X
> 4. Sample from pa+1 by starting from sample Y (in parallel, for n = I,.....N)
11 = MALA(pe+1 AE-1; L, L,8yg.1 see Algorithm 5
X
X
Output: Sequence samples X1-8 such that X lid
pk (approximately)
Recently, Phillips et al. (204)) proposed extension of the SMC algorithm, Particle Denoising Diffusion Sampler (PDDS). where they consider the specific case where the intermediate distribu- tions {pk) k= are defined as the approximate marginals of denoising diffusion process. Below, we explain the fundamentals of PDDS.
K
Let E (RR) be a hard-to tord-to-wae to-samlle distribution. With the same notations as in Section 2.1 consider noising process (X* ")eeo,,y induced by the following SDE
ax; = d() (etedd+B((())l)) VB())W,, XØ 7
(43)
Then, assuming that X base the time-reversal of SDE is given by
--T t) 3(T t) log log B(t)d Y huse
(44)
where pt the density of the marginal distribution associated to X². which intractable in gen- eral. Assume that we are given a sequence of intermediate unnormalized densities (p<)ue(o.r that approximate the target densities (pī )ejo,
In partieular, the time-reversed SDE (44) can be approximated by
dy, --( t)t + ((T (Yi)p + VB(T- - f)dB, Yo n base
(45)
Consider a time discretization { such that b 0 and tx = T, and the sequence of target distributions {pe { defined by Ps = pr- such that pr. is expected to harder and harder to sample from l increases.
34

---

Published as a conference paper at ICLR 2025
linear, it can be exactly Let & .O.....K 1}. define dx tx+1 tg. Since SDE 43
integrated on time interval [T t tgl, see Lemma We denote by P the correspond- ing transition density. On the other hand, SDE (45) can only be approximately integrated on time interval [tx, tx+1] (due to the score in the drift). For instance, this can be done using the transition kerne
Pk-1jk (1ya)))))))))))))))) Ny f(T tu)ŌRyk + 3(T t4 logpT- (vx), -
To sample sequentially from (p=) Phillips et al. (2024) suggest apply the SMC procedure detailed in Algorithm 7
where the k-th resampling step (1-3) is now defined as follows.
Sample Xk+1 Pk- ( X (in parallel, for n 1,..., N);
1. Compute the importance weights (in parallel, for I....,, N)
Pk+1 (X k+ 1+1)
(X k
W3:
A+1
Pk (x P X
xx)
2. Apply the resampling procedure (based on the ESS value) to the particles X weights given by (WT+1)m)
with importance k+1
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 sampling from the annealed distributions { with K + 1 parallel Markov chains, i.e., the k-th chain targets Ps. For the sake of pedagogy,
assume here that K is even.
At each level, RE sampling is done via local MCMC sampler, such as MALA, which is expected have poor performance for large k. To alleviate this issue, every L MALA steps, RE randomly performs swapping between Markov chains. More specifically, each consecutive pair (k, k + 1) of Markov chains (i.e.. either (ks k+ E {(0, 1). .... (K-2, K --)) or (k, k+11 E ((1,2).... 1, K)) K)), with respective current states X and X is swapped with probability
Pk,++1 = min (1, (Xx x 1
Pı (XA (Xx+1)
We provide the pseudo-code of this algorithm in Algorithm
Algorithm 8: Swapping step of Replica Exchange
Input: Sequence of target distributions {pk k with even K, current samples XI.K, swapping state 8 (0,1)
if = 0 then
Indexes = (0,2,..., K- - 2}
else
L Indexes (1,3....,
Compute the swapping probability between & and & 1 (in parallel, for Indexes) (P
Pe,++1 min (1, Pr((xpix (1]
& Swap depending on P4,641 (in parallel, for k Indexes)
U ~ 1])
U < Ph.k+1 then
L (Xe, (X++1.
Output: Randomly swapped sequence Xo-K
'In the original paper, the authors consider the VP noising scheme and compute P4 the transition
kernel obtained Exponential Integration of SDE 45 [fe, f++1]-
35

---

Published as a conference paper at ICLR 2025
Algorithm 9: Replica Exchange (RE) algorithm
Input:
Sequence of target distributions {ps J with even K, step-sizes {4 E (0, oo) total
number of MCMC steps M 1, swap frequency $ 2
Initialization lid
.o.
for m
................................ M do
(m + mod == then
& Swap Markov chains indexes determined by
x. = swapRE( (X)...s) see Algorithm
(s + 1) mod2
else
Sample locally (in parallel, for ..... K) m++1
MALA(pe, Ae; 0xgx) see Algorithm
Sequence samples Xa such xg. idd pk (approximately)
Output:
Known limitations of annealed MCMC samplers. In general, annealed MCMC samplers usu- ally target sequence of distributions { which forms geometric bridge between an easy- to-sample distribution huse and sard-to-sample distribution म. i.e., Pk (T (K/) ///K However, this scheme suffer from mode-switching for multi-modal distribution и (Phillips et al.) 2024]. Le., the mode weights of the intermediate distributions may vary fore, transitions between annealing levels is hard in practice, samples can easily get stuck in high-probability modes. Moreover, due to the use reweighting schemes, SMC and RE are highly sensitive the overlap between consecutive distributions, which requires to take K large in complex multi-modal settings, hence increasing the complexity of the methods.
increases. There-
F ENERGY-BASED MODELS
This section is dedicated to describing mhergy-Based Models.
F.I SINGLE-LEVEL EBMs
Let p E (RR). We assume having access to samples from p but not to its density. Energy- based Models (EBM) perform density estimation task by modeling the density of p with (z) exp(-E*(()) Zz, where E : Rd +
is an unknown normalizing constant. This model is trained by maximizing the log-likelihood of the model ((p) = Ellog pd (X) where X p, whose gradients can be computed as
neural network and Z = fRv exp(-E*
v.D(p) =
-V,,,""()) X~
X+ ~p~p.
(46)
Given samples from p (referred to as positive samples) and samples from p (referred to as negative samples), can compute Monte-Carke estimator of (yp(()) using (46). In practice, the negative samples are usually obtained by running a MCMC sampler on pi (which is possible because it is known up to a normalizing constant). However, obtaining truthful negative samples the main challenge of the EBM training procedure as pd is expected to be multi-modal as p.
F.2 MULTI-LEVEL EBMs
Let ((K+1), with K 2 1. For = L...., K, we assume having access to samples from pk but not to its density. Regarding Po, assume to be known entirely. Our goal is to approximate the density of Pı at every level E f1....,AA K}. A naive approach would consist in defining single- level EBMs for each k. Assuming that the distributions Pk are not completely unrelated, doing this would be however very inefficient. This section presents alternatives to this approach, which are all implemented in code.
36

---

Published as a conference paper at ICLR 2025
F.2.1 PRIOR WORKS
Diffusion Recovery Likelihood (DRL). Recently, Gao et al. (2021): Zhu et (2024) proposed the DRL framework, based on the extra assumption that for any k E (0. ....... K - 1]. the condi- tional distribution Pijk-1 known and given by P&|a++ (- 2k+1) = N(a++174+1, 0 Ia) for some ak+1 E I and 0k+1 > O. This is typically the case when the joint distribution p is de- fined via a discrete time approximation of a linear SDE, see Appendix C In the following, we denote Yg= Xx where X Ph- In DRI, the authors suggest to use a multi-level EBM p = exp(= (h: 2/ where EYP (1...... K) x R + IR. is a neural network and the normalizing constant 2% = -sg cxp(-E" 2: 2))dr is intractable. approximated, up to a multiplicative constant, by exp(- (B,(,,)) Using Bayes rule, they define for any k E {O,... ......K. K 1} the following conditional EBM
this case, the density of ps
Pi+1n(Uk+ æx) x PKA (Irk /1+1//04+11 1)p (Uk+1)
= (Z (k + [x
exp
| 3b+11l2 E 1. shott)
20
dy. Note that if 02 is small enough, for instance when K is large, then the quadratic term dominates in the expression of log +1je .g), which means that this conditional distribution localizes on Th- Therefore, by ap- plying a first order Taylor expansion of the energy term zg. the following Gaussian approximation holds
where ( + 2,,)) fRr= exp (20
||2
yll Ex + 1,9) +1,)
Pk+ (vk+1 (x)) N (Ix of V + 1, o
Hence, if T 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
K-1
CDRL (4) = Σ (PR)(p), with (p)) = [og Pi (ax+1 Xx) k=0
(47)
where X Pk+1 and Xx ~ Pk|++1 1 (11X++1) In particular, the gradient of the single-level DRL loss defined in (47) can be expressed as
Vec (x) EP EP (ki Y - E Xt)]
where X₃ ~ ph. Y P (- Xx) and X pu+1- Up multiplicative constant, the density of Ps can then be approximated by xy + ph (o)
Algorithm 10: Gibbs-within-Langeee MCMC transition kernel
Input: Previous state (k),z) {O..... K) x TRd, step-size À 0
p Sample ē given with standard MCMC
MALA(pp(( k)), X) where p²( k)) x p-))
Sample & given म with multinomial distribution
h~ M (exp(-E* (0, (0.))) #)))
Output: Next state (k, F)
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 lule EBM expressed as pf (z) exp(-E* ZZ,, where E {O,.... ......,K K} x gd + R is neural network and the normalizing constant g = Σt=o.f See exp(-E* (k:, x))d, which can be learn by maximizing the following maximum likelihood objective
CDALBM/ (6) [log pi (Xxx)]
k ~ U(FO,.... K}), X Pk
(48)
The gradient of objective (48 can be written as
V CDABBM
E E (KT (A-, - x
where (k- X pt, A+ U(1O. K}) and X Pk+ In the same fashion Kim & Ye (2023), the authors suggest to perform the negative sampling from pd by doing a Gibbs withim- Langevin procedure, whose transition kernel is summarized in Algorithm T0
37

---

Published as a conference paper at ICLR 2025
F.2.2 OUR APPROACH: USING AN ANNEALED MCMC SAMPLER AS BACKBONE
In this paper, we also suggest learn a multi-level EBM approximating p. defined by P (2) exp (-E*(kxx)) / where E {0,..., K} I + I neural network and the normalizing constant 2 Jaa exp (k) dr. We propose to learn yp by maximizing simple yet novel maximum likelihood joint objective
K
oours
(p))
= (p) = E llg pf (X+)., X Pk.
(49)
O
with single-level gradient given by
V cours (p) = -
X5 ~ p X ~ Pk-
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 En the sequence of EBM densities (p: )) This allows us 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 default annealed MCMC sampler because of its massively parallel capabilities. We summarize our training procedure in Algorithm| [1] 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 traing iterations N 2 1, initial parameter you batch size B 2 number of MCMC steps L 1, predefined annealed MCMC sampler
for 0........ . do
b Get samples from Pa (in parallel,
0,...,K) K)
{X k,6
1 lid
p&
b Get B3 samples from each p using annealed MCMC sampler
(x annealedMCMC (ipl" "1p"1K fk B) > Apply a stochastic gradient descent on şon
BB.K
Compute the MC estimator 9n of the gradient of the loss L defined in 449
9.. = B E (Ei-, V.EEFn (k, Xx.) Σi= E. (k, XT,)
Update fh. to 9 with Adam optimizer based on 9.. Output: Optimized parameter 1pN
Multi-level EBM parameterization as GMM tilting. To reduce the computational footprint of multi-level EBMs in the context of LRDS, suggest define our EBM model tilting of learned GMM. Let ş (wy. my. Σ, 1 be the parameters learned by GMM-LRDS (Algorithm
to approximate gel and denote by of the f-th marginal density of the resulting reference process. Building on this, we define pi tilting of ps by
log P. (x) log p (xx NN
where NN [0, T'] X Rd + R is a neural network with parameter p. In this scenario, we design the initial (Po to ensure that N 0. In practice, in the case where the GMM shares very few similarities with re close to f = O, we rather consider
log (x) = Lextu (1) x log pf (x) NN
where 4mm E [0. For most of our numerical experiments, we found that plim = 0.2 brought satisfying results.
38

---

Published as a conference paper at ICLR 2025
G ADDITIONAL METRICS FOR VARIATIONAL IEFFUSION-BASER 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 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 2 1, these performance criteria require the evaluation of variational importance weight function variational approximation, that is specific to each variational setting.
(IZ)) + R,. assessing the quality of the
For RDS (including PIS and DDS), we define
zel dpo. (Ho.k) (No.K (
w(yok) = 2 dy (yo:k)
dp (yo:k) )7(k))))))))))))))))))))))))))
where p" E (K+1) and prf E (K++)) are respectively defined as discrete time approximations of the path measures ppe and (poety R. Appendix D.3 for more details.
For DIS, define
w(yu dp (Vo:k)
(po.k(y):k)
where p² E (1) and p" (K++)) respectively defined discrete time approximations of the path measures IPo, induced by SDE (1)) and (P*) see Appendix D.4 for more details.
For CMCD, we define
dp (Bo.k)
w(yo:k) = Z dg (Vo:k)
where pd E (((K)) and qa E (K+1) are respectively defined as discrete time approximations of the path measures 10 induced by SDE 233 and (Q where Qo is induced by SDE (2) 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, Z simplifies in the DIS and CMCD expressions. This also stands for g in the RDS setting.
G.1 RREVERSE' PERFORMANCE CRITERIA
Definition. al. 2024).
Standard variational metrics, referred to as 'reverse' performance criteria by Blessing defined as expectations with respect the variational distribution.
(a) The Evidence Lower Bound (ELBO), expected to be maximized, defined by
ELBO-
-Bloo w/ (Yo.)) Yo.k ~ Po.K
For any RDS setting with intractable gıef, 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 2, expected to be close to Z, defined by
2, = E (w((G.))) Yo.k ~ P
For any RDS setting with intractable 2re this criterion cannot be used. (c) The normalized Explained Sum of Squares (nESS,)) expected to be close 1, defined by
E E[ (Y )1 ([w U.K )
DESS
Yo.k ~ Po.
39

---

Published as a conference paper at ICLR 2025
Computation. Based on the results from Appendix D.2 detail the simulation of Yo and the computation of log w( o:K) for each setting.
RDS setting: Yo x hase and for any l E {O...., ......KK- -1-1,,
+ = agYK ( t4 (Y,") + -7 42 (Vh)) + volka Zn ~ N(O,Ia),
and the importance weight function is given by
K
K Σ
(VP) Z5 log 70e (YR)
wk log w(Yok) Σ
1
(22)
2
where {wk, ags bb: Ck
K- depend on the noising scheme and the discretization setting. k=u
DIS setting: yo hasee and for any & KK.....K -
Y Y b (X2) VolZa. Zn ~ N(O,Ia)
and the importance weight function is given by
K-
π hase (Y
log w(Yokk)) Σ
wk
0 TT
(12)
v Vurs 0 (82) Zz + log
2
>(gg)
R=0
where {wk, aks by, ck k=u come from the EI discretization of the VP noising scheme.
CMCD setting: Yo ~ T hase and for any k E {O....,K' 1}.
bh V log Te (YL) + hlg (Y + JoZza Ze ~ N(O,La),
r641 r + and the importance weight function is given by
2
K
K-1
base (Y
log w w(Yo.k) =
Σ
wk
CMCD
||
+
w
CMCD (u²) ZE + log
7(YR)
k=0
k=0
with ux V log Te, (V:) + V log Teat1 (Y++1) n6. (Y n1 (Y
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 'orwwrdd' 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
EUBO- E|log w wYo.k)), Yo.k Pb.k
For any RDS setting with intractable zref, 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 2,. expected be close to 2, defined by
2, = (E(w<Yo.k)~1))
Yo.k ~ Pó.k
For any RDS setting with intractable nee this criterion cannot be used.
(c) The normalized Explained Sum of Squares (nESS,), expected to be close to 1, defined by
nESS, (E(w<Yo:k)) ((B[x))))][[((xxxx)))) Yo.k Pó.K
40

---

Published as a conference paper at ICLR 2025
Based on the results from Appendix D.2 detail the simulation of Yo:k and the computation of log w (Y..)) for each setting.
Computation.
RDS setting with EM scheme: Yx n, and for any (0....; -)
Y = St+1(t)) + Sk+1 (1s (R+ 2x Zx ~ N(),
where Sk(t) = exp (f- f())du). o2(1) o = - Se du, and the importance weight function given by
K
{ (Tr
log (Yox)))) =
Σ
Wk JT - (Yx)
EM a
88 -1 1 (Vt) + 2 20r-04 (Ya)
k=
K-1 Σ (t) (Yt) (x) Zx
K-1
Σ
+ f(T tx)ok
gt (Yi)) TY + log no (Yk),
Sk41(tz)
as consequence of Lemma
10
RDS setting with El scheme (PBM): YK ~ n,, and for any h E KK...., 11. 1}.
PBM
Y
Y
Zk Zz ~ N(O,Ta),
PBM
and the importance weight function is given by
K
w - (Yx) { (Yk) + 1
log w (Yo.k)
PBM
2 2-- to (Ya)
k=
w 9T (Yx) Zx log se (Yk)
PBM 0
7
k=0
as a consequence of Lemma
RDS setting with El scheme (VP): Yx ~ n, and for any k E (0...... .......K - 1}.
Y
Zz Zx N(O,Ia),
VP
VP
4
0k
and the importance weight function is given by
K -
(Y)) + 1 -0 (Y+)
log w/ (Yo.k) =
wk 91-44 (Yx) (t (Yr)
VP
k=0
K-1
w --- (Yk) Zx log n00
y (Yk),
k=0
as a consequence of Lemma 16
DIS setting:
Yx cu n, and for any k {0...., K- 1}.
Y++
V ck
Y8
+
Zx 7/ N(O,Ta).
VP
VP
dg
and the importance weight function is given by
K-1
K-1
base (X))
log w/ (o.k)
Σ
W
-4
(Y+)
Σ
V
w ST (Vx) Zx + log
VP 0
n:
(Vx)
k==
0=
consequence of Lemma 19
41

---

Published as a conference paper at ICLR 2025
CMCD setting: Yк n,, and for any (o...., K -1},
Yx = Yl+1 + o² bk (V log 1 (Y+11) - Mi.s (V++1)} + Vo 2x ~ N(O,Ta).
and the importance weight function is given by
K 2
K-1
huse (Yo)
log =
CMCD
|
Wk (OCCDD( T Zs + loo
k
>(xx)
with ug Vlog r (Yk) + V log Teat (Yk+1) + h1g (Y&) h (Yx+1),
as consequence of Lemma 22
H IMPLEMENTATION DETAILS
H.1 DETAILS ON TARGET DISTRIBUTIONS
Gaussian mixture from Figure 2 The distribution from Figure 2 16-dimensional mixture between N(m,, Σn) and ((m Mg)
with weights 0.7 and 0.3 respectively. Let 1a be the d- dimensional vector with all components equal to 1, then we have
mi = 0.66 x 116. ma = +0.6 x 116
Moreover, denote Rg rotation matrix with angle 0 between the first and last axes, then
Σ = Rx/4 diag(10' 110-
10 , ))= E₂ = Rz/6 diag(1001 1000 ...., 10 2
In practice, we estimate the weight of the first mode wn by computing a Monte Carlo estimator of
I Ni) (I))
(50)
Bi-modal Gaussian mixture. This Gaussian mixture is used in Figure also in Section is d-dimensional mixture between N(-1a, Σ) and N(1a, )) with weights 2/3 and 1/3 respectively. The different values of 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 number of mixture components
2, we define d-dimensional Gaussian mixture x E IRd ~ Ete waN(<;; my, 0.5 where the mean locations {me)/.1 are independently distributed according to U ([-L, Lid) and the mixture weights [wr} are strictly increasing, defined with a constant geometrical increment such that wL/ws /w1 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 TR+ 10, 2 be the transformation from polar cartesian coordinates. Let 0 ~ U( 2T]) and R ~ N(r, with r, a > O, we say that (RR,8 0) is distributed ring of radius and width σ. The Rings distribution a mixture of rings with radiuses and width a = 0.1. The weights of each ring within the mixture are 1/9, 3/9 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.
1,3,
Checkerboard distribution. We divide the square [-4, into 16 squares of size x 2. This distribution is defined as mixture of uniform distributions on those squares in interleaving fashion. We assign a weight 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.
42

---

Published as a conference paper at ICLR 2025
Table 3: Covariance settings for the bi-modal Gaussian mixture. We denote by logdiag b) with 0 < < b the diagonal matrix in Rdxd whose diagonal is a log-linear interpolation between and b. Moreover, Qa & Raxd is a random orthogonal matrix built by doing the QR decomposition of random matrix distributed according to U([0, 5](xd). The seed used to generate this last matrix fixed at 42.
Type Diagonal Diagonal
Difficulty
Cond. number
Covariance
Isorroicc
1
(5 X 10-
Medium
102
(5 x 10 10-2)2 x 1ogdiag(10-3.1
Diagonal Hard
10a
(5 x 10 10-2)2 x logdiag(10 (10-1.1)
Full
Medium
102
(5 x 10-2)² Qa logdiag(10 (10-2,1) Q7
Full
Hard
10a
(5 x 10- Qa logdiag(10
Q4
The o² field system. The o model is a simplified framework often used as 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 oie-dimensional grid with size d 32. Each configuration is represented
d-dimensional vector (d:) To further
constrain the system, we the field to zero at both ends by setting Φo daa1
The negative log-density of the distribution is then expressed as follows
In (T) -3 ad d-1
Σ 0 2.1)
4ad Σ (1-07)²
1
hop
(51)
We selected parameter values that make the system bimodal, setting 0.1 and the inverse tem- perature 3 = 20, while adjusting the value h. We define W4 as the statistical frequency of configurations where 6a/2 O, and as the frequency where da/2 O. When h O, the system symmetric under the transformation & + Φ. so we anticipate w/4 = w- However, for > O, the negative mode becomes more prevalent.
For large values of d, the relative likelihood of the modes can be approximated using Laplace ex- pansions at the Oth and 2nd orders. Letting di and df represent the local maxima of (51). these approximations provide the following results respectively
(Øth Th (Ø4)
(x() ) x det Hy, (1 Th (04) X det Hí, (04))
(52)
W+
W4
where H2 is the Hessian of the function @ - In IT. (Ø). We compute the modes df and df using gradient descent on the distribution's 's potential. Since we do not have access to ground truth samples in practice, we compare the Laplace approximations defined in (22 considered our ground truth, to Monte Carlo estimation of w. w//w..
Bayesian logistic regression models. Finally, we evaluate the performance of Bayesian lo- gistic model, defined for any pair (s. y) E RRdim x {0, 1} by the likelihood p(y|æ; w, b) = Bernoulliy;; a x + 6)), where w idim weight vector, 6 E
an intercept and the sigmoid function. Given a dataset D = { 1y) C IRdim x {0, 1] of size M, we aim to sample from the posterior distribution p(w, b/D) x p(w, b) II p(y, w, b) where p(w, 6) = p(w)((())))))) (bb) a fixed prior distribution. Following Blessing et al. (2024),
h. we consider four real-world settings of binary classification problem: Ionosphere (dim 35. M 351), Sonar (dim 61, M 208), German Credit (dim 23, M = 1000), Breast Cancer (dim = 31, M = 569). Each of these datasets is randomly split into training subset Duaain of size 0.8.M and test subset Diest of size 0.2M. In this setting, we define the target distribution श = p(w, b/ d/Dimin) and evaluate the sampling quality by computing the (unnormalized) average predictive log-likelihood log p(w, b(Dues), which expected to maximized.
43

---

Published as a conference paper at ICLR 2025
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 is as expensive computing V log æ. We assume an infinite parallel computational capabilities. Note that the cost of the training procedure can be amortized sample multiple times.
Method
Number of target evaluations
Number of neural network evaluations
SMC
O(KM)
0
RE
((M) O(NK) O(A(N + KM)) O(AN)
PIS/DDS/DIS/CMCD
O(NK) OAA(N + KM)) (AAN))
PDDS
iDEM
LRDS
O(N)
O(NK)
In each case, the prior carefully chosen as a Gaussian distribution with the following parameters:
lonosphere: p(w) = (0, 5.25 la1) and p(b) = N44.25, 0.25²),
Sonar: p(w)= N(0, 4.5 Teo) and p(b) N( ( -25, 0.5²),
German Credit: p(w) = N(0, 1.25 and p(b) - N(3.25, 0.5*),
Breast Cancer: p(w) N(0, 3.75 130) and p(6) N331,
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 the target distribution than a GMM approximation. The initial point of the MCMC samplers 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 dis- cretization 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 Note that, in LRDS, the complexity of MCMC sampling (to compute (nd) or the complexity of learning the reference process were ignored
as they are negligeable compared the training and sampling budgets; see Appendix H.3 for details.
H.3 HYPER-PARAMETER OF EACH SAMPLING ALGORITHM
In this section, we detail every the results in Section The computationally-rreee parameters were chosen to ensure a comparable exccution clock-tima between the different algorithms the same hardware.
for each algorithm involved in the computations of
Construction of the MCMC dataset for fre We build the dataset of reference samples by run- ning MCMC samplers initialized in the modes of the target distribution. We run 4 chains mode. We use the Metropolis-Adjuster Langevin algorithm (MALA), see Algorithm pler, except for the Checkerboard distribution where (RWMH) algorithm as the score is constant. In both cases, we adapt the step size automatically for 8192 warmup steps aim a 70% acceptance rate. The datasets 60000 samples long.
as default sam- the Random Walk Metropolis Hastings
44

---

Published as a conference paper at ICLR 2025
Annealed MCMC methods (SMC, RE). For both algorithms, incorporate prior knowledge into the base distribution by taking π huse = [(m, Σ) where m and Σ are computed using maxi- mum likelihood samples from nh In each case, the annealed density path is a linearly-spaced
geometric interpolation, where Th x ( / 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
= 1. For RE (see Algorithm 9.. use K = 128 annealing levels with 256 independent chains per level and a swapping frequency every S = 8 steps for 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 Riichter et al. (2023), while we use the implementation provided by Blessing al. (2024) for DIS and CMCD. For PIS, DDS and DIS, set the arpar-prraette as advised by Appendix 1.2 by computing the mean and the variance of the samples from 사 to add prior knowledge. For CMCD, we sample a version of ı which is shifted by m, defined the mean computed on samples from and,, and we set a 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 Fourier MLP, as in Zhang& & Chen (2022). with 4 layers of width 64. In the case of PIS, DDS and DIS, we use target-informed parameterization by adding (N) x V log (xx (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 caareeprrametrizzationn in CMCD, since the drift of the generative process is already informed by n. We highlight that, by default, we do not use this target-informed parameterization LRDS, since the reference pro- cess 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 9o O. 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, use the implementation pro- vided by the original paper Bakhound-Sadege et al., 2024). In order to provide prior information leveraging gref, (i) we standardize the target distribution by using the empirical mean and variance of an (ii) we preload the buffer with samples from भnर Following the design choices of the original paper, we use Variance Exploding noising scheme with geometric variance schedule o,- Because we standardized the target distribution, we decided to take = 5 for each target distribution to en- sure 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 f-t we standardized the target distribution by leveraging the empirical mean and variance of ree thus bypassing the original VI step of their algorithm. We also significantly increased the default computational budget by using 128 discretization steps, 2048 par- ticles and training for 100000 steps with a batch size of 512. The rest of the apperpparameters are the default ones.
LRDS reference fitting details. For GMM-LRDS, the EM algorithm is taken from Pedregosa et al. (201T). We use diagonal covariance for Rings distribution, Checkerboard distribution, di- agonal 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| II using the Replica Exchange annealed sampler (see Al- gorithm 99 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
45

---

Published as a conference paper at ICLR 2025
keep the negative samples but update the positive ones. We perform 32 MCMC steps to sample a batch of negative samples with a swap happening every steps. We keep 16 MCMC steps out of 32 to compute the expectations. To compensate the short length of our MCMC chains, 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 parame- terized as + NN r)) " as suggested by Salimans & Ho (2021). We leverage the GMM tilting initialization of the EBM perform exact MCMC sampling the very first gradient step. More- over, we preconditioned the network as advised by Karras et al. (2022) by leveraging the Gaussian approximation of gre
ADDITIONAL RESULTS
1.1 LINK BETWEEN TARGET AND REFERENCE PROCESSES
Let gee E (RR) be an arbitrary distribution. We recall that under mild assumptions on n² see e.g., Cattiaux al. 2003). the time-reversal of pref is associated the SDE
dyfs f(T tf dt B(T
V f)dB, Yfo Pp
(53)
Link with Schrödinger Bridge. Since Tp* and ppeef are based on the same noising diffusion pro- cess, these path measures are linked by the relation P = オ @ pref 10 Hence, using the KI. chain
rule (Léonard 2014). it is clear that P is solution to the following minimization problem over path measure space
arg (KL poo) PE E (CT) =
which is often referred to as a Schrödigrr bridge problem in stochastic control literature.
Link with Doob's s h-transform. This relation path measure space can also be written as IP* de peer, see Appendix D.1| Therefore, solely based on the SDE (53) describing (Iree)), R (IP*) dred may be expressed as Doob's h-transform of (Fref) R via the SDE
dY; = - t)dt + B(T +)(s (XG) + ( + VB(T - t)dB, Yo PT
where the Doob's control function h : [0, x Rd + Rd is defined, for any E Idd and any E [0, T), by
dm dr (Yr)) Y = V1
dr dapter (Xo) X, =
(x)) = V log E R
V log Egyer OIT
Additionally, we have for any E Idd
1 pi (x.) / P
plat (xt)
dm Pijo (1 (d
I dr (NO)
1t (5re)
dr dhyeel
Epst
(Xo) X1 Zt
By combining previous computations, we obtain h4 = V log / pp. for any t E T).
I.2 OPTIMAL SETTING OF ISOTROPIC GAUSSIAN REFERENCE DISTRIBUTION
The goal of this section is to explain how to set the hyperparameter o E (0, oo) in PIS and DDS set- tings to 'optimize' the reference distribution म.. = N(0, o² Ia). when targeting a multi-moda distribution. The same reasoning can also be applied in the DIS setting to set the base distribution of base chosen an isotropic Gaussian distribution.
Let E (RR) be our target distribution. Assume that we are given a diagonal Gaussian approx- imation के N(m, Σ) of with mean Rdd and covariance Σ diag(T2) E Rdxd where
46

---

Published as a conference paper at ICLR 2025
I' (0, o0)) We propose set nt 'close' as possible i in PIS and DDS by solving the following variational problem
arg (KL N(0, o² (a)) (0,
Let 0. We have
Ta
KL N(0, ('))))))))))))))))))
= log
- d+o -2 2m m Tr(o 2
d
log( o
Σ
log 12 - d -3m33mm m
==1
Moreover, the derivative with respect 2 is given by
ag² KL
N(O. rl)) = 2 do 0-1
m m r
In particular, the optimal solution of this variational problem is
o = dd
m mmm m + ΣE
which motivates us to set nd N(O, a Ta).
In practice, when do not have samples from the target distribution, rather rely diagonal Gaussian approximation of the empirical distribution fre (obtained with local MCMC samplers) and set ro N(O, o2 La).
Application to Gaussian mixtures. In the particular case where π is Gaussian mixture, we have an analytical formula for
Let ग Σj=1 WyNNm, Σj) be our target distribution with wi.r E As, where the j-th component has mean m, E I& and covariance Σ₁ diag(Aý) Rdxd with R3 E (0, co) Consider X ~ and Z M(W,...,,.). We first have
E[X] = E Z]| Σ wymy
We define the diagonal variance of by diagVar( = (Var[X] Var(X]a) (0, oo) where Var[X], [X, [] for any E {I...., d). By the law of total variance, we have
Var[X] = E [Var [X
Z1y1 + Var [ Zlk
wy (23)
Σ wy (m 2
By characterisation of Gaussian distributions, N(E(X), diag Var(m)) is the closest diagonal Gaussian distribution to π in the Kullback-Leibler sense. Therefore, in this particular setting, т may be defined by m E[X] and 2 diag Var(m). We use this approach when computing o in the numerical experiment commented in Figure
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 MMerroooli et al. 1953), the Metropolis-adjus/ee Langevin Algorithm (MALA)
Roberts & Tweedie 1996). the Hamiltonian Monte Carlo (HMC) algorithm (Duane et al. 1987;
Brooks et al. [2[1] or the No-U 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.
47

---

Published as a conference paper at ICLR 2025
RWMH
MALA
HMC
NUTS
21
30
3
300
Chain
Chain
Chain #
Chain #4
Figure 7: Samples from different MCMC samplers when sampling from bi-modal Gaussian mix- ture 2 dimensions There are different chains (in different colors). The MCMC samplers ran for 4096 warmup steps before producing those 1024 samples. We only display sample every 4 steps.
- -
- 8
J - 32
64
J J 128
Target
Figure 8: Increasingly expressive GMM approximations of the rings distribution (Top) Approxi-
mation using diagonal covariances (Bottom) Approximation using full covariances
1.4 LIMITATIONS 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, this case, getting a good approximation of the distribution requires significantly more components that the number of modes (here, 3). This highlights the limitation the approximation power of Gaussian Mixture models.
1.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 large diversity of settings. When ground truth samples are available (i.e., in all set- tings except Ø² and Bayesian regression), we display statistical estimations of the integral probabil- ity metrics defined in Appendix B namely, the regularized Wasserstein distance with regularization E = 10 (W2,)) the Maximum Mean Discrepancy (MMD) and the sliced Kolmogorov-Sminn Smimnov 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 modal Gaussian mixtures defined in Appendix H.I with increasing dimension d {8, 16, 32, 64). To assess the sampling quality at a global level, we compute for each sampler the absolute enror "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,
48

---

Published as a conference paper at ICLR 2025
Table 5: Results for bi-modal Gaussian mixtures with diagonal covariance and medium condi- tioning, averaged over 16 sampling runs (same setting as T 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 where results could not be obtained due to numerical issues. The MMD and KS results are displayed with 100-factor rescaling.
36
=2=22
Alascithm
Ert.
XBID
8S.
3DID
MDID
RS
SMC
11:45amm
1200
10.50.7.44
10 15
0.70cmm
18 18
1453
10
1.36a0
17.87
14.44-771
RE
18.3%
0:38aase
20.57emm
15.20am
18.03 444
0.20axae
19 20
14.55a..
17 46
20 20:364:::
13.66-1.0
LX-MS
6.8Suau
7.85uun
6.802
48.18
54.0
38
1.83au
29.3-amm
LN-DES
156
0.40saoe
31
16
0.77ax
21
76
16
16 650500
1.38amm
59 7177
15:00an0
LX-CMCD
12.25
53.09
20.15
78.07-117
45. 48.79us0
16.35
3.47amm
58 B8.15a2as
41.79
IDEM
83,35 mmm
a1.52
36,80uns
56. 7754488
4.26a445
50-021833
0..10.r..
117.52um44
50.1318-22
NIA
PIDS
0.40seve
1.86one
2 2.50
11:30rr
17
NIA
NDA
N7A
N/A
NIA
GMDA-LRDS
0.38aare L.75atoth
2.7
0.71axx8
4.1%oo
Table 6: Results for bi-modal Gaussian mixtures with diagonal covariance and isotropic condition- ing, 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 where results could not be obtained due to numerical issues. The MMD and KS results are displayed with 100-factor rescaling.
61
Algurihm
Made Fr.
MXDD 1
KS
Eis |
MMD 1
RS
Modic Brr 1
We.. L
XBID
RS L
15.2%
18.72-1/4
18:42:143
LiBusve
19,50.1.0
14.44.1
1941838
15,46.1.1
14.30.1.a
LX-PIS
1.96
2.63e
1.35
1.90ese
2.35mn.
28
1.88210
7.9"2mm
9.51.mm
LX-DDS
0.755...55
0.621800
0,56.0.7
1.50.0.30
0.8%
1.09.0.0
1.56.031
1951440
2.13.035
2.33:0.6
LX-DIS
16.45611335
10.29-1.3
14:15115:
16776137
1220e00
20.12:1.0
14.4..07
1007
1951448
15.70,00
14.55_0v6
LX-CMCD
6.503 amd
0.6
28
2.92a 49
1.7% 15
12tesue
2.
3.7onca
3.45am
08.756
2.2333 17
83.92-000
182561108
35.02008
07.17a1
NIA
PDDS
1:25s0r
3.500997
19658.00
2.000000
GMM-LRDS
1.65
010089 69
9.77gmm
2: 24
2.85
1.10eae
3.79
3.75
4.00
3.45
4.Menn
Table 7: Results for bi-modal Gaussian mixtures with diagonal covariance and hard conditioning, 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 where results could not be obtained due to numerical issues. The MMD and KS results are displayed with 100-factor rescaling.
Algerilhm
Mise Er:
Wark
SEID
KS
Mode Err
MMD
KS
Mode
3DDD
XS
SMC
15.00
17.00ee7 2n
16:56
0-47ax
10. m
13-712a4
13.98111
LV-PIS
8.451144
D:46.006
18.521559
12.45154i
10.090.0.0
1.251442
53,801.00
23.52.052
32.508
234139.
81.01115
30.401011
LX-DOS
26 73878888
0.32.000
25,80144
22.921729
0.781x42
41.971231
17.2011154
86.12 12
1.00aa:
43.64.3951
20.451
LV-DOS
23.00LE37
19.091x4
0.651xm
28.0111.
19.38_600
16.88 67
1,061092
22.04:14:
18.07
LX-CMCD
4.
02.331m44
NIA
1.34aa3
30022553
31:30-12a
41.18
41.33em611
NIA
PDDS
NA
N/A
NAA
NiA
N/A
N/A
N/A
NAA
GMM-LRRD
211511488
11.5.000
3.41.1.0
3,461050
1.P2
0,471E40
3.00.1:32
2.99.0.0
4.6%L.77
0.12...1
5.26.144
Table 5 corresponds to Gaussian mixtures with diagonal covariance and medium conditioning. completing the results of Table presented in the main paper,
Table
6
corresponds to Gaussian mixtures with diagonal covariance and isotropic conditioning,
Table
corresponds to Gaussian mixtures with diagonal covariance and hard conditioning,
Table
8
corresponds to Gaussian mixtures with full covariance and medium conditioning,
Table
corresponds to Gaussian mixtures with full covariance and hard conditioning.
For all of these bi-modal settings, we observe that the mode weight estimation error is consistent with the values of probability metrics (except for the regularized Wasserstein distance, which is not discriminative between the methods). In particular, GMM-LRDS is on par superior to all com- petitors in each setting, except the least challenging setting (" ("Isotropic' conditioning, see Table (). where LV-DDS is more performant.
49

---

Published as a conference paper at ICLR 2025
Table 8: Results for bi-modal Gaussian mixtures with full covariance and medium conditioning, 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 where results could not be obtained due to numerical issues. The MMD and KS results are displayed with 100-factor rescaling.
Alprritu
SMC
0.
21.56ee.4
16.7
18.8%
3.50
16.5
12.52
a0.5R
0.30
12.18exar
0.72aa.1
RK:
15.8523188
0:3739.0
21. 21.34t:34 14
15 15.80s141
16 20476a1.0
0:38
20,0311 58
15.3 15.5111:34 34
16.38623 25
0,700030
19.00s:34
15.0421.
LX-P0S
4.8%
0,20aa.0
8.80
8.0
26.4%
2.52...
28.21433:31
28.37
R1.184528
35.024139
27.85
1.6%0om
0:17am
2.68:0
338
0.3t2axe
1.83s m
4:27a9
0.75em
6.85x
5.021m
0.22
22:00:134
17.331044
10.37%a119
3.601802
21,651155
15.0018:3
18.794148
0.52.m
21.561181
17.18.6m
EX-CMCD
66.77
2.0010.9
85.07.144
61.421106
84,681853
02:58.1:39
05.77
6,80.431
85,011341
62.77111
PDDS
1.5%
0:
1.52
2.05
1.9
0.90a91e
2.63e1
3.09.8:77
0.7%
0.78e133
1.58
2 2.81
GMM-LRDS
1.5%
0.37
1.07.6.a5
2:36.0.3
1.8%u4.5
0.384940
2.05.sar
2.73:855
2.47
0,71.030
3,681139
8.21.66
Table 9; Results for bi-modal Gaussian mixtures with full covariance and hard conditioning, av- eraged 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 where results could not be obtained due to numerical issues. The MMD and KS results are displayed with 100-factor rescaling.
36
di 32
Alginihm
SMC
14225
D. 0.10aeae
17.60
13,85.11.0
10. 10.15a
0.24acu
12.0
+0.00
17.0
0.47a440
21.01a13
10 18.96.s.e
RE
16.15621
0.11aso
15 15.12a3
16 15
21 21.582311
15.99 33.81:113 32
16.35c1 mm
0:5729.0
30.25a1.0
15.30=1.0
EX-PIS
16.96411
15.78-1.m
1. 13.33
10
10.002979
30 756
14342
31.25om
1.25811.8
5.58-1140
5.901166
0.2uca
11.33:144
0821143
14:10.406
LV-CMUD
18.65
18.70.0%
36.4211.1
807.04
10.8518:2
6001112
65.88.000
IDEM
06.8
02211.0
84.50.0.2
60.95
11.2
96.75
0.4516.3
56.811..7
39.7111.9
19
06.75 JUM
10.3311.4
8LSTT
14.70.1.
PROS
N/A
NA
N7A
NAA
3A
N/A
NA
N/A
N6A
N/A
N/A
NIA
L85
2.55..100
3.22.6.4
1.0%..7%
0.25ace
2.22:091
2.00 02
2.2%
0.474xx3
9.13:1.7
Execution time. To demonstrate the applicability of LRDS in practice, we provide in Figure the execution clock-time of each sampling method for several Gaussian mixture sampling settings. The computations were all ran on the same V100 GPU. We notice that (i) the initial MCMC sampling step to build the reference dataset, then used initialize each sampler, and the reference fitting time from GMM-LRDS are completely negligeable compared to the training time, and that (ii) the training time of GMM-LRDS on par with previous variational diffusion- approaches.
16; type: diag
d 32; type diag
16; type full
= 32; type full
2000
Clock time (in s)
1500
500
15705
INtes
1838
1400
18406
GLOULRDS
Inference
Training
MCMC sampling
Reference fitiing
Figure 9: Execution clock-time for bi-modal Gaussian mixtures with diagonal covariance (left) and full covariance (right), averaged over the conditioning settings.
Ablation study on GMM-LRDS: effect of the reference distribution. To asses the robust- ness of GMM-LRDS with respect to the setting of the reference distribution, we conduct an ablation study that reveals the impact of modifying the location and the entropy of the refer- ence modes. To do so, we set d 8 and define the target distribution (2/3)N(c: 11, 0.05) + (1/3)N(z: 1 0.05), as detailed in Appendix H.l Instead of defining the reference distribution in GMM-LRDS propose to consider a flexible reference distribution given by the Gaussian mixture x E Rd H (2/3)N(z; m, 0.0 (1/3)N(r; m2, 0. where (m), m2) Rdd Rdd and 0.
E Rd +
Gaussian mixture fitted on MCMC samples,
50

---

Published as a conference paper at ICLR 2025
0.29
0.6
Mode Err.
40
0.15
MMD
W2.5
0.4
0.10
20
0.05
0.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
Perturb. intensity
Perturb. intensity
Perturb. intensity
Perturb. intensity
Figure 10: Results of GMM-LRDS for an 8-dimensional bi-modal Gaussian mixture when varying the location of the reference modes: the performance degrades as soon as the reference modes and the target modes further from each other.
0.4
0.3
0.8
Mode Err.
0.3
Wa.
MMD
0.2
0.1
0.1
0.0
0.0
Perturb. intensity
Perturb. intensity
Perturb. intensity
Perturb, intensity
Figure 11: Results of GMM-LRDS for an 8-dimensional bi-modal Gaussian mixture when varying the variance of the reference modes: except for small reference mode variance, the performance is unchanged.
2.5
0.2980
Mode Err.
2.0
0.2975
0.03
0.030
MMD
0.2970
KS
0.025
0.2965
0.01
0.2960
0.015
5
10
15
5
10
15
10
15
10
15
Mode gap
Mode gap
Mode gap
Mode gap
Figure 12: Results of GMM-LRDS for 8-dimensional bi-modal Gaussian mixture when varying the distance between the modes of the target distribution: the performance does not significantly degrades, despite the increasing sampling complexity.
To observe the impact of the location of the reference modes on GMM-LRDS, we take a = (same variance setting the target) and (m₃, my) N(-1a, oa N(1 lala) where om (0, 0.25, 0.5, 0.75, 1.0} is a mean perturbation parameter. Hence, the higher om, the further the reference modes might be with respect to the target modes. We display Figure the the results of GMM-LRDS for this setting, averaged over 16 sampling runs. This ablation study notably reveals the major need of precision on the location of the target modes (in practice, brought by MCMC sampling) to build the reference distribution.
On the other hand, to understand the effect of the entropy of the reference modes, we take (-1a, 1a) (same mean setting as the target), and {0.5, 0.75, 1.5, 2). which can
(m) m2) be seen as a variance perturbation parameter. We display in Figure II the results of GMM-LRDS for this setting, averaged over 16 sampling runs. Interestingly, GMM-LRDS works well as soon as the support of the reference distribution includes the support of the target distribution. In practice, this is verified by MCMC sampling.
ablation study demonstrates that
51

---

Published as a conference paper at ICLR 2025
Table 10: Results of GMM-LRDS for bi-modal Gaussian mixtures with diagonal covariance and various conditioning settings (d = 128), averaged over 16 sampling runs. The MMD and KS results are displayed with 100-factor rescaling.
Isotrpic conditionng
Medium conditioning
Hard conditioning
Mode Em: 1
War 4
MMD
KS
Mode Ert
Wa,
MMD 4
KS
Mode Err.
Wax
MMD
KS
6.0%
3.0aan 7.84a0.22
6.16amm7
6.8%
8.33e1 1.8000.00
0.2%6.3.%%
1:32erpp
8.08
Ablation study on GMM-LRDS: effect of the distance between the target modes. To asses the robustness of GMM-LRDS with respect the complexity of the target distribution, we con- duct a second ablation study that illustrates the behaviour of GMM-LRDS for target with higher energy barrier. To do so, we set (2/3)N[(z; a 1a, +
E R + 1a. 0.05), where E {1, 2, 4, 8, 16) indicates the complexity level of the target. For each target, we conduct 16 sampling runs of GMM-LRDS, and display the results in Figure [2] Our numerics show that the performance GMM-LRDS remains significantly consistent with increasing a, proving the interest of our method for complex multi-modal targets.
and define the target distribution π
Limitation LRDS: high dimension. Although LRDS outperforms its competitors when sam- pling from challenging bi-modal Gaussian mixtures, its performance tends to decrease when the dimension is large, like most samplers. To highlight this limitation, we provide in Table I0 the results of GMM-LRDS for various mixture settings with d = 128, where we observe that it fails to recover good estimation of the strongest mode weight (roughly, 10% of relative error).
1.5.2 MULTI-MODAL GAUSSIAN MIXTURES
We display in Figure 13 the results of all considered samplers, including LRDS, when targeting the multi-modal Gaussian mixtures defined in Appendix H.l with fixed dimension 8 and in- creasing number of modes L E {4, 8, 16, 32, 64). To assess the sampling quality at a global level, we dispense for each sampler the total variation distance between the true weight histogram and the weight histogram obtained via Monte Carlo approximation. All metrics are computed based on 8192 samples. We notably observe that GMM-LRDS has the best performance compared to other sampling methods, independently of the number of modes in the target distribution. In particular, Figure 13 demonstrates that GMM-LRDS is able to recover both global features and local features for complex multi-modal distribution.
Execution time. To demonstrate the applicability of LRDS in practice, we provide in Figure 14 the execution clock-time of each sampling method for all multi-modal Gaussian mixture sampling settings. The conclusions are the same in the bi-modal setting: GMM-LRDS has equivalent training time to previous variational inffusion-basee approaches, with negligeable extra cost due to reference fitting.
Limitation of LRDS: high number of modes. As depicted in Figure 13 (top left), the recovery of the mode weights for a multi-modal target distribution gets worse as soon as the number of modes is large. Although GMM-LRDS performs better than its competitors, it suffers from this limitation.
I.5.3 RINGS AND CHECKERBOARD DISTRIBUTIONS
In Figure 15 and Figure I6 we illustrate the impact of the number of Gaussian mixture compo- nents when running GMM-LRDS to sample from Rings and Checkerboard distributions, detailed in Appendix by taking E {1,8, 8, 16, 32, 64). particular, we consider mixture models with diagonal covariance (first two rows) and full covariance (last two rows), fitted on MCMC samples from the target via EM algorithm. For both settings, we observe that setting J large enables to recover with more precision the support of the target distribution, leading consequently to better performance of GMM-LRDS. In the special Rings setting, where the target energy landscape is nar- row and complex, support of the distribution. On the other hand, we remark that GMM-LRDS struggles to recover the geometry of the Checkerboard distribution, while EBM-LRDS is more performant, see Figure
also notice that mixtures with full covariance provide better estimation of the
52

---

Published as a conference paper at ICLR 2025
Total Variation
1.5
W2.e
1.0
0.5
0.0
10
0.6
0.6
0.5
0.4
0.4
MMD
0.3
KS
0.3
0.2
0.2
0.1
0.0
10
20
30
40
50
60
10
20
30
40
50
60
Number modes
Number modes
SMC
LV-PIS
LV-DIS CMCD
PDDS IDEM
GMM-LRDS
RE
LV-DDS
Figure 13: Results for multi-modal Gaussian mixtures, averaged over 16 sampling runs. (Top Left): Total variation distance between weight histograms. (Top Right): Results with regularized Wasser- stein distance. (Bottom Left): Results with MMD distance. (Bottom Right): Results with KS distance. Incomplete curves refer to settings with missing results due to numerical issues.
J 8
J == 16
32
J = 64
2000
Clock time (in s)
1000
500
22
GOMN LEDN
12
Drans
IMALLCSS
22
1308
11000
GAOCLEDS
Inference
Training
MCMC sampling
Reference fitingg
Figure 14: Execution clock-time for multi-modal Gaussian mixtures with increasing number of mixture components (from left to right).
1.5.4 Φ1 FIELD SYSTEM
Table 11 shows that all competing algorithms fail at recovering the weight ratio of the of system as they are highly prone to mode collapse either on the 'negative' mode or on the 'positive' mode.
I.5.5 BAYESIAN LOGISTIC REGRESSION MODELS
Finally, we display in Table 12 the results of all considered samplers, including LRDS, when target- ing the Bayesian posterior distributions obtained from real-world data and defined in Appendix [.]] To assess the sampling quality, we compute the average predictive log-likelihood (i.e. the expected posterior distribution on the test dataset) via Monte Carlo approximation with 8192 samples. Since these distributions do not exhibit explicit multi-modal characteristics, we observe that standard MCMC-based techniques such as SMC and RE perform better than deep learning- based approaches. Interestingly, LRDS is on par or superior previous diffusion-ber -baedd approaches. For completeness, also provide in Figure 17 the execution clock-time of all samplers in this setting.
53

---

Published as a conference paper at ICLR 2025
==1
==
J 16
J 32
J = 64
Target
Ref. GMM
GMM-LRD LRDDS
Ref. GMM
GMM-LRDS
Figure 15: Results of GMM-LRDS for Rings distribution with increasingly expressive GMM for reference distribution (Top two rows) GMM density and GMM-LRDS samples for diagonal co- variance parameterization (Bottom two rows) GMM density and GMM-LRDS samples for full covariance parameterization.
Table 11: Estimation the weight ratio for the d² system with different values of h, averaged over 16 runs. The Laplace approximations, that stand for ground truth, are displayed the bottom of the table. Orange cells refer to settings with uninformative mode weight estimation (i.e., uniform mixture), red cells denote mode collapse. N/A denotes settings where results could be obtained due to numerical issues.
Algorithm
h
9 x 10-4
h= 2 10-3
2.5 x 10-3
3.5 x 10
SMC
1.28=1119
1.71+1.47
1.51+138
1.60=1.52
1.57=1.43
RE
1.00=0.12
1.01
1.04
1.11
1.07
LV-PIS
00
00
0C
00
00
LV-DDS
0.00=0.00
0.00=2.000 00
LV-DIS
1.0720.03
1.05++000
1.0544.06
1.1620.12
1.17=0.07
CMCD
1.01 20.05
0.98
1.01
0.96
0.96=0.02
iDEM
0.59
35.9+1.43
98 6 10.6
1.81 20.03
170
PDDS
N/A
N/A
N/A
N/A
N/A
Laplace thh order
1.00
1.35
1.95
2.30
3.22
Laplace 2nd order
L.00
1.34
1.90
2.23
3.08
54

---

Published as a conference paper at ICLR 2025
==1
=
16
J 32
J 64
Target
Ref. GMM
GMM- LIM-LRDDS
Ref. GMM
GMM-LRDS
Figure 16: Results of GMM-LRDS for Checkerboard distribution with increasingly expressive GMM for reference distribution (Top two rows) GMM density and GMM-LRDS samples for diagonal covariance parameterization (Bottom two rows) GMM density and GMM-LRDS samples for full covariance parameterization.
Table 12:
Average predictive log-likelihood for Bayesian posterior distributions obtained from logistic regression tasks, averaged over 16 runs. Bold font indicates best performance.
Algorithm
lonosphere 36) f
Sonar (d 62) f
German Credit 26)
Breast Cancer (d 32)
SMC
139. 3
191.
122 210.1
240.s
RE
139.5
-191.3+0.1
122 2
-
LV-PIS
156.011.
202.210.2
153.7
203.913.3
LV-DDS
-147.710.5
-195.210.7
151.2=1.3
- 128.410.9
LV-DIS
-221.3+22.6
-318.7+47.0
-206. -206.6=1.0
- -152.1+1.
CMCD
970.612.
1645.416.2
171..10.0
175.3
153. 7.0.2
4634.111299
132. 13..410.1
320.
PDDS
139.620.3
-191.3a0.1
-122 2
-102.2+0.1
GMM-LRDS
138. 0.0.55
190.5110
120.0.2..7
100.112.
EBM-LRDS
148.0+11.4
191.8
1331.205.2
-9..744.
Cancer
Credit
Sonar
Ionosphere
3000
Clock time (in s)
2000
1000
1577
Lvo0s
10mm
Inference
Training
MCMC sampling
Reference fitting
Figure 17: Execution clock-time for Bayesian posterior distributions obtained from logistic regres- sion tasks 'onosphhre', 'Sonar', 'German Credit' and BBraast Cancer' (from left to right).
55


Reviews:
[{"review_id": "FX41f9DpAq", "reviewer": "Reviewer_NXKL", "summary": "This paper introduces a new approach to sampling from multi-modal distributions called Learned Reference-based Diffusion Sampler (LRDS). The proposed method is an extension of previous variational diffusion-based samplers and tries to address their limitations related to their sensitivity to small perturbations of hyperparameters. The authors propose two approaches. The first one GMM-LRDS is based on Gaussian Mixture Models, a more lightweight approach targeting simpler target distributions. The second one EBM-LRDS uses a more computationally expensive Energy-Based Model, allowing one to sample from harder target distributions.", "questions": "I would like to suggest authors to \n1. add the comparison on real-world datasets.\n2. add the exact timing comparison of the discussed methods for training and sampling stages.\n3. expand the experiments section with an ablation study of the LRDS limitations.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The paper's structure is logical and well-organized, and the writing is clear and concise, making it easy to follow and understand.\n2. The paper demonstrates a strong mathematical foundation that supports the authors' claims. The inclusion of detailed mathematical proofs in the appendix adds to the paper's technical depth.\n3. In contrast to previous methods, which rely on hyperparameter tuning for reference distributions, the proposed approach learns these hyperparameters.", "weaknesses": "1. While the figures are generally well-presented, some may benefit from additional context or explanation in the main text. For example, the importance of parameter $h$ in Figure 4 is not obvious without checking the appendix.\n2. The authors present the experiments only on the synthetic datasets. The paper would greatly benefit from the inclusion of experiments on real-world datasets (e.g., images), which would help to demonstrate the practical applicability.\n3. The authors provide the computational complexity of the overall sampling procedure for different methods in the appendix; however, the memory complexity and exact timings might help the reader.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "NH9aARFXq3", "reviewer": "Reviewer_GBCd", "summary": "The authors propose a sampling method (named as LRDS) which extends existing works to sample from a probability density under the assumptions that they know the location of modes of the distribution while do not have access to the ground truth samples. Inspired by the observations that the existing works are sensitive to hyperparameters, the authors propose two methods (GMM-LRDS and EBM-LRDS) that surpasses the performance of the baselines on sampling from multimodal distributions.", "questions": "1. Fig. 4. Is GMM-LRDS as fast as the Laplace approximations?\n2. It is very difficult to find strengths of the proposed methods in Fig. 5. LV-DDS looks better, actually.", "limitations": "", "rating": 3, "confidence": 3, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "- Trying to solve interesting and fundamental problems.\n- Performance improvement over the baselines for estimating mode weight is shown.", "weaknesses": "- Please revise the introduction entirely from the high-level perspective. Currently, Introduction is written as if sort of a problem definition and preliminary knowledge (which are usually put in related works or method section).\n    1. What is the task to solve?     \n    2. What are the real-world applications?\n    3. What other baselines have been proposed? And what are the limitations of them?\n    4. What is the proposed method concisely and what is the improvement that is theoretically and empirically shown in the paper?\n    - It should be understandable to some extent for non-expert readers in this domain.\n    - I believe this is a basic rule of thumb in paper writing.\n- Weak motivation\n    - Why do we need to care about this paper’s assumptions (L40-L41)? Is it a realistic scenario? What are the examples in practice?\n    - Why do we need to sample from multi modal distributions? (L35-L36) Why not simply sampling from a single-mode distribution?\n- Weak experiments\n    - Higher dimensional dataset experiments are needed (e.g., MNIST or Cifar10) to show that the proposed methods can be applicable beyond the toy-level experiments.\n    - Quantitative experiments\n        - The measure for computing the distance/divergence with the target distribution is required.\n        - How is Table 2 (mode weight estimation error) connected to the actual task which is approximating the target distribution?\n        - It seems like each column of Table 2 is almost identical. Is there something meaningful to report? \n\n- Overall, the proposed method itself is interesting. However, its presentation (paper writing) needs to be significantly revised and the experiments are not good enough to show the performance improvements in sampling and to verify that the proposed method can be applicable in the real world applications.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "BDWubXx0YL", "reviewer": "Reviewer_FUSS", "summary": "This paper introduces the Learned Reference-based Diffusion Sampler (LRDS), a novel approach for sampling from multi-modal distributions. The authors identify a fundamental issue with existing variational diffusion-based methods (their high sensitivity to hyperparameters) and address this by proposing a methodology that leverages prior knowledge of target mode locations. LRDS is implemented in two variants: (1) GMM-LRDS, which is well-suited for a variety of target distributions while being relatively lightweight, and (2) EBM-LRDS, a more computationally intensive approach for harder sampling problems. Through extensive experiments on challenging multi-modal settings, the authors demonstrate that LRDS outperforms previous diffusion-based samplers.", "questions": "1. what the impact would be if random sampling was used to obtain reference samples rather than sampling that reflects a target distribution such as mcmc?\n\n2. I understand that challenge mode collapse settings are rare in general real-world datasets. Although these experiments and analyzes may not be appropriate for the contribution of this paper, I wonder whether this sampling strategy can be applied to general real-world applications.", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper demonstrates a critical limitation in existing diffusion-based sampling methods. To solve the hyperparameter sensitivity issues in previous approaches, the authors provide a well-grounded theoretical method for their solution. The mathematical framework is developed and presented clearly, making their contributions theoretically sound and academically valuable.\n\n2. The experimental evaluation is encompassing a diverse range of challenging distributions. The authors' comparative analysis against multiple baseline methods across different settings provides evidence for their method's superiority.\n\n3. The paper offers two complementary approaches (GMM-LRDS and EBM-LRDS) that provide flexibility in handling different types of multi-modal distributions. These two versions of the method demonstrate the validity of the reference-based sampling approach proposed in this paper. Figures 5 and 6 show that the proposed method generally works well for various distributions.", "weaknesses": "1. A limitation of the proposed method lies in its computational overhead, particularly in the pre-training phase of the reference process models. The authors could have provided a more detailed analysis of the computational costs and optimization strategies as the complexity and scale of the target distribution increases.\n\n2. The method's reliance on prior knowledge of mode locations (target distribution) represents a substantial practical limitation. While the authors acknowledge this requirement, the authors could address how the method might be adapted for scenarios where such information is incomplete.\n\n3. There are concerns about the method's applicability to real-world data distributions, despite its impressive performance on synthetic examples. While the paper demonstrates effectiveness in controlled settings, the approach of learning mode information from reference samples may not translate well to real-world scenarios where modes are less well-defined and data structures are more complex and noisy. The lack of validation on real-world datasets leaves uncertainty about the method's practical utility and generalizability beyond carefully constructed example distributions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "KJyAlXUQLP", "reviewer": "Reviewer_yQsx", "summary": "Unlike previous score-based diffusion sampling methods, this paper takes a step further by learning a reference process based on GMMs or energy-based models (EBMs) to facilitate sampling from multi-modal distributions. The energy-based models (EMMs) presented are novel compared to recent literature, enabling the learning of references for various challenging distributions that GMM-based references cannot handle. Additionally, the authors summarize the key differences and innovations of their method compared to other diffusion sampling techniques in Table 1.\n\nThe authors also compare their approach with other state-of-the-art algorithms, including LV-PIS, LV-DDS, and LV-DIS, on various challenging distributions. And the results are satisfactory.", "questions": "Regarding the weaknesses, my questions are as follows:\n\n1. How sensitive is your algorithm to the choice of different reference distributions (e.g., in GMMs, the variances/means of the modes; and in EBMs, the choice of hyper-parameters)?\n\n2. How does the model perform when the modes of the target distributions are widely separated or when the number of modes is higher? What about the training and sampling complexities?\n\n3. Can the authors present additional evaluation metrics to provide more straightforward assessments of sample quality?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The paper proposes using GMMs or EBMs, based on prior knowledge of the modes, as a reference for sampling multi-modal densities. The EBM framework is particularly novel, as it helps learn initial references through a neural network for sampling a range of distributions that standard GMMs reference cannot handle.\n\n2. The paper is technically sound, presenting solid theories, lemmas, and propositions in the Appendix, which provides a good supplement and well support their algorithms.\n\n3. The paper compares their algorithm with various other state-of-the-art algorithms, and the results are satisfactory.", "weaknesses": "1. While the paper highlights leveraging prior knowledge of mode locations to avoid hyper-parameter tuning challenges, I believe the sensitivity still depends on the choice of reference distributions. The authors may illustrate how different choices of references (e.g., varying the number, means, variances of GMM components, or EBM architectures) impact the sampling results.\n\n2. In the experiments, the modes of the target distributions (such as the high-dimensional Gaussian, rings, and checkerboard) are relatively close together. Have the authors tried target distributions with varying number of modes (e.g., 10, 100, 500) with increasing separation, and how does the algorithm perform in these cases?\n\n3. The authors primarily evaluate mode weight recovery error, but there are other important evaluation metrics, such as Maximum Mean Discrepancy (MMD) and Wasserstein distance. Since it seems possible to obtain true samples for each of the experiments, I believe calculating MMD and Wasserstein distance would be feasible.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)