Sampling by averaging: A multiscale approach to score estimation
Authors: Paula Cordero Encinar
Andrew B. Duncan
Sebastian Reich
O. Deniz Akyildiz
Venue: NeurIPS 2025 poster ()
Keywords: sampling from multi-modal densities
diffusion process
multiscale SDEs
Primary area: probabilistic_methods
TL;DR: We propose a training-free framework for sampling by exploiting multiscale dynamics and general diffusion paths.

Abstract: We introduce a novel framework for efficient sampling from complex, unnormalised target distributions by exploiting multiscale dynamics. Traditional score-based sampling methods either rely on learned approximations of the score function or involve computationally expensive nested Markov chain Monte Carlo (MCMC) loops. In contrast, the proposed approach leverages stochastic averaging within a slow-fast system of stochastic differential equations (SDEs) to estimate intermediate scores along a diffusion path without training or inner-loop MCMC. Two algorithms are developed under this framework: *MultALMC*, which uses multiscale annealed Langevin dynamics, and *MultCDiff*, based on multiscale controlled diffusions for the reverse-time Ornstein-Uhlenbeck process. Both overdamped and underdamped variants are considered, with theoretical guarantees of convergence to the desired diffusion path. The framework is extended to handle heavy-tailed target distributions using Student’s t-based noise models and tailored fast-process dynamics. Empirical results across synthetic and real-world benchmarks, including multimodal and high-dimensional distributions, demonstrate that the proposed methods are competitive with existing samplers in terms of accuracy and efficiency, without the need for learned models.

Sampling by averaging: A multiscale approach to score estimation
Paula Cordero-Eneinar Department of Mathematics Imperial College London
Andrew B. Duncan Department of Mathematics Imperial College London
paula.comm www.eroo-enciaar2221nnn..eeeeeee ac.uk .uk
duncan@imperian ac. uk
Sebastian Reich
O. Deniz Akyildiz Department of Mathematics Imperial College London
Institut für Mathematik Universität Potsdam sereich@uni- pot.iiii..dddmmm..m de
deniz. aimiiiaayyyddiz2@mmrraaal ac uk
Abstract
We introduce a novel framework for efficient sampling from complex, unnormalised target distributions by exploiting multiscale dynamics. Traditional score-based sampling methods either rely on learned approximations of the score function or in- volve computationally expensive nested Markov chain Monte Carlo (MCMC) loops. contrast, the proposed upproach leverages stochastic averaging within slow- fast system of stochastic differential equations (SDEs) to estimate intermediate scores along a diffusion path without training or inner-loop MCMC. Two algo- rithms are developed under this framework: MULTALMM, which uses multiscale annealed Langevin dynamics, and MULTCDIFF, based on multiscale controlled diffusions for the reverse-time underdamped variants are considered, with theoretical guarantees of convergence the desired diffusion path. The framework is extended to handle heavy-tailed target distributions using Student's 's t- t-based noise models and tailored fast process dynamics. Empirical results across synthetic and real- weal-world benchmarks, including multimodal and eigh-dimensional distributions, demonstrate that the proposed methods are competitive with existing samplers in terms of accuracy and efficiency, without the need for learned models.
-hblebbee process. Both overdamped and
Introduction
Efficiently sampling from complex unnormalised probability distributions is fundamental challenge across many scientific domains, including statistics, chemistry, computational physics, and biology. see, e.g. Gelman al. [2013] Lelievre al. [2010]. Liu (2088.. Classical Markov Chain Monte
Carlo (MCMC) methods provide asymptotically unbiased samples under mild assumptions on the
al. 2023 However, they often become computationally inefficient in practice, especially for digenssimnsinnl.. multimodal targets, due to the need for long Markov chains to ensure adequate mixing and convergence. To alleviate these issues, annealing- strategies construct sequence of smoother intermediate distributions bridging simple base distribution and the target. This principle underlies algorithms such as Parallel Tempering [Geyer and Thompson 1995. Swendsen and Wang. 1986), Annealed Importance Sampling NNal
target [Durmus and Moulines 2017 Aeausavi--osseinn
2001], and Sequential Monte Carlo [Del Moral et al. 2006
Doucet et al. 200T].
Recently, score-based diffusion methods [Hyvärinen] 2005
Song and Ermon 2020 have emerged as
powerful models capable of generating high-quality samples. These approaches simulate reverse
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

diffusion process guided by time-dependent score functions, While successful in generative modelling, they rely on access to training samples setting that differs from sampling problems, where only the unnormalised target density is available. Therefore, traditional score-matching techniques are not applicable. To address this, growing body of work has explored ways to estimate the score function using only the target density, either through neural networks using variational objectives (Chen et al. 2025 Richter and Berner 2024. Vargas et al. 2024] or training-free alternatives based on MCMC Chen et al. 2024 Grenioux et al. 2024 Huang al. 2024 Phillips et al., [2024] Saremi ct al.
2024). In this work, we present novel training-free sampling framework that eliminates the
need for the inner MCMC loop required in previous works to approximate the score function. Our
main contributions are as follows:
We demonstrate how multiscale system of SDEs can be leveraged enable efficient sampling along a diffusion path, obviating the need to estimate the score function.
In particular, we propose two algorithms: MULTALMC: Multiscale Annealed Langevin Monte Carlo (Section B.I and MULTCDIFF: Multiscale Controlled Diffusions (Section 3.2) based on discretisations of two different multiscale SDEs: annealed Langevin dynamics for general noise schedules and the reverse SDE associated with an Omstein-UUllenbeck (OU) noising process, respectively. We also explore accelerated versions using underdamped dynamics.
We provide a theoretical analysis of the proposed methods (Section 4p and demonstrate their effectiveness across different synthetic and real-world benchmarks (Section 5
2 Background and related work
2.1 Diffusion paths
The reverse process in diffusion models consists in sampling along path of probability distributions (po)teto.)). which starts at simple distribution and ends at the target distribution π x e Vv on RR Following Chehab and Korba [2024), the intermediate distributions can be expressed as follows
=(s/VA.) (</VVI-x) 2/4/2
[x(x)
(1)
(1-2.)472
where * denotes the convolution operation, v/ describes the base or noising distribution, and A, an increasing function called schedule, such that, d E [0, 1] and λ = 1. We refer to the path (Ne)tejo.1) in [T] as the diffusion path. By choosing X
(7 t) for some fixed T,, we recover the path corresponding to a forward OU noising process, a widely used approach in diffusion models [Benton Song et al. 202T]. In this case, the reverse- reverse-time dynamics can be characterised [Anderson] 1982]. requiring only to estimate the score functions V log μ4 along the path. However, for general diffusion paths, the reverse process cannot be described by dlosed-form SDE due to intractable control terms in the drift, which can be estimated using neural networks as in Allergoo et al. [2023]. Alternatively, sampling from the diffusion path can be achieved through annealed Langevin dymamics (Cordero-Eininarr et al. 2025 Song and Ermon 2019]. which also relies on score estimation but avoids dealing with intractable drift terms.
al. 2024 Chen et al., 2023
The diffusion path has demonstrated very good performance in the generative modelling literature where the score of the intermediate distributions can be estimated empirically using score matching techniques [Song and Ermon, 2019]. In the context of sampling, however, estimating the score is not available and instead we only have
more challenging, as samples from the data distribution access an unnormalised target probability distribution. particular, under the assumption that v, n are bounded and have finite second-order moments, and that V log / and VV=l are bounded, the expression of the score of the intermediate distributions is given by
V Vlogue(>)) VI d
E Bp),z())
V Vlogy
VI (
x' Y
)
vxe E
(v)
VV=
) Vx
Y
(2)
V=(w/VAe) and we have used that V * g) (Vf))g=ff (Vg) when f and y are differentiable functions. Notably, in the limits X, = and X, = 1, the conditional distribution Pt.z converges to Dirac delta function. The score in ② involves an expectation over Pt. which typically requires MCMC sampling and results in a computationally expensive nested-loop structure. We avoid this by exploiting multiscale dynamics, as described in the following section.
where (v) XV

---

2.2 Multiscale dynamics and stochastic averaging
Consider the slow-fast system of SDEs
dxx= Xi, V.)dt + v2dB. dYi ft, X1, Yijct / /dd.
where 0 < < 1 controls the scale separation between the slow X, and the fast component Y, We assume that, for = and fixed t and X,, Y, is ergodic. As E + 0, the fast component Y, evolves on much shorter timescale than X, due to the time-rescaling property of Brownian motion. Intuitively, this means that Y, rapidly reaches its stationary distribution while X, remains substantially unchanged. In the limit - the slow-process X₁ converges the averaged dynamics X, which is given by
dx. = b(t, X.)dt Vad.
b(t,z)) = By- [b xr V)
where Pr.z is the invariant measure of the frozen process dY, f(l Y.)ds + v2dB.. In computa- tional statistics, stochastic slow-fast systems have been leveraged simulate the averaged dynamics efficiently [Akyildiz et al. [2024] Harvey et al. 2011 Pavliotis and Stuart [008 WWeinan et al. 2005]] In our context of sampling from diffusion paths, the averaged drift 6(t, x) corresponds the score v log pe (z), which is defined as an expectation (see Eq. (2. . By leveraging multiscale dynamics, we can approximate this expectation without relying on a nested-loop sampling structure, key insight
that underpins our method, as summarised in Section
2.3 Related work
Diffusion-based samplers Several recent works have proposed sampling algorithms based on the
diffusion path. We divide them into two categories: neural samplers and learning-free samplers.
Neural samplers such as those in Chen et al. [2025), Noble et al. 22025, Richter and Berner (2024). Vargas al. (2024), Zhang and Chen [022] estimate the time-dependent drift function of the diffusion process using a parametrised model, typically a neural network, by solving a variational inference problem defined over the space of path measures. As in diffusion generative models, their performance is limited by the expressiveness of the chosen model class. In contrast, learning-free
samplers [Chen et al., 2024 Grenioux al. 2024 He al.) 2024b Huang et al. 2024 Phillips et al.) 2024] Saremi et al. [2024] do not require any traming similar to our work. Instead, they estimate the score function log #4 using MCMC samples from the conditional distribution Pt.z. That is, these approaches involve a bi-level structure, with an inner MCMC loop used at each step of the outer diffusion- -based sampler. this sense, they implement a Langevin-within siieii-wttti--LLngeeee strategy. Our approach avoids nested loops by leveraging multiscale dynamics along the diffusion path. This results in a Langevin-by-LLngeein sampler that is both simpler and more computationally efficient.
Besides, while acceleration techniques based on underdamped dynamics [Duncan et al., 2017 Monmarché 2020] have shown promise in the generative modelling literature [Dockhorn al. 2022a b. Holderrieth et al., 2024], and recent work by [llsssing et al. (2025) introduces neural underdamped diffusion samplers, no training-free accelerated samplers based on diffusion processes have been proposed. Addressing this gap is one of the key contributions of our work.
Our approach is also closely related to the class of proximal samplers CCen al. 2022 Lee al. 20221b. These methods define an auxiliary joint distribution of the form T [(, 2) ग. (xx) exp( - (22(22)) and perform Gibbs sampling by alternating between the convolutional distribution (=z[²) and the denoising posterior (=]=)) In practice, sampling from the denoising posterior #xx22) involves finding a mode using gradient- optimisation techniques, followed by rejection sampling around that mode. This step can be computationally expensive and inefficient, especially in dinhe-iimnsional settings. Our method improves upon these approaches by posterior sampling step with a fast-timescale dynamics, enabling
Proximal samplers
replacing the proximal-point-spy more efficient exploration of the denoising posterior. Furthermore, unlike proximal samplers that operate at a fixed noise level, our approach targets the entire diffusion path, smoothly interpolating between the base distribution and the target.
3

---

3 Sampling using multiscale dynamics
A key challenge in sampling from the diffusion path in T is accurately estimating the score function when we do not have direct access to samples. Since the score of the intermediate distributions along the diffusion path is computed an expectation over the conditional distribution Pt,.: (see Eq. (2)) we propose to sample from the time-dependent conditional distribution Pt.æ using an alternative diffusion process with shorter timescale than the original dynamics following the path in (T) This enables both processes to run in parallel, unlike existing methods that execute them sequentially. In those approaches, each iteration of the original diffusion process (outer loop) requires multiple time steps of the process targeting the conditional distribution Pt.z: (inner loop) to learn the score function, leading to increased computational cost and higher number of evaluations of the target potential. By employing different timescales for each diffusion process, our method improves sampling efficiency while ensuring convergence the correct target distribution,
justified by stochastic averaging theory (Liu et al. 2020]. We explore this approach in two settings: (i) when the original diffusion is a Langevin dynamics driven by the scores of the intermediate probability distribution V log frr (Section 3.1. and (ii) when we use the reverse SDE corresponding to an OU noising process (Section 3.2)
3.1 MULTALMC: Multiscale Annealed Langevin Monte Carlo
Following Cordero- Cordero-Encinar et al. 2022 Song and Ermon 2019). we use time inhomogeneous Langevin SDE to sample from the diltusion path 1 given by
dx,- Vlog (x, VadB, t [0, 1/m],
(3)
where je denotes the reparametrised probability distributions from the diffusion path (ū, (xz))eju.1/n) for some 0 < < 1, Xo no
and (Be)r20 is Brownian motion. We consider both Gaussian diffusion paths where the base distribution v is a Normal distribution and heavy-tailed diffusions corresponding to a Student's base distribution. Since the scores V log Pu are intractable and computed expectations over Pt.. (see Eq. (2). we adopt multiscale system to approximate the averaged dynamics B). The fast process targetting the conditional distribu- tion Pr.y can follow any fast-mixing dynamics. In particular, we explore the use of overdamped Langevin dynamics or modified Itō diffusion (He et al. 2024a], but it is important to remark that our framework is more general. When using overdamped Langevin dynamics, the fast process is given by
dY; 1 - log Pr.x. (YG)dt +
2 Lan. iE 1/k],
Ye Vent
1
Yi X1
v Vlog Pr.x (v) =
VV= Vhnt
v Vlogv
V1
Aus
V1 Axet
where pt.x, is a reparametrised version of Pt.x, and (Bz) Brownian motion on Rd suggests sampling from the following stochastic slow-fast system
This
V logu
Xe-X dt + v2dB, if Aut < (V=Aw
dX.
vT
Ane
1
VV
Yo
di + Vdd.
if Aut 2
(4)
Vinc
dY;
VV
VI-ANe Vogy (iT
dt + / /ddB.
VAni
VAni
This system converges the averaged dynamics described in as 0, since taking the expectation of the drift term in the X, dynamics with respect to Y, ~ Pr.x, recovers the score function V log ft, (see Equation (2) and Section 4). For numerical implementation, we use an adaptive scheme that alternates between the two expressions for X, depending on the value of Aet in order to avoid numerical instabilities as ARt approaches 0 or 1. The resulting system of SDEs can be discretised using different numerical schemes, however, standard integrators such as Euler- Maruyama (EM) will be unstable when is small, due the large-scale separation between the slow and fast processes. To mitigate this, we leverage discretisations which remain stable, despite the stiffness of the SDE, such as the exponential integrator [Hochbruck and Ostermann 2010 and the SROCK method (Abdulle et al. 2018). Furthermore, since Pt.z. converges to Dirac delta distribution centred at when Aut and centred at & when ARe = 1, we apply further slight modification to the slow-fast system to improve numerical stability (see App. B.I.I for details).

---

Accelerated MULTALMC When implemented in practice, the overdamped system in requires very small value of (which corresponds to slowly varying dynamics driven by V log (uy to accu- rately follow the marginals of the path (pe))e(0.1/ and ultimately sample from the target distribution. This results in large number of discretisation steps, thus making the method computationally expen- sive, see App. E.] for more details. In contrast, underdamped dynamies exhibit faster mixing times [Bou-Rabee and Eberle 2023, Eberle and Lörler [024], which intuitively helps the system explore the space more efficiently and better track the intermediate distributions along the path (pe))c(0.1// These dynamics have also shown strong empirical performance in the neural sampling literature [Blessing et al., 2025]. requiring fewer discretisation steps. Motivated by these advantages, we propose augmenting the state space with auxiliary velocity variables V,, and adopting the following underdamped Langevin diffusion to sample from the path (00))c(011/4)
f [xx, Vidt dv; V log (x))dt - -.M-VVdd + /21dB,
for t E [0, 1/v]..
(5)
The mass parameter M determines the coupling between position X, and velocity Vi, while the friction coefficient I controls the strength noise injected into the velocity component, Both M and I are symmetric positive definite matrices, in practice we consider M = m.fd and T = TTd. The behaviour of the system is governed by the interplay between M and I' [McCall [2010). Based on this, we propose to use the following underdamped slow-fast system sample from the target distribution
dxx; M Vidtt
{
x V
Vloge
- IM- di + vrrdB if Aut < TM VV) dt + v2TdB,
dV. =
(6)
VV
if ARt 2
VI-Amt log (T- ) dt + V2dB.
dY;
VV
Vdst
Note that in the underdamped setting, both X,, V1 are treated as slow variables, and as s + 0. this system converges to the averaged dynamics in (5) To implement the system in practice, we combine a symmetric splitting scheme (OBABO) [Mnnaachhe] 2020 for the slow-dynamics, and discretise the fast-dynamics using SROCK method [Abdulle et al.] 2018 We will evaluate the performance of the proposed sampler based on this underdamped slow-fast system (Al (Algorithm 33) referred to as MULTALMC, under two choices of the reference distribution standard Gaussian and Student's t distribution. A simplified version of the implementation is presented in Alg. D while a more detailed description is provided in App. B.1.2
Heavy-tailed diffusions When the target distribution has heavy tails, standard Gaussian diffusions often fail to capture the correct tail behaviour (Pandey et al. [2025] Shariatian al.) 2025 This motivates the use of heavy-tailed noising processes, which have been shown to offer theoretical guarantees in such settings [Cordero-Encinar et al. 2025].
We propose sampling along heavy-tailed diffusion path (#e))ejo,1//)) where the reference distribu- tion v is chosen Student's t-distribution with tail index a, i.e., v oc (1 (1+1) (z(12/a)-(a+d)/// This can be implemented using the underdamped slow-fast system defined in (6) However, when the target distribution is heavy-tailed, the use of overdamped Langevin dynamics for the fast component results in slow convergence Moossavi-Hosseiii et al., 22033, Wang 2006). This motivates the consideration of alternative diffusion processes for the fast dynamics that mix more efficiently under such conditions. Specifically, we explore a modified Itô diffusion proposed in He et al. [2024a). defined as
dY. = - + d- VU,.., (v)dt +
/2 AU dB
Pr.x.()) (3) (v++ || - 211²
(a++))
U,,,x (2) (a+d))
a(1 Ant)
This modified Itô diffusion for heave-tailed distributions is shown to have faster convergence guaran- tees when the target has finite variance (He et al.) 2024a]. By replacing the fast process in (6) with the diffusion defined in 77. we obtain samplers better suited for heavy-tailed targets, with improved convergence properties.
3.2 MULFCDIFF: Multiscale Controlled Diffusions
Annealed Langevin dynamics (3) offers convenient approach to sampling from the target distribution for general schedules X; that interpolate between the base and target distributions in finite time.
5

---

Algorithm 1 MULTALMC sampler: accelerated version
Require:
Schedule function As, value for As and X, friction coefficient 1, mass parameter M, number of Tr. 1/1s, step size Ti+1 Ti. Constants
sampling steps L. time discretisation 0 = To
J., M; & for SROCK step from 14). 15
i0 Aut As (1o
Initial samples Xo (),, ((O, ~V(0,)). Xat
Ant if S Act 1.
for - 0 to do
< , (()) (0, I)
ifo < ArT, < 1 As then
stee- for velocity component
1/
hi
V logu
Xi-Mi
if MeT, Ā
V/' 1 ( Vi hiTe) 2
VI"
3/1 Vi1-Nam
Vi' 2//ien,
VV=
if Ne 2 2
V
Full EM-step for position component X1+1 X, + 1u-M-10
2h (2)
Full step for fast component
V441 f(V,, Xi+1, Nen E, +
Half Htep velocity component
V/"
by
V
Vloge
X411-141 X
if NeT, < <x,
V41 (1 (1- -rM V/"" - "0"+vhil VhiTg
Vi
Vi"
by VV=
Y1+1 Mer
ir ir 2 x,
if As S Art, S then
Half step velocity component
V/' - 1 rM-
Vi + V
1/" Vi¹ byg Vv. (Xi)
Full EM-step for position component
X1+1 X1 hy 'V"
step for velocity component
Vith = Vi"
hy VV, (X4+1), V441 =
(1 11-M)M-) vi" + V
2
end for
return (Xz, VT.)
However, a key drawback of annealed Langevin dynamics is that it introduces bias, even perfectly simulated (see Section 4.2). This bias can be corrected by incorporating control terms. Although such control terms are generally intractable for arbitrary schedules, they can be computed explicitly in the case of an OU noising process. The multiscale approach introduced in the previous section naturally extends to this controlled setting, in both overdamped and underdamped regimes. Given the improved efficiency of underdamped dynamics, we focus our discussion on the underdamped formulation (see App. B.2. for details on the overdamped case). When using an OU process as the forward noising model, the underdamped diffusion, initialised at X and V 0 NN0, I), is defined
dx. M dt
dv. = (-x---x---a-- di + /2r dB(.
In this case, to ensure efficient and stable dynamics, we adopt the critical damping regime, where r 4M. This choice leads to fast convergence without oscillations [McCall 2010]. The corresponding time-reversed SDE is
Tax. -pp - 2 dt av. (xx - 2 log o V.)) di + v2T dB.,
(8)
where de is the marginal distribution of the forward process which has the following expression
q: (z. v) = ((v)p(vo))((z e/y, vo) dy deo.

---

with & = NN I) and conditional distribution Pz(x, vly, vo) normally distributed with mean mily vb) and covariance Σ, given by
mely, 0o) = 8
y Ub
E.
1/- Alt
(
21 c [) d A= A-
4 -4r-1
We observe that vy can analytically integrated out the expression for pe(a, uly, vo). Moreover, the conditional distribution pe ((la,,) y) remains Gaussian, with mean m2 y) and covariance of I, where ( y) is linear y. This structure allows us to express the score V log qi (x, v) as
V. log (e v) = V. qc(l)()))) v pr. (x, Y)
v S. (x. Eyla.((M
(9)
of
where fa a linear function of its arguments, see App. the expressions above, we define multiscale SDE system that enables sampling from the target distribution via the reverse SDE (8. without requiring prior estimation of the denoiser Eyja,"lY that appears in the expression for v log ( v) (9 In this formulation, the fast, process is modelled by an overdamped Langevin SDE that targets the conditional distribution Y (X,, V.), given by
B.2 for a detailed derivation. Using
dY; 1 =Yi) Vy, log V.lxa) P.B..
Combining this fast process with the reverse SDE, we obtain the following stochastic slow-fast system, which we use for sampling from the target distribution
dx. -
dV. 1
2 (V.--v- (x)) di /2r dB, 02-
(10)
dY, = (Vlog= d + Vy, log pr (Xx, V V.lx)) dt 2dd..
To implement this sampler, referred to as MULTCDIFF,
construct novel integrator inspired by Dockhorn et al. [2022b]. Specifically, we leverage the symmetric Trotter splitting and the Baker-Campiln Campbell-Hausdortt Hausdorff formula (Strang] 1968 Trotter, [959 Tuckerman 2010) to design stable and efficient symmetric splitting scheme, see App. B.2 for details. We also outline in the appendix the difficulties of extending controlled diffusions to the heavy-tailed scenario.
4 Theoretical guarantees
In this section, we identify the different sources of error the proposed sampling algorithms. For MUL.ALLMC, the total error consists of three components: the discretisation error from numerically solving the slow-fast system, the convergence error of the slow-fast system to its averaged dynamics, and the bias arising from the time-inhomogeneous averaged system, whose marginals differ from those of the diffusion path in [ this this last component discussed in 4.2 For MULTCDIFF, the error bound includes initialisation error present in the error bounds of diffusion models [Benton] et al. 2024| Chen et al. 2023), the discretisation error of the slow-fast system, and the convergence of the multiscale system the averaged dynamics. The discretisation error depends the choice of numerical integrator, underscoring the importance of an accurate and stable integrator scheme. We do not analyse this error in detail as it follows established methods (Leimkuhler et al.) 2024 Monmarché 2020). We now study the convergence properties of the slow component X₂ to the solution of the corresponding averaged dynamics that evolves along the time-inhomogeneous diffusion path in (T).
4.1 Convergence of the slow-fast system to the averaged dynamics
Building on results from stochastic averaging theory [Liu et al., 2020 , our goal is to derive sufficient conditions on the coefficients of the multiscale systems namely, y, 4 (6) and (10) that guarantee strong convergence of the slow component to the averaged process as e - O. We focus on the case where the base distribution is Gaussian where v a heavy-tailes distribution, such as a Student's distribution, presents significant technical challenges, which fall outside of the scope of this work. We outline these challenges in App. To establish convergence guarantees, we assume the following regularity conditions.
~N00,) 7). Extending the convergence analysis to the case
A1. The targel distribution n has density with respect to the Lebesgue measure, which we write т X e V=, The potential Vz has Lipschitz continuous gradients, with Lipschitz constant I In addition, Vn is strongly convex with convexity parameter Ma > 0.
7

---

A2. The schedule X : R + [0, 10 is non-decreasing function of weakly differentiable and Hölder contimuous with exponent 7/1 in (0, 1].
Under these assumptions, obtain the following result.
Theorem 4.1. Let the base distribution v ~ NN. D). Suppose the target distribution x and the schedule function x satisfy assumptions A and A2 respectively. Then, for any E (0. where Eo specified in App. C.I and any given initial conditions, there exists unique solutions
{(X;, ) 1 0}. [(x, Vi, Y.), 2 0} and {(x,, V ), t. 2 03 the slow-fast stochastic systems
and T0P. respectively. Furthermore, for any > 0, it holds that
lim E c-
XxIl") 0, || re(o,T]
sup
where X denotes the solution of the averaged system and X is the corresponding slow component of the multiscale systems 4
6 and (TO).
The proof is provided in App.
The theorem applies to both the annealed Langevin dynamics and the controlled diffusion. We note that the strong convexity condition in A which guarantees exponential convergence the fast process to its unique invariant measure when & = 1 is fixed restrictive assumption. However, shown in Section [5] our proposed algorithms show strong empirical performance on benchmark distributions that do not satisfy this condition. Extending the theoretical analysis cover non-log-concave targets remains an important future direction.
4.2 Bias of the annealed Langevin dynamics
The bias of the overdamped annealed Langevin dynamics 3 has been studied in prior work (Cordero Encinar et al. [2025] Guo et al. 20225. Here, we focus on quantifying the bias introduced by the underdamped averaged system (5) relative to the true diffusion path in the augmented state space
with the velocity. At fixed time 7, the invariant distribution of the underdamped averaged system takes the form
(a v) exp (- 1 - (og
(11)
Note that the Hamiltonian is separable, meaning that and v are independent. It is important to emphasise that, even when simulated exactly, diffusion annealed Langevin dynamics introduces bias, as the marginal distributions of the solution of (5) do not exactly match the prescribed distributions (pr):- A key quantity for characterising this bias will be the action of the curve of probability measures J (14 kejo. defined in ( interpolating between the base distribution v and the target distribution π, denoted by A(u). As noted by Guo et al. 2025 the action serves as measure of the cost of transporting v to along the specitied path. Formally, for an absolutely continuous curve of probability measures [isinn 2007 with finite second-order moment, the action is given by
A(p) I-" [im | mi) de.
Based on Theorem from et al. (2025), we have the following result.
Theorem 4.2. Let QU-ALD
(...-ALL))E(0o1/// be the path measure of the diffusion annealed Langevin dynamics 5). and Q = (P:)tejo,1/vl that of a reference SDE such that the marginals at each time have distribution Pr [1). go.U-ALD = Po. the KL divergence between the path measures is given by
KL (Q [[QUALD) m((u)/4
The action A(μ) is bounded when the target π has bounded second order moment and under mild conditions on the schedule function, see Cordero--Encinee et al. 2025 Lemmas 3.3 and 4.2]. Full details and complete proof are provided in App. C.3
5 Numerical experiments
We now evaluate the performance of both proposed underdamped samplers, based on the annealed
Langevin dynamics (Section 3.J and the controlled diffusions formulations (Section 3.2)
across
range of benchmark distributions. Full details of cach benchmark are provided in App. D
8

---

Mixture of Gaussians (MoG) and mixture of Student'st t (MoS) distributions.
Rings: Two-dimensional distribution supported on a circular manifold.
Funnel: 10-dimensional distribution defined by x(x) x MM n²) Na;; for (z:)) RR with 1j = [Neal 2003).
Double well potential (DW): dd-dimensional distribution given by exp( ( 5)2
Lizm+1 x3) with m E N and E (0, oo). Ground truth samples are obtained using rejection sampling with Gaussian mixture proposal distribution [Midgley et al.. 2023).
Examples from Bayesian statistics: Posterior distributions arising from Bayesian logistic regres- sion tasks on the Ionosphere (dimension 35) and Sonar (dimension 61) datasets.
Statistical physics model: Sampling metastable states of the stochastic Allen-Cahn equation o² model (dimension 100) [Albergo et al. 2019 Gabrié et al., [0022]. At the chosen temperature, the distribution has two well distinct modes with relative weight controlled by parameter h. Following Grenioux et al. 2024], we estimate the relative weight between the two modes.
While we focus here on the underdamped versions of our algorithms. we have also evaluated their overdamped counterparts. These approximately require an order of magnitude more steps to achieve comparable performance. Detailed results and comparisons are presented App. E. Besides, we have set the base distribution, denoted of the mixture of Student's t benchmark distributions, we compare the performance of using either a standard Gaussian or a Student's s base distribution. For the latter, we further examine the impact of using a modified Itô diffusion for the fast dynamics 7 as opposed to a standard Langevin diffusion.
in our algorithms, to be a standard Gaussian. In the case
We compare our approach against a representative selection of related sampling methods: Sequential Monte Carlo (SMC) [Del Moral al., 2006 Doucet et al.] 2001], Annealed Importance Sampling
(AIS) [Neal] 2001]. Underdamped Langevin Monte Carlo (ULMC) [Cheng et al.) 2018 Neal et al.
2011 Parallel Tempering (PT) [Geyer and Thompson 1995 Swendsen and Wang 1986 Diffusive
Gibbs Sampling (DiGS) [Chen et al., [2024], Reverse Diffusion Monte Carlo (RDMC) (Huang
al.
2024] and Stochastic Localisation via lterative Posterior Sampling (SLIPS) [Grenioxx et al.
2024 We ensure a fair comparison by using the same number of energy evaluations and the same initialisation based on a standard Gaussian distribution for all baselines. For completeness, we also report the performance of SLIPS with optimal initialisation (App.D)) DD)-hichh assumes knowledge of the distribution's 's scalar variance, unavailable in practice. For the di² model, we only compare results against Laplace approximation and SLIPS algorithm, as the other methods provide degenerate samples. Since our approach is learning-free free, we exclude comparisons with ncural-based samplers as it would be difficult to make comparisons at equal computational budget.
Our methods consistently achieve strong performance across all benchmarks, outperforning all baselines. They also exhibit low variance in results (see Tables and Among our samplers, the one based on controlled diffusions, MULTCDIFF, generally achieves slightly better results than the annealed Langevin-based sampler, MULTALMC. Additionally, Figure [1a illustrates the estimation of the relative weight between the two modes for the d model, where our algorithms demonstrate
Table 1: Metrics for different benchmarks averaged across 30 seeds. The metric for the mixture of Gaussian (MoG) and Rings is the entropy regularised Wasserstein-2 distance (with regularisation parameter 0.05) and the metric for the Funnel is the sliced Kolmogoroy-Smimov distance.
8-MoG (L)
40-MoG (L)
40-MoG (L)
Rings (L) (d [d=3) 2)
Funnel (L) (d 10)
Algorithm
(d= 2)
(d= 2)
(d = 50)
SMC
5.26 + 0.19
4.79 + 0.16
52.17 + 1.32
0.29 = 0.05
0.034 + 0.005
AIS
5.53 = 0.23
5.01 0.20
65.83 + 1.66
0.20 = 0.03
0.040 + 0.006
ULMC
6.90 = 0.28
0.72
109.26 + 2.90
0.35 : 0.06
0.123 + 0.011
PT
3.91 + 0.10
0.92
2.551
0.20 0.04
0.033 4 0.004
DIGS
1.22 + 0.09
0.95 0.03
21.87 0.97
0.037 + 0.005
RDMC
1.01 0.21
0.98 0.07
36.10
0.33 0.01
0.080 + 0.007
SLIPS
1.15 1.15
1.04
23.71 0.65
0.26 0.01
0.039 + 0.006
MULTALMC
0.65 + 0.07
0.91 + 0.04
20.58 + 0.71
0.18 + 0.02
0.032 + 0.005
MULTCDIFF
0.62 + 0.10
0.93 0.06
20.13 + 0.59
0.19 + 0.03
0.081 + 0.005

---

Table 2: Performance metrics for different sampling benchmarks averaged across 30 seeds. The metric for the double well potential (DW) is the entropy regularised Wasserstein distance (with regularisation parameter 0.05) and the metric for the Bayesian logistic regression on lonosphere and Sonar datasets is the average predictive posterior log-likelihood on test dataset.
5-dim DW ()) (m = 5,6 4)
10-dim DW (L)
50-dim DW (L)
lonosphere (1)
Sonar (+) (d = 61)
Algorithm
(m 5,8 3)
(m = 5, 2)
(d = 35)
SMC
4.06 1 0.13
11.86 + 0.70
28.92 0.99
-87.74 0.10
-111.00 0.11
AIS
4.11 0.17
12.30 + 0.61
39.20 + 1.07
-88.11 = 0.13
-111.15 0.17
ULMC
7.87 0.51
25.21 1.00 1.00
62.74 + 1.52
-116.37 1.85
-173.61 2.48
PT
+.20
4.97
16.13 + 0.89
-87.91
01.1.
DiGS
2.51 + 0.18
5.02 +
17.04 0.95
-88.02
-110.53
RDMC
3.68 + 0.32
9.57 + 0.75
25.21 + 0.73
-108.44
1
SLIPS
2.46 + 0.21
5.09 0.39
18.15 + 0.50
87.34
-110.144
MULTALMC
2.08 + 0.16
4.48 +: 0.40
14.03 + 0.63
86.85 - 0.09
-009.05 0.13
MULTCDIFF
1.95 + 0.24
4.23 + 0.37
13.98 0.56
86.33 + 0.10
-09.60
PUICALMC:
Heavy-tazled with
-X
Lap lace upproxination
approximation (2nd order)
PUTCALMC: Hoavy-tailee
Mul LAUNG
Pullctiiff:
Saussian
Mul ECOLFF
SLIPS
3-
2-
1-
Dimension
(a) Estimated mode weight ratio of ot different
(b) Regularised Wa for MoS in different dimensions.
Figure I: Results for sampling benchmarks.
competitive performance. For the mixture of Student's t-distributions, Figure Ib shows that using the modified Itô diffusion for the fast dynamics, combined with a Student's base distribution (red boxplots) outperforms both our other proposed algorithms and the baselines. Further experimental details and additional results are provided in Appendices D and
Diseussion
In this work, introduce framework based multiscale diffusions for sampling from an normalised target density. In particular, we propose two samplers, MULTALMC and MULrCDIFF depending the dynamics used for the slow process: annealed Langevin dynamics controlled diffusions, respectively.
establish theoretical guarantees for the convergence of these sampling algorithms and illustrate their performance
range of high- nigh-dimensional benchmark distributions.
Our approach has certain limitations. Notably, the theoretical guarantees rely on stringent assumptions, which we aim to relax in future work. Additionally, the current method requires manual tuning of hyperparameters such as step size 8, scale separation parameters & and friction coefficient r. Automating this tuning process remains an important direction for future research. Further research could explore extending the controlled diffusion framework to heavy-tailed target distributions, which pose additional challenges developing more efficient numerical schemes for implementing the proposed multiscale samplers.
10

---

Acknowledgements
PCE gratefully acknowledges support from the EPSRC through the Centre for Doctoral Training in Modern Statistics and Statistical Machine Learning (StatML),, grant no. EP/S023151/1. SR's work was partially funded by the Deutsche Forschungsgemminschaft (DFG) under Project- 318763901 SFB1294.
References
Assyr Abdulle, Ibrahim Almuslimani, and Gilles Vilmart. Optimal Explicit Stabilized Integrator of Weak Order for Stiff and Ergodic Stochastic Differential Equations. SIAM/ASA Journal on Uncertainty Quantification, 6(2):937-964. 2018.
O Deniz Akyildiz, Michela Ottobre, and lain Souttar. A multiscale perspective maximum marginal likelihood estimation. arXiv preprint SEGi:2406.0418877 2024.
M. Albergo, G. Kanwar, and P. E. Shanahan. Flow-based generative models for Markov chain Monte Carlo in lattice field theory. Phys. Rev. D, 100:034515, 15, Aug 2019.
Michael Albergo, Nicholas Boffi, and Eric Vanden-Eijnden. Stochastic interpolants: unifying framework for flows and diffusions. arXiv preprint urXiv:2303.08797, 2023.
Luigi Ambrosio and Bernd Kirchheim. Rectifiable sets in metric and Banach spaces. Mathematische Annalen, 318(3):577-555, 2000.
Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Grudient Flows: In Metric Spaces and in the Space of Probability Measures. Lectures in Mathematics. ETH Zürich. Birkhäuser Basel, 2008
Brian DO Anderson. Reverse-time diffusion equation models. Stochastic Processes and their Applications, 12(3):313-326, 1982.
Joe Benton, Valentin De Bortoli, Arnaud Doucet, and George Deligiannidis. Nearly d-linear con- vergence bounds for diffusion models via stochastic localization. In The Twelfth International Conference on Learning Representations, 2024.
Richard Bertram and Jonathan E. Rubin. Multi-timescale systems and fast-slow analysis. Mathemati- cal Biosciences, 287: 2017.
Denis Blessing, Xiaogang Jia, Johannes Esslinger, Francisco Vargas, and Gerhard Neumann. Beyond ELBOs: A large-scale evaluation of variational methods for sampling. In Forty-first International Conference on Machine Learning. 2024.
Denis Blessing. Julius Berner, Lorenz Richter, and Gerhard Neumann. Underdamped Diffusion Bridges with Applications to Sampling. In The Thirleenth International Conference on Learning Representations, 2025.
Nawaf Bou-Rabee and Andreas Eberle. Mixing time guarantees for unadjusted Hamiltonian Monte Carlo. Bernoulli, 29(1):75 104, 2023.
Omar Chehab and Anna Korba. A practical diffusion path for sampling. In ICMI. 2024 Workshop on Structured Probabilistic Inference & Generative Modeling, 2024.
Junhua Chen, Lorenz Richter, Julius Berner, Denis Blessing, Gerhard Neumann, and Anima Anand- kumar. Sequential Controlled Langevin Diffusions. In The Thirteenth International Conference on Learning Representations, 2025.
Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru Zhang. Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. In The Eleventh International Conference Learning Representations, 2023.
Wenlin Chen, Mingtian Zhang, Brooks Paige, José Miguel Hemmández-Lobato, and David Barber. Diffusive Gibbs sampling. In Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp. editors, Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 7731-7747. PMLR, 21-27 Jul 2024.

---

Yongxin Chen, Sinho Chewi, Adil Salim, and Andre Wibisono. Improved analysis for a proximal algorithm for sampling. In Po-Ling Loh and Maxim Raginsky, editors, Proceedings of Thirty Fifth Conference on Learning Theory, volume 178 of Proceedings of Machine Learning Research, pages 2984-3014, PMLR, 02-05 Jul 2022.
Xiang Cheng, Niladri S. Chatterji, Peter L.. Bartlett, and Michael I. Jordan. Underdamped L Langevin MCMC: A non-asymptotic analysis. In Sébastien Bubeck, Vianney Perchet, and Philippe Rigollet, editors, Proceedings of the 31st Conference On Learning Theory, volume 75 of Proceedings of Machine Learning Research, pages 300-323, PMLR, 06-09 Jul 2018.
Sinho Chewi, Chen Lu, Kwangjun Ahn, Xiang Cheng, Thibaut Le Gouic, and Philippe Rigollet. Optimal dimension dependence of the Metropolis-Adjusee Langevin Algorithm. In Mikhail Belkin and Samory Kpotufe, editors, Proceedings of Thirty Fourth Conference on Learning Theory, volume 134 of Proceedings of Machine Learning Research, pages 1260-1300. PMLR, 15-19 Aug 2021.
Paula Cordero-Encinat O. Deniz Akyildiz, and Andrew B. Duncan. Non-asymptotic Analy- sis of Diffusion Annealed Langevin Monte Carlo for Generative Modelling. arXiv preprint arXiv:2502.09336, 2025,
Marco Cuturi, Lactitia Panagaapxxan,hos, Yingtao Tian, Charlotte Bunne, Geoff Davis, and Olivier Teboul. Optimal Transport Tools (OTT): A JAX Toolbox for all things Wasserstein. arXiv preprint arXiv:2201.12324, 2022.
Pierre Del Moral, Arnaud Doucet, and Ajay Jasra. Sequential Monte Carlo samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3):41 1-166, 2006.
Tim Dockhorn, Arash Vahdat, and Karsten Kreis. GENIE: Higher-Order Order Denoising Diffusion Solvers. In Advances in Neural Information Processing Systems, 2022a.
Tim Dockhorn, Arash Vahdat, and Karsten Kreis. Score- -osre-Based Generative Modeling with Critically- Damped Langevin Diffusion. In International Conference on Learning Representations (ICLR), 2022b.
Arnaud Doucet, Nando de Freitas, and Neil Gordon. Sequential Monte Carlo Methods in Practice. Springer New York, 2001.
A. B. Duncan, N. Nüsken, and G. A. Pavliotis. Using Perturbed Underdamped Langevin Dynamics to Efficiently Sample from Probability Distributions. Journal of Statistical Physics, 169(6): 1098- 1131, 2017.
Alain Durmus and Éric Moulines. Nonasymptotic convergence analysis for the unadjusted Langevin Algorithm. The Annals of Applied Probability, 27(3):151 2017.
Andreas Eberle and Francis Lörler. Non-reversible lifts of reversible diffusion processes and relaxation times. Probability Theory and Related Fields, 2024.
Marylou Gabrié, Grant M. Rotskoff, and Eric Vanden Eijnden. Adaptive Monte Carlo augmented with normalizing flows. Proceedings of the National Academy of Sciences, 119(10):2109422011 19, 2022.
Andrew Gelman, John B. Carlin, Hal S. Stern. David B. Dunson, Aki Vehtari, and Donald B. Rubin. Bayesian Data Analysis. Chapman and Hall/CRC., 3rd edition, 2013.
Charles J. Geyer and Elizabeth A. Thompson. Annealing Markov Chain Monte Carlo with Applica- tions to Ancestral Inference. Journal of the American Statistical Association, 90(431):909-925 920, 1995.
Louis Grenioux, Alain Oliviero Durmus, Eric Moulines, and Marylou Gabrié. On sampling with approximate transport maps. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 1 11698-188 11733. PMLR, 23-29 Jul 2023.
12

---

Louis Grenioux, Maxence Noble, Marylou Gabrié, and Alain Oliviero Durmus. Stochastic localization via iterative posterior sampling. In Forty-first International Conference on Machine Learning. 2024.
Wei Guo, Molei Tao, and Yongxin Chen. Provable Benefit of Annealed Langevin Monte Carlo for Non- log-concave Sampling. In The Thirteenth International Conference on Learning Representations, 2025.
Emily Harvey, Vivien Kirk, Martin Wechselberger, and James Sneyd. Multiple Timescales, Mixed Mode Oscillations and Canards in Models of Intracellular Calcium Dynamics. Journal of Nonlinear Science, 21(5):639 683, 2011.
Ye He, Tyler Farghly, Krishnakumar Balasubramanian, and Murat A. Erdogdu. Mean-Square Analysis of Discretized Itô Diffusions for Heavy-tailed Sampling. Journal of Machine Learning Research, 25(43):1 44,, 2024a.
Ye He, Kevin Rojas, and Molei Tao. Zeroth-Order Sampling Methods for Non-Log-Concayy Distributions: Alleviating Metastability by Denoising Diffusion. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024b.
Marlis Hochbruck and Alexander Ostermann. Exponential integrators. Acta Numerica, 19:209-286 2010.
Peter Holderrieth, Yilun Xu, and Tommi Jaakkola. Hamiltonian Score Matching and Generative Flows. In The eighth- Annual Conference on Neural Information Processing Systems, 2024.
Xunpeng Huang, Hanze Dong, Yifan Hao, Yian Ma, and Tong Zhang. Reverse Diffusion Monte Carlo. In The Twelfth International Conference on Learning Representations, 2024.
Aapo Hyvärinen. Estimation of Non-Normalized Statistical Models by Score Matching. Journal of Machine Learning Research, 6(24):695-709, 2005.
loannis Karatzas and Steven E. Shreve. Brownian Motion and Stochastic Calculus. Springer New York, NY, 1991.
Yin Tat Lee, Ruoqi Shen, and Kevin Tian, Lower Bounds on Metropolized Sampling Methods for Well-Conditioned Distributions. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, 2021a.
Yin Tat Lec, Ruoqi Shen, and Kevin Tian. Structured Logeoncave Sampling with a Restricted Gaussian Oracle. In Mikhail Belkin and Samory Kpotufe, editors, Proceedings of Thirty Fourth Conference on Learning Theory, volume 134 Proceedings of Machine Learning Research, pages 3993-3050 PMLR, 15-19 Aug 2021b.
Benedict Leimkuhler and Charles Matthews. Rational construction stochastic numerical methods for molecular sampling. Applied Mathematics Research eXpress, 2013(1):34 2013.
Benedict Leimkuhler, Daniel Paulin, and Peter A. Whalley. Contraction rate estimates of stochastic gradient kinetic Langevin integrators. ESAIM: M2AN, 58(6):2255-222) 2024.
Tony Lelièvre, Mathias Rousset, and Gabriel Stoltz. Free Energy Computations. Imperial College Press, 2010.
Stefano Lisini. Characterization of absolutely continuous curves in Wasserstein spaces. Calculus of Variations and Partial Differential Equations, 28(1):85-120, 2007.
Jun S. Liu. Monte Carlo Strategies in Scientific Computing. Springer New York, NY, 2008
Wei Liu, Michael Röckner, Xiaobin Sun, and Yingchao Xie. Averaging principle for slow-fast stochastic differential equations with time dependent locally Lipschitz coefficients. Journal of Differential Equations, 268(6):290--294) 2020.
Martin W. McCall. Classical mechanics : from Newton to Einstein : modern introduction. Wiley,
Hoboken, N.J., 2010.
13

---

Laurence Illing Midgley, Vincent Stimper, Gregor N. C. Simm, Bernhard Schölkopf, and José Miguel Hienández-Lobatte Flow annealed importance sampling bootstrap. In The Eleventh International Conference on Learning Representations, 2023.
Pierre Monmarché. High-dimensional MCMC with standard splitting scheme for the underdamped Langevin diffusion. Electronic Journal of Statistics, 2020.
Alireza
Tyler K. Farghly, Ye He, Krishna Balasubramanian, and Murat A. Erdogdu. Towards a Complete Analysis of Langevin Monte Carlo: Beyond Poincaré Inequality. In Gergely Neu and Lorenzo Rosasco, editors, Proceedings of Thirty Sixth Conference Learning Theory, volume 195 of Proceedings of Machine Learning Research, pages 1-35. PMLR, 12-15 Jul 2023.
Radford M. Neal. Annealed importance sampling. Statistics and Computing, 11(2):125-13, 2001.
Radford M. Neal. Slice sampling. The Annals of Statistics, 31(3):705 767, 2003.
Radford M Neal et al. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 2(11):2, 2011.
Maxence Noble, Louis Grenioux, Marylou Gabrié, and Alain Oliviero Durmus. Learned reference- based diffusion sampler for multi-modal distributions. In The Thirteenth International Conference on Learning Representations, 2025.
Kushagra Pandey, Jaideep Pathak, Yilun Xu, Stephan Mandt, Michael Pritchard, Arash Vahdat, and Morteza Mardani. Heavy-Tailed Diffusion Models. In The Thirteenth International Conference on Learning Representations, 2025.
Grigorios A Pavliotis and Andrew Stuart. Multiscale methods: averaging and homogenization, volume 53, Springer Science & Business Media, 2008.
Gabriel Peyré and Marco Cuturi. Computational Optimal Transport: With Applications to Data Science. Foundations and Trends@ in Machine Learning, 11(5-6)3355 607, 2019.
Angus Phillips, Hai-Dang Dau, Michael John Hutchinson, Valentin De Bortoli, George Deligiannidis, and Arnaud Doucet. Particle Denoising Diffusion Sampler. In Forty-first International Conference on Machine Learning. 2024.
Lorenz Richter and Julius Berner. Improved sampling via learned diffusions. In The Twelfth International Conference on Learning Representations, 2024.
Saced Saremi, Ji Won Park, and Francis Bach. Chain of Log-Concant Markov Chains. The Twelfth International Conference Learning Representations, 2024.
Katharina Schuh and lain Souttar. Conditions for uniform in time convergence: applications to averaging, numerical discretisations and mean-field systems. arXiv preprint arXiv:2412.05233, 2024.
Dario Shariatian, Umut Simsekli, and Alain Oliviero Durmus. Denoising Lévy Probabilistic Models. In The Thirteenth International Conference on Learning Representations, 2025.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems, 32, 2019.
Yang Song and Stefano Ermon. Improved techniques for training score-based generative models. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 12438-12448 Curran Associates, Inc., 2020.
Yang Song, Jascha Sohl-Dickstein. Diederik P. Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score- -asee Generative Modeling through Stochastic Differential Equations. In International Conference Learning Representations, 2021.
Gilbert Strang. On the Construction and Comparison of Difference Schemes. SIAM Journal on Numerical Analysis, 503):506-517, 1968.
14

---

Robert H. Swendsen and Jian-Sheng Wang. Replica Monte Carlo Simulation of Spin-Glasses. Phys. Rev. Lett., 57:2607-2609, Nov 1986.
H. Trotter. On Product of Semi-Groups of Operators. Proceedings of the American Mathemati- cal Society, 10(4):555-551 1959.
Mark E Tuckerman. Statistical mechanics: Theory and Molecular Simulation. Oxford university press, 2010.
Francisco Vargas, Shreyas Padhy, Denis Blessing, and Nikolas Nüsken. Transport meets Variational Inference: Controlled Monte Carlo Diffusions. In The Twelfth International Conference on Learning Representations, 2024.
Fengyu Wang. Functional inequalities Markov semigroups and spectral theory. Elsevier, 2006.
E Weinan, Di Liu, and Eric Vanden-Eijnden. Analysis of multiscale methods for stochastic differential equations. Communications on Pure and Applied Mahhmatics, 58(1 (5(11))15441188, 2005.
Keru Wu, Scott Schmidler, and Yuansi Chen. Minimax Mixing Time of the Metropolis-Adjuster Langevin Algorithm for Log-Concars Concave Sampling. Journal of Machine Learning Research, 23(270): 1-63, 2022.
Qinsheng Zhang and Yongxin Chen. Path integral sampler: A stochastic control approach for sampling. In International Conference on Learning Representations, 2022.
15

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract summarises the main contributions. In addition, we describe these contributions in more detail at the end of Section T Rach of our claims is clearly linked to the relevant sections of the paper for further detail.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
It fine to include aspirational goals motivation long it clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We outline the limitations of our methodology Section and provide more detailed discussion in Appendix F
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "I Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken low lighting. Or a speech-to-text system might be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed not penalize honesty concerning limitations.
3. Theory assumptions and prool's
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
16

---

Answer: [Yes]
Justification: We provide full assumptions for all the theoretical results. Furthermore, the complete proofs of our results included in Appendix C
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: Detailed information necessary to reproduce the main experimental re- sults sults-including the methodology and hyperparameter choices for both our algorithms and the baselines is provided in Appendices and D
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution
dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
be accomplished in various ways. For example, if the contribution a novel architecture, describing the architecture fully might suffice, if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, provide access the model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., the case of a large language model), releasing of model checkpoint, or other means that are appropriate to the research performed.
Depending the contribution, reproducibility
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily the architecture clearly and fully. new model (e.g., a large language model), then there should
model architecture, the paper should describe
(c) If the contribution either be way to access this model for reproducing the results way to reproduce the model (c.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e-g.. to registered users), but should be possible for other researchers to have some path to reproducing or verifying the results.
17

---

5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: The code is available at the following anonymised link https: huttps:/aaoonymouss //nonymus 4open science/ r/ samssimmnnngg barigg_bb___eeeee www.aavrrgging-44
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips... cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/ http:ss.cc/publigguuueeee
for more details.
The authors should provide instructions on data access and preparation, including how to access the data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The experimental setup for our algorithms is detailed in the main paper, specifically in Sections 3 and 5 Additional implementation details, including those for baseline methods, are provided in Section of if the Appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: All numerical experiments report eтror estimates in the form of standard deviations computed over multiple runs with different random seeds.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
18

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: The compute resources and runtimes for the experiments are provided in Appendix [D.4
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (c.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics httiissstttt::///tt
Answer: [Yes]
Justification:
We have read the code of conduct and ensure that the paper complies with it.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
the authors answer No, they should explain special circumstances that require deviation from the Code of Ethics.
The authors should make sure preserve anonymity (e.g., if there is special considera- tion due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification:
We address broader impacts in Appendix G
Guidelines:
The answer NA means that there is societal impact of the work performed.
19

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mecha- nisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
This paper proposes a general purpose algorithm for sampling from a probabil- ity distribution, which does not have potential for direct misuse or dual-use. We do use any scraped datasets in the numerical experiments which could lead to potential misuse.
Justification:
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
We properly acknowledge and cite all assets and resources used the paper.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset.
20

---

For scraped data from particular source (e.g., website), the copyright and terms of service of that source should be provided.
assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi patasets.papersvi thcode. dontsstsss has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged reach out the asset's & creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper introduces a sampling method for unnormalised probability distribu- tions and does not release assets.
Guidelines:
The answer NA means that the paper does not release assets.
Researchers should communicate the details of the dataset/code/model as part of their sub- missions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL include anonymized zip file.
14. Crowdsourcing and research with human subjects
For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Question:
Answer: [NA]
Justification: The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
21

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in this research does not involve LLMS as any important, original, non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs any important, original, or non-standard components.
Please refer to our LLM policy hattp:::/eerrips..cc/o/oeeeeeeee what should or should not be descmbbe..
/LMM) for
22

---

A Preliminaries
Multiscale methods and stochastic averaging. Before presenting a formal definition, we first provide the intuition behind multiscale methods, following Schuh and Souttar 2024] In many areas of science and engineering [Bertram and Rubin 2017 Weinan et al.) 2005 I. one often encounters systems that are difficult to analyse directly. However, there may exist related system that is more tractable, though it does not produce exactly the same behaviour. In such cases, the simpler system can be viewed as approximation of the more complex one. In our setting, the simpler system is slow-fast system, while the complex system corresponds to the averaged dynamics. We develop these concepts in the following.
To conduct theoretical analysis of the convergence of the proposed sampling algorithms, we examine the general stochastic slow-fast system studied in Liu et al. 20200]
dx; b(t, Xt, df + att, Xf)dBB dyy [f, Y/ )d + l Xt. xt)d
X = x E R Yo y E R.
(12)
where e is a small positive parameter describing the ratio of time scales between the slow component X and the fast component YŁ, and Bg and B, are mutually independent standard Brownian motions
on complete probability space (22, P) and (Fi, 2 0} is the natural filtration generated by Bg
and Bg. Let X, denote the solution of the following averaged equation
dxx- b(t, X,)dt + ott, X,)dBr, Xo 2,
where b(1,x)) = Sg blt, pt. (dy) and Pt, denotes the unique invariant measure for the tran- sition semigroup of the corresponding frozen process, which can be informally defined as the fast process with E = 1 and fixed slow dynamics X;- = F,
dy, f(t, x, Y.)ds xl,, Y.)dB' Yo y.
where B3 is Brownian motion another complete probability space.
We begin by presenting the assumptions underlying the analysis in [iu et al. (2020) and stating their main results, which will form the basis of our study in Appendix C.I
The assumptions the coefficients of the stochastic slow-fast system 122 are the following.
A3 (Coefficients of the slow process). (i) There exists 0 2 0 such that for any t, R R 0, E IR",
y E n with || < R,
211b((, v) b(L, 22, 1)ll æ2ll (10(4,21,)) 8) 22), - K. + " x2 1?
where (R) is an R4 -valued F; Fa-adapted process satisfying for all R, T, p E (0)
ar aT(RR) = L-" (((Hd) < co, on S2 (or(1)) <0,,
sup eE|o,T]
(1)I < cc.
Furthermore, there exists Ro 0, such that for any R > Ro. T20 0.
E S [K]( dt 80.
(ii) There exist constants 0z, 0g 2 and 7 E (0, 1] such that for any x E R", y. 1..11.10 U/ E Rm and T > O with 1,8 E [O,
[b(t,z, v1) - x,, v2)l < Crl - well [1 + I1 02 + K.(1) + [
and
[bbt, b(s w)ll s Cr 7|17 [||e] 1|=00100 2r Ω
where Cr 0 and Ze is some random variable satisfying EZG < oc.
23

---

(iii) There exist X1 2 0,C 0, 0, 2 2 and 0s. 0. 2 such that for any t> 0. z x E R", V
2(xz, b(t, x, v)) < (.(1(1 + || + A1 \1lloll
and
2,,t,,,,)), & (.()) + C[| llyll
[|o Kc(1)+ C||
4 (Coefficients of the fast process). (i) There exists 3 2 such that for any 2 O.
y1 U2 E R.
((, - x(1,,, - 12) (lg(t, n. v1) g(t, x, v2)l - --11 vall.
(ii) For T > O, there exists 72 [0. CT > 0, 2 1, 1, 4 such that for any 1, (0 and 24 E R", V/+ E Rm, i 1,2,
(L, v1) (aa)))((saaa v2)) )I 1 Cr((i - د/ + || - x2)l) + 21 -| + || + 11 2)),
(1,(1, 0(s x2. 02)l Cr a|T + || 221 + llve vell),
| - Cr(1 + || + |l)11"
so(t,z1,uill S Cr(1 + || + llval))-
(ii)) For some fixed 2 2 and any > O, there exist Cr, Sy such that for any & E [0, T.. x E IR", IRm"
2(y, f(t, a. y)) (k- 1)]lg((, x, w)1l S -lllll1 Aellvll Cr (ael(a+1)).
where X2 = 0f A1 = 0, and x2 otherwise.
Note that assumption A4 on the coefficients of the fast process guarantees that the frozen process has a unique stationary measure and that the fast process converges to exponentially fast. The following theorems establish the existence and uniqueness of solutions to the system in 12 as well as its convergence to the corresponding averaged dynamics.
Theorem A.1 (Theorem 2.2 Liu et al. 2020)). If assumptions A3 and A hold with A > 0, let EO = + Then for any E E 50) and any given initial values x E R", E RT, there exists
a unique solution {(x, ),), 2 0} to the system (12) and for all T >> O. we have (X; Y,,YF)) E C([0, x T): R²),, PP-almost surely. Furthermore, for all t E [0, T}. the solution (X,,Yf) given by
x + So b,, bs, ds So XJ)BB,
Yf y - f(s X Yf dss te So gls X dd
Theorem A.2 (Theorem 2.3 [Liu et al. 2020)). If dasssupptions AB and A4 hoold with X > 0 and k > 0g where 02 = max 20 + 2, 20, 40z, 010s, 201 01). Then, for any 0 <
2k B4
we have
lim E
sup [E(O.T]
|X -v,1=x111- = 0.
Optimal transport. Let (vr : R2 Edd) be a vector field and μ = (He)et(a,t) be a curve of probability measures on Rd with finite second-order moments. f is generated by the vector field v the continuity equation
dep4 + V ((Mr)) = holds for all E [a, b). The metric derivative of at E [a, b] is then defined as
Wa (a(me+..
|W : lim
[6]
If [A exists and is finite for all E Ja, b]. we say that μ is an absolutely continuous curve of probability measures. Ammbrosio and Kirchheim (2000) establish weak conditions under which curve of probability measures with finite second-ordd moments is absolutely continuous.
By A Ambrosio et al. 2008 Theorem 8. we have that among all velocity fields v, which produce the same flow ja, there a unique optimal one with smallest T.P (pt; X) X)-norm. This is summarised in the following lemma
24

---

Lemma A.3 (Lemma 2 from Guo et al. measures μ (p.))) any vector field (el,)c that generates i satisfies lid, S [|eelleanp) for almost every E (a, b]. Moreover, there exists a unique vector field ot generating μ such that |id, || (((()) almost everywhere.
2025 For an absolutely continuous curve of probability
We also introduce the action of the absolutely continuous curve (J4)EE(aM) since it will play key role in our convergence results. In particular, we define the action A(u) as
A(u) = S. 1A dt.
Girsanov's theorem. Consider the SDE
dX, b(X,, f]dt + (X.,,)d for t E [0, T), where (Be)e((oT)) is a standard Brownian motion in Rd. Denote by Px the path measure of the solution X (X:)rejo.7) of the SDE, which characterises the distribution of X over the sample space Ω.
The KL divergence between two path measures can be characterised as consequence of Girsanov's 's theorem Karatzas and Shreve 1991]. In particular, the following result will be central in our analysis. Lemma A.4. Consider the following two SDEs defined on common probability space F,, dX, ax(X)dt v2dB.,
dY, = bY )dt v2dB., E [0,1]]
with the same initial conditions Xo. Yo ~ po- processes X and Y, respectively. It follows that
Denote by px and p the path measures of the
KL(F* (P") 1 [xxEX (['lla())) - (x
B Sampling using multiscale dynamics
In this section, we provide a detailed description of the implementation of our proposed algorithms including the numerical discretisation schemes used.
B.1 MULTALMC: Multiscale Annealed Langevin Monte Carlo
Here, we present further details of the overdamped and underdamped versions of the sampler.
B.1.1
Overdamped system
In the overdamped setting, we propose a new sampling strategy based on the following stochastic slow-fast system
V
Voog
X VT-
dt v2dBe
if Aut < <
dx,
=
Aut
VV=
dt
2dB,
if ARa 2 A
Vhnt
Vdxt
dy; =
=
VV=
-ーー V Vo v
(x>ン*=] dt /
As discussed in the main text, the conditional distribution at. converges to a Dirac delta centred at 0 when Aet O, and to a Dirac at z when Aur = 1. These degenerate limits can cause numerical instabilities the endpoints of the diffusion schedule. To mitigate this, we define 0 < As < 1 and consider the modified dynamics
log
di + 2dB.
0 ARi As
VI-Aa
dXe
VI-Am
Vlog
VI-hxe
X
2dB
if As $ Ant <
Vhnt VV= Vine
di + +2dBB
if S Aet < 1 A8
VVa (Xi) dt + v VdB.
1 As S Ana $1
1 VV
v logg
Yi-Xi
+ / LaB. 0 Ame < As
dY.
V As
As
v
1->s
V1-A,
1
VVa
Yr
V log
Y=Xx,
dl + v Pab.. if As < An. $ 1.
VAni
y1 \T
(13)
25

---

Note that if the slsw-fast system converges to its corresponding averaged dynamics, then running the modified multiscale system up to time tg. where Axt, Ht, that is close to the target x, provided As is sufficiently small. For /. 2 ls. the dynamics in (13)
Ag. we can approximate the distribution
reduce to standard overdamped Langevin dynamies targetting я, initialised fúr, This warm start known to significantly improve convergence rates to the target distribution [Chewi et al. 2021 Lee
cet al. 202
Wu et al. 2022].
We now analyse different numerical schemes for implementing the sampler in practice. As we have just mentioned, when xxt dynamics, which can be discretised using any preferred numerical method. Therefore, we focus on the regime S Axt
A5. the dynamics in (13) reduce to standard overdamped Langevin
As.
For general base distributions v, we propose to use the following numerical scheme that combines Fuler-Maruyama for the slow dynamics with the SROCK method [Abdulle al.. 2018 for the fast dynamics. Before presenting the complete algorithm Algorithm 2 we introduce the necessary notation and coefficients for the SROCK update.
Denote by T, the Chebychev polynomials of first kind and define the coefficients in the SROCK update as follows
"1
Ta(wu)
1
(14)
g2
(w0)
and for all i 2, 2.....s.
2w0 1w Ti Ti(wo)
2w0 2
T:(wo)
Ks = P4.
(15)
Additionally, when (O, I), further exploit the linear structure of the score function for a Gaussian distribution to design an efficient exponential integrator scheme. Below, we derive the corresponding update rules, distinguishing between the regimes 0 I Axt < and À I Ant S 1.
Regime 1: 0 S Ane < x. Define the modified schedule
AAs
iro S Aus < As
Not
Anne
is S ARt $ 1.
The slow-fast system
13)
simplifies to
dx,
Ye Not
dl + 2dB
dY
VV=
Ye-Xe I-Nas
de + VzdBe.
Vx.e
The exponential integrator scheme is then expressed as
0a
X₂ dt
v2 v
0 dB. d
d
Y. - X
-
Ye
/v
Vv.
Y
dt
VE-T
Ia
where given time discretisation < To < < Ta of the corresponding time interval, we define t_ T--1 when E [T)-1, Ti). The explicit update rule is then
X++1 V1+1
Xi
A (11,)
)
Ao(Ti, Ti+1) & la
@ Ta
x Az(Ti, Ti+1) @ Ta &
VV=
where & N/0, I2a) and
(-(-)- --) 1' (1 Noww -du du
A)(a,b) = exp
Ar(a,b) = 1.' Vie Ao(u, b)du
b2(a,b) = 2/ Ao(u, b)
(0= 5-1) Ao(u, b)
26

---

Algorithm 2 MULTALMC sampler: overdamped version
Require:
Schedule function A₁ value for As and À, number sampling steps L, time discretisation To <... < T' 1/k, step size hy Ti+1 -T, Constants for SROCK step from (14). 15
Initial samples Xo ~ (N,,,) r), Yo ~ NO0 T). Define the schedule
Aa Autt
i Aut As
Xse
it Ao S Anee 1.
for 1 0 L do
(1) (2) (O I)
0 < AR, Ag then
EM for slow dynamica
X1 V--Xe
V logu
N.-Y h + 221
if Man
X1+1
RT
X1
VV=
Yi Mir,
lu + (2 (1)
if i 2 Ā
SROCK for fast dynamics
K1,0 =Y;
K1,1 = KLD + c
h1
(- 1 VV
Kio0 Vinzz
+ Vi-
V Vlogo
K1oo-1111 VI-
Kis
ly (- VVr
Ki
+
V Vlogv
[Ku-1-Xi411
))
+is R:Kei-2
Vor,
Y141 = K₁ + v 2h (2)
if - ds S ART. - then
X14 X VV (Xi) + 22hie")
end for
return Xz
Regime 2: á S Ant < As-
Here, the dynamics become
dX,
V
Yr
di
+ v2dB,
Annt VAni
Vhat
dY;
VV
Y--Xe dt + /dd.. T-
In this case, the exponential integrator scheme can be formulated as follows
VV=
X1
0 0
X1 Y:
+ v2 6
ddB. dB,
Y.
-
di
Vinc
VV
Y /
VE-I @ I.
The corresponding update rule is
VV.
Yi VATi
Xx+1 Au(TT, Ti41) @ la
Xi
AA Tim1) 11 1a
Au(Ti Ti41) 8 Ia Es.
Yi41
Yi
VV
Yi (Ar
where & N/O, Ied) and
(- (-(-0-
0
1" (1
Ào(a,b) = exp
g-1
- du
27

---

Ai(a,b) = S-" Vina
__bo(u, b)du
As(a,b) - / Ao(u, 6) (o ( Ao,u, b)rdu.
B.1.2 Underdamped system
Similarly in the overdamped case, to mitigate numerical instabilities arising from the degeneracy of Pt. when ARt and Aue 1, consider the following modified dynamics.
dXe = M Vidt
Vloge
-V di VTdd if 0 S Aut < As
1-28
VI-
V Vloge
-
1V- dt VTTdB.
if A Aut <)
dV.
vI- Auc
VV=
IM (M-'v) di + v2TdB,
if Á I Aut < 1 20
Vince
VAgt
VV= (X+)
rM
Vo) dt d2IdB.
if 1 As Ant $1
VV
v Vogg
dt v Ldd.. if 0 Aat As
dY
V
V log
dt + PaB. if Ao Aut 1. Xx
(16)
Note that when Ane 2 1 As. the dynamics in (66) reduce to standard underdamped Langevin dy- namics with a warm start, which can be discretised using any preferred numerical scheme. Therefore, we focus on the regime 0 < Axt
As. We propose use hybrid discretisation method that combines the OBABO splitting scheme [Monmarche] 2020] with an SROCK update [Abdulle et al. 2018) for the fast-dynamic the matrix square root of the friction coefficient, which is well defined since a symmetric positive definite matrix. That is, VT is any matrix satisfying VTVT
The algorithm is explictly defined in Algorithm 3 where vT denotes
It is important to mention that for implementation, we set the mass parameter to M = I and consider a dimendependent friction coefficient I' to control the degree of acceleration. Additionally, the same discretisation scheme applies to the heavy-tailed diffusion, with a modified fast process defined by Eq. 77)
B.2 MULFCDIFF: Multiscale Controlled Diffusions
B.2.1 Overdamped system
Diffusion models typically consider a stochastic process (X, (ejoo constructed by initialising X the target distribution æ and then evolving according to the OU SDE
dx,
- dt + v2 dB;.
(17)
Under mild regularity conditions are described by the following SDE
Anderson 1982], the dynamics of the reverse process (x
()
dx.
(x 2V log (xx))-((X))) dt + /2 dB.,
(18)
where q denotes the marginal distribution of the solution of the forward process given by
(z) / ((v))r
with Pr (zly) being the Gaussian transition density Nie 'y, where of = - e 2t, Using this, the score function log (Xx, can be expressed as
X -
x [Y]
v log qT -
Since the score term involves an intractable expectation, we approximate the averaged dynamics in 18) using multiscale SDE system. The fast process of the system is modelled by an overdamped
28

---

Algorithm 3 MULTALMC sampler: accelerated version
Require:
Schedule function A₁ value for A5 and A, friction coefficient (scalar), mass parameter M (scalar), number of sampling steps L, time discretisation 0 = To ... Tz 1/k, step size
h; Ti+1 - Ti. Constants for SROCK step from (4). (15). Initial samples Xg NNO,,), Vo (((O,TT) Yo NN 1). Define the schedule
[25 A8 if ARc < A5 Aut
Xxt
As - Ar $
for 1 = 0 to do
41 52 51
(2) (3)
1)
0 & ART < - As then
Salf-step for velocity component
1 h T' Vi VhTTe
V
Vi + V/i V logg
Xi-Yz Vi-
if Men, < X
V,"
Vi
VV
if Xem, 2 Ā
2 2./le,,
Full EX-step for position component
X(+1 Xi h,M h V11"
Full step for fast component
Ki,0 Yi
K1,1 = KLD + uI hi
(-VA.o,
(-1VV,
VVa
Kio Vir,
+
v log
VXe,,
VI-Mo
VI-Mas
Kiss = u
b
VV-
+
V Pogg
Ki 1
))
+ MiKii-i + R:Ki--2
Y44 = KN +
(2)
Half-step for velocity component
V1"
h
V Vlgg
if io
Vi"
Xen
V."
VV
M441
if if 2
hi 1'M
Vi+1 =
1
vi"h + Vhil'e!
2
if - As S AxTi then
Half-step for velocity component
V/' = [ 2rM Vi (hr
V1" V/ VV. (Xi)
Full EM-step for position component X X1 + h,M
Half-step for velocity component
Vi"! V." 21 VV, (Xi+1) Vi+1 = (1 4 vin + (())
end for
return (X,, Vi.
29

---

Langevin SDE targetting the conditional distribution Y| Y|X,, which takes the form
dy, = _ (log + Vy, log (X.1))))))) dt + v 2 ZaB
where
Vv, log pr- Xx. V.Ik. Y: e" 2(7-0)
Putting everything together, we obtain the following multiscale system of SDEs
2 dxx = (x.- og-- (x--e-1T-11)) dt + v2 dB.
(19)
dY dY = Vlogg -15-- -7 (Vi (V--e"- xx)
+ v dd..
This system of SDEs can be discretised using either an exponential integrator scheme or a numerical scheme that combines Euler-Maruyama for the slow dynamics with the SROCK method (Abdulle et al. 2018) for the fast component, similarly to the discretisation of MULTALMC in the overdamped case Appendix [B...I.
B.2.2 Underdamped system
When using an OU noising process, the forward underdamped diffusion has the following form
dx', M M-VV.d. dt dv. (-x, V. di + /2r aB,
which can be written compactly
d d A La 2, dt B La dB",
where Z,=(X, (X,, Vi) and
M-1 -
0
=
B =
V2T
where, as in Appendix 1.2 VI denotes the matrix square root the friction coefficient I, which is a symmetric and positive deñnitt. The solution of this SDE is given by [Karatzas and Shreve 1991]
7. A 2 J,¹ (cAitt=))) B)
@ Ta dB
As discussed in the main text, there crucial balance between the mass M and friction coefficient I. We focus on the critically damping regime by setting 1 4M. This provides an ideal balance between the Hamiltonian and the Shnstein-Uhlenbee components of the dynamics, which leads to faster convergence without oscillations [McCal] 2010). Under this setting, the matrix A can be written
A (-
4 -
Then, the matrix exponential cA has the following form
Al ( (1
2
4r-2 1-21
mi. m2,t ma,1
e
m4.1
On the other hand, the reverse SDE is given by
-
aV. (x.+4r-7 V 21' log 1o (X.V.) di + /2TdB.
where ai is the marginal distribution of the solution of the forward noising process which has the following expression
qe (x v) = m())o(ow))) aly, ou) dy dou,
(20)
30

---

with & [0, I) and conditional distribution Pr(x, uly, vo) given by a Normal distribution with mean
m4ly, vo) = c¹ 8 Ia 10 =
(mI.y + m2,
m3,ey m4,,000
and covariance matrix
E11. E12,3 Exr
80 Il (fr
el(t-s) BB [ealt-o]] ds @ la.
Σ, 8
E2a,t)
which results into
E11,1 (cur-) - 411 8 -er-
E12,1 = E21, = 4
2
411- 1) [T 21² -40-1
Σ22,1
We note that to in Eq. (20) can be integrated out analytically. Since both pe(ar, ulu, vo) and (p(1)) Gaussian, follows that (a(x,, ely) is also Gaussian with mean
mm.yy ma.1y)
nk(y)
and covariance matrix
Eun Ee1
Ei2 Σ22,
[Eii, E12,0
m2,,
m2,574,4 ) Ia.
Σ, ② Ia =
@ Ia =
E21, E22,1)
+
m22
We can rewrite Eq.
200 in terms of (=( vly) as
q((,,) v) x x(y) (z, vly) dy,
We note that the conditional distribution pc(u (u y) remains Gaussian, with mean ma v) 133,11 + E1g, E (x m).,9))
which is linear in y and covariance of where
of S22,t E12, 24,, $1,,2111
This provides the following expression for the score
V. ge(o, v) = V. log (c)))) V By | (z, Y)
V- ft (2, (Yln.
o?
of
where
fe (r (Y) ma [Y]+ Enn (x (y,
Using the derivations above. we define multiscale SDE system that enables sampling from the target distribution via the reverse SDE, without requiring prior estimation of the denoiser Ey [Y] We consider fast process governed by an overdamped Langevin SDE targeting the conditional distribution Y given the current state (Xe, V.). which takes the form
dY = 1 (V (Vlog((H)) + Vy, log poogr---x xx,, v.1r)) dt V dd..
where
Vy, Pr Xx, V,jx) = =- [m(
my,T (
(x., V.) LT' mL.T-()
-13
Combining altogether, we arrive at the following multiscale system of SDEs
(ax, = -1 dt
dv,- -r 17-1
(V. (Ma.T.YY E12T 8--1 (X. mir-
v2T dB.
dY ? (VlogT(Y.). + Vx, loglogge , F.(x.)) dt + VidB..
(21)
31

---

Inspired by Dockhorn et al. [2022b] we propose novel discretisation scheme which leverages
symmetric splitting techniques [eiikuhllrr and Matthews [013 justified by the symmetric Trotter splitting and the Caap-ealmbbbllaaae The 1-th iteration of the proposed numerical scheme consists the following steps, with by denoting the step size.
formula [Strang] 1968 Trotter 1959] [Tuckerman] 2010
1. Evolve the following system of SDEs exactly for half-step h
dxxx - dt, a = (Xx + 4x.+4r-117 - V.) di + 22r dB..
(22)
2. Evolve the velocity component (full step) under the following ODE using Euler's method
dV.
o (V.-
2r
(ma.r-ei + E1a.- (x (X.- mi.r-iXi)) dt.
3. Evolve the fast dynamics (full step) using SROCK method
2 dB. dy, (v (Y)) + Vv, log pr (x,, V V.11a)) dt +
4. Evolve the system Eq. (22 exactly for half-step h2/2.
We summarise the steps above in the form of concise algorithm. To do so, we first derive the exact solution of the SDE system (22 used Step Define the matrices
41 -1
0
0
1r-1
B
D
2T
The matrix exponential e At takes the form
2
-
e 12
t
+ 2
The solution to the system (22 at time t, given an initial condition (Xo, Vo), can be expressed as
(X,, Vi) (Xo. Vo), Ia)
where the mean and covariance matrix are defined by m... Vo) e @
Xu Vo
Σ EEn, Σ12, = L.' f(a Ait BBT [cdle-o)] ds. (23) Eea,))
La
Em.t
The explicit expressions for the entries of the covariance matrix Ex are
En.t (-ear +1 41 + 8t21 -e--
112 = S211 = 41.p-11 -eer
L02,
2
1 (ar
+
26 a0r-
4
The full algorithm is presented in Algorithm
Challenges for heavy-tailed controlled diffusions. When using heavy-tailed noising process given by the convolution with Student's distribution with tail index a, the forward underdamped diffusion takes the form
dxx = M V dt dv.
atd
VV. dt /2T dB, 1X
X
The non -onelrncar drift term in the velocity dynamics prevents from obtaining analytic solutions via standard techniques for linear SDEs. As result, we cannot directly reproduce the computations used in the case of an OU noising process. detailed investigation of this regime is left for future work.
C Theoretical guarantees
In this section, we provide detailed proofs of the theoretical results presented the paper and discuss the additional challenges posed by heavy-tailed diffusions.
32

---

Algorithm 4 MULTCDIFF sampler
Require: Friction coefficient I (scalar), M r2/4, number of sampling steps I,, time discretisation
= To < ... < TT = T, step size by = Ti-1 T. Constants for SROCK step from (1). (15).
Functions (:)) and M, from 23
Initial samples Xo ~ (OO, Vo ~ MT), ~N0,
for / = to L do
(X1+1/2, Vi+1/2) Niiiii Vi), Enn
2r
(Vi+1/2 (a.r-n Σ12,7 Sui.-n (X++1/2 mi mn))) h
V441/2 V/+1/2 -p
2
KL, Yi
KL,1 Ki.0
+ Mn (-
Vr ( (81))1122//+))))
Kis
=M by (-VV, (Ku,-1) log (or
+ K1K--1 + K,Ki
:
Yiu1 = Ko + t & ~ (0,11)
(X441/2 V/4122) ~ Ja Vi-1/2), (2)
end for return (Xz, Vt)
C.1 Convergence of the slow-fast system to the averaged dynamics
Based on the stochastic averaging results from Liu et al. 2020 presented in Appendix A we restate our convergence result, Theorem 4.1 and provide the proof.
Theorem C.1. [Theorem 4.1 restated] Let the base distribution v ~ N I). Suppose the target distribution æ and the schedule function A satisfy assumptions 47 and AQ respectively. Then, for any E Eo), where En is specified in the proof, and any given initial conditions, there exists unique solutions {(X;, Yf), t
0}. [(x, VI, Y,Y, 2 0} and (Xx,, V (Y)) 2 0} to the slow-fast 16 and (10). respectively. Furthermore, for any p > 0, it holds that
stochastic systems
lim
sup [E|O,]]
||X
-x11") =0,
where X, denotes the solution of the averaged system and X the corresponding slow component of the multiscale systems
and (10).
Proof. If the coefficients of the proposed slow-fast systems 6 and
10 satisfy assumptions
3 and A4. then the result follows directly from Theorems A.1 and A.2 Thus, we verify below that
these assumptions hold in all cases.
It important to mention that although the theorem specifically addresses the case where A I). assumption A only requires that V log v is Lipschitz continuous with constant L, and has maximiser x². Without loss of generality, we assume x* = 0, i.e., V log v0) = 0. Therefore, we show that assumption A3 holds for general base distributions satisfying these properties.
Besides, for the systems based on annealed Langevin diffusion (4) and (6). we instead consider their modified versions (13) and (16), respectively, these reflect the dynamics used in practice during implementation. We observed in Appendix B.1 that these systems reduced to
A multiscale system with time- indeieddeenddnt coefficients for 0 S Axt < Ag-
A multiscale system with time-dependent coefficients for As S Awt < 1 - As- Here, we distinguish
between two different dynamics depending whether Ane Ā or Are 2 X.
Standard Langevin dynamics for - \ Aut S 1.
33

---

We analyse the first two cases, as standard Langevin dynamics is known to converge to its invariant
distribution when initialised with warm start [Chewi al. 2021 Lee et al. 2021a Wu et 2022]
See Appendix B.I for further discussion.
We denote by m² minimiser of the target potential V
Overdamped MULTALMC
13
The coefficients of the slow-fast system for Aut [0, 1
As] have the following expressions
logu
x-n Vi-A5 x-v vT-Ane
if Aae < As
b(t, x, 2)
Vlogv
As Axx Rx if v Aat < 1 A8.
-I-
1 VV
VAnc
Axrc
o(t,z) = 22,
V logu Vi-Aa
1 VV
if0 S Ant < As
f(1,2,))
Vi-Aa
A
V
(24)
Vi-
V log
(vT
Vial
VVa
is S 1
g(t, V2.
We consider two regimes based on the value of ARr= Regime Ane E [0, ÀJ and Regime 2 Axa [} >s). We analyse the convergence within each regime separately.
Substituting the expressions of the coefficients, it follows that assumption AB holds with 01
02 = 0g- = 0g- = Og = 1, 0g = 2 and Zr ||m= X/2 is constant random variable in both regimes, shown below.
2 (bb(t, v) b(t, ,)|||> x2l [l((1 x(2,22))1
(05Am
24z ||+1 - xall + (AsSAR <5) T-Ane
26e ||-] x2ll²
Regime
0
Regime
2Lv
$
|::1 x2ll 11
(ussrrsss <ds) 1-2o
1111 vall + (ASSA <d) x-x 11 yall
Regime 1 Regime 2
16(t,,,v1) - b(t,,ttt,..z w2)1 <
1l [
Ily1
22ll
Aus
Ly
Lx
< max
11-X
Illn s2ll-
X
For the following, assume without loss of generality that [bb(,, - x,
Ly
V log
VI
VI-2a
Vi-x
(vi
Ln | + || logg
(x= )l)
Hlsssteedenc vinnne
Regime
(iSharnrooe--o) Tane
La Iill
VV, (otn
Regime
Vhcs
Cle-s/1 In (||]| + (!yl)) (1-1)2
Regime 1
0pe-pm1 La || + Im 1 23
Regime 2
< C' max {
Ly Lo }
|l II VX
(1 XI= X2 } le s
|| + ||a|||| +
2
2(x, b(t, x) v))
34

---

(osAe<As) v1-2s 2|z| VV=
2||=|
V log
v1-2s
2|1r|| 1l (ASSAarc)) vi-An!
V Vloge (VT-Anc
Regime 1 Regime 2
(ASActEl-As) Vicc
VAni
1-
Lu
+ (2]1a|[++|1b11²)
Regime Regime 2
(21/z| + ||ull + Al1m=*13 11²
Lu Lx
{
Lu Lx 1
$ max
max AlImn 11') (1 + |--l²) + max
11 11- 3 ||a1l?
1- X
which leads to the constant Å assumption A3 being given by A = max { { (x }.
Regime Regime
[]ttt 2., v/ll 1
lbll + VA m')l
L m|m" I| + max
L 1- [[ + Ilvll),
L
[t](1, 2.
On the other hand, to verify that assumption A4 we recall that for a standard Normal distribution V log v(x)= = =r Besides, note that we can express the drift term of the fast process, f(t, x, y), in compact form by rewriting (24)
VV, ( ())
f(t,,,
V
Vog
VI- X
(1
Nai
V
(Xai
where the schedule Xxt is defined by
A if < Ane < As ARt
Xxe
if As S Aut 1--a Aa-
Therefore, we adopt this compact form in the following analysis and omit explicit distinctions between regimes to simplify notation. Using this, we have
2 (1, v1) 11, M2), 91 wz) + (19(1, - (I, N2)11
1 X/t
11y1 3²11²
VNxs
VV,
111
VV=
3/2
21
8/2
s (1 1 Net
2Mx
+
1181
32ll2 < Bblly1 Mell²
No.
1
1
(llvll + || + ||:11 - No
xzll
[f v) -f(s, F2, vi1l 5
1- N-xa
1- No.
1
Txlllvll
VV=
y
VNa
VXaa
VXxa
VXae
<C1.slt s/m
[llvll || + C2. x2ll + Ca.slt - s| lbvll
<Co(It s|"h
||
m2ll))
(1 + || + Ilvl)),
L
1
La Amllm"]
al) max
(1 + ||-||| +1+1|=ll
+ (bll - Cs (1 + || Ilvl))
(1e Xxa)* 1
Nar
Xxe
(1t 21, v1) - g(s, 22, 02)1l 00 2o(t, v)ll = V2.
Since the potential V, is strongly convex with constant Mr, then V. satisfies the dissipativity inequality (VV)(x), 2 a |2:1 - bm. Therefore, we have that for any k 2
2(y, f(t, x, v)) (k- 1)llg(t, w)1 S-2 5
V: Not
y x
VV.
y
+ 1)
Vxa
VX
2
2a.
(llvl12 - (y, x))
|///l
+ by 2(k
X
X
35

---

20x
110112 +
1 - Nat ||//12 + by + 2(A
1
1- Xxx Xxx
S delly/l + Cu.o (11xr1 + 1),
where X2 is given by
20z
A2 1 As
This concludes the proof for overdamped MULTALMC, where the constant en the theorem statement is explicitly given by
20x
A2
60
max
1-X
Lx
Underdamped MULTALMC or 16
In this case, since we introduce an auxiliary velocity vari- able, the slow component consists of both position and velocity, (z, v). Accordingly, the coefficients of the slow-fast system are given by the following expressions
M 0
if S Ant < As
log
TM-
1-
1-18
M M
b(t, (2)
is S Ae
V Vlogv
v
vI-ANE
x-N
-M
M-
1
it < - As.
VVz
#)
IM-
Vike
Vhxt
(0
0
0
ott (2, f)) =
la.
0
V2r
v logv
y-z
VV
if S Ane < As
[(1. (2,)))))
V1-25
VT-A8
Tx VAnt
log
y-z VI-Axt
VV=
As S Ant < -- 2d,
-1-
git, (x, v). v) = 22.
Since the coefficients of the fast process the same in the overdamped setting, then assumption A44
also holds the underdamped case. To verify assumption A3] we adopt the same two-regime framework defined in the overdamped case and analyse convergence within each regime separately.
2ljb(s, (11, 21).2) dér, (22, a) #1) (#2, #2)1l [l0((t a2t)
A -
v2)
2
(osAnesas) vi-sa (Vloge
-リーツ Vi-As
logu
<イ=))a ))
TM- 122)
[[(E,, vi) (02, v2)ll Rogime
+ (((AsSAnes)) vi-zne
V log
11-8 VI-Ane
log"
12-y VI-Ane
Regime
213
$ M 2 + 21 llv1 - 1//12
(1 x)2 || sall [((2))0)) (2)) v2)
$ 2max M-1 V vZLu
[[( (22""2)1"
[[b(t, (z, v), b(t, (2, v)) y2)1l
10oSAweeae VI-Aa
VVlog
m- VI-ho
logv
1-12 VI-Ma
0
1 VI-Aer
Regime
Vlog
2-
V llog
2m12 VI-Aot
VI-Aw
(VV,
-
VV=
v2 Aiui
Regime 2
36

---

Lx) ) l1 yell- 11 X
Lv
max
For the following, we can assume without loss of generality that S t []b(t, (z, v), b(s, (2, v),
(n البلبريع
V logu (vi-
Vi-aas
vi-No
vi-X
+1 (ANSSAnai vi-hee
vi-der
In العا + | Vloge (vi
Regime
V1-A
I (ASAatiio Vice
Thce
Lell
VV= (Ttee)) )l)
Regime
C]e-a|71 Ly (|1/| + (1vl))
Regime 1
(1-2)²
Opt-al11 La ||bll + lm
Regime 2
32
by
Tix
||=" / V
s C'mxx (1 (1- 2/22 x2
||v|| + ||-|||| +
2
2((x, v), 61, (x, (c,0),)))
2M -1 (xr v) - -2 "TMM- + IOSAn
2||l|L
v logg
z-ツ
V1-26
Vi-28
+1 (A -x V log
2|vl
z-y
Regime
2M (x, v) - 21' -0 ||||| 2 + VAni
2101
VV
Regime 2
M M (1>=]] + 2 + Lo
1-x
(2||z|l + ||]||| lbll
Regime 1
M (|==) + ||||| - 2MM 2 4 [21101P ||v|| + X||m I1²)
Regime 2
max MM- 2Lμ 2L.
L
Lu
-1-x max Xl1m=111 (1 + |||| + [ell + max
[]6(, (x. will
Vlog v ( 1|"
M M-2 4-2(1+213]11111 + 1Ossnn 1-
+ (A 1-3
V logy (Vi-Aet
Rogime I
M-2(1 + 2|-)||=1112 +
vv.
Regime 2
2| < max M 1u-ivi+ 2T
1-X
2Ln
11(2, w)ll
2L
V
max
1- x
Lu La ]
livll
[|m max ss
1/11+
21
2hr
21 } (1l( w)1l llol))
1-X'
(t) = 2T.
This concludes that assumption AB holds with 0 0, 02 03 = 03 = 0, = 1, 0₂ = 2 and
ZT = |l VT/2 is a constant random variable in both regimes. Besides, the constant Eg in the theorem statement is given by
Aa
24
As
A1
L
max
Underdamped MULTCDIFF (10) or 21 From Appendix system are given by the following expressions
B.2 the coefficients of the slow-fast
-
b(t, (z,, v).v) x x 11
2r (u- (ma.T- (w-(mar-11) + S - m..T-1))
(0)
0
alt, (z,
@ Ia.
0
v V2F
37

---

(ML,T-- mg
(z, V (v) mms. (x,, u)]
g(t, (a, n), v) V2.
where m1,7-t+ -3, Er- are defined Appendix B.2 Substituting these coefficients, it follows
that assumption holds with 01 = 0, 82 = 03 0g = 0g = 1. 01 = 2 and Zr is a constant random variable, as shown below.
(), (m11, 21).2) ber, (m2 12), ( (ill - (z2. ا2)د4 + ||or(e, (21, at, (#2, m2))]l
(nr - x2) (1+ - - -1 - 12)
2
+ (e) vz) (11-111-11 (1-1- - _3_1/
[|(21, (22.
1+ 2r Star S1ir )) |lo1 2212+ 16T
2
))
2
-r-1
o7
1|01 -e212 (w1, v1) e2]1l
$ S.ll (1) - (2)1=
[b(t, (z,
- bb( (z,
o. (Ma.T- + E 1 ir m..T..) -. (y) - wa)) s Collyn vell,
[|b(t, (x, v), v) - (x, v),
$
1
||
E E -11
- 55
11,T
||x|||
xf or-
of
07-1
ma,7-1
-12T-1
E12, 11, o7-
M1,7-
|l
o7
I Cr/l 3/at (llvll [1(x,)
where we have used the expressions in Appendix B.2 In 10 such that
which imply the existence of constant
E1ar- 11.7
o
o7_1
07_
m3 112
-
ma.T
E
5 11,
T
m1.T-5
S Crlt 3/31,
o7-
of
2((z, v), (1t,
(x,, v).8)1
22
2 ++
21
S12T 1 (2.0) +2 24
2r
||n11
o3_.
of_.)
41
(MaTT-1 E1ZT - -)) (0:0)
o
S max
| 1--
2r
112r (
2 4r
21
2r
(ma.r-1- - E.2. - - :) 11(r, 1111
of
a7_e/
21
(Ma.T-1- E12 5 m1
|1al 2 C. ((,, w)ll A.llal
[b(0, (x. v).
$ 3
1
+
21
S12.T-
lall2
443 16r +
|1
12r
(ma,r-e
B12. S mi -m 11
o__
o7
of
5
V3
2r
S
L12T-
1
0-1
ar
v3
4r
2r
[1 [[(x,]]] +
2v3r
ma,r- 21 - llvl
07
< []](,,0]]
||vl])
38

---

2)) -rr.
We finally check that assumption A holds.
(ff(t, (c,, v), ML) - ftt, (2, v), 32), yn 82) + llg(t, (xt, v1) g(t, (2, v), 82)l
-V2 (v1) VV, (U2) 11 v2) -/ ma, -. ) (mi.T-t
ma.r-1 11.1
82ll 9221
- 5 + (mirr- + m Amit (ET 11vr - y2ll S - B B1y
wy
where we have used that, for all I, the precision matrix Σ is positive definite and Amin (2an(())))))) denotes its smallest eigenvalue.
m1,T-
(Mi,T-0)
[I, (11, v) f((: (122 02),
11vll
m3,T-1/
+ v1) (1z, v2) - (m)r=r) [1 v2)l
E (m1.T-)
(mI.T- m3,T--
(M..T-))
I Crlt- - 8 22 + Cr [[(21,, 01) - v2)) + Cr [ u2)1l
s/²²
Cr (It [t- 8| + 11 01) (xz, ua) (1 + x 02)ll),
where the existence of f/2 E (0, 1] follows from the expressions for -,, ET- provided in Appendix B.2
[f(L,, (z, vill & V log (v)]l (m.,r- m, .)5---1 (m).T=) ||
(c, u) (mm))=))
s L=((llll + Imm III + //mmirrr + m Amax (ST-t) -m, + m + ([(x,))]]
Cr (1+ |vvl + v)),,
(lg(t, (m) 01), 31) g(s, (2), v2)ll = 0, 11 (zz, v). v2.
Recalling that the strong convexity of V2 implies dissipativity inequality, we have that for any k
2(y. (4, (z, (),,)),)) m)) + (A= 1)]ll( v). -2 (v. VV) m....-
m13.T-1
+ (z, 0) (mur-:)) .) - 2(k- 1)
-2aaallvll* ba) - (m1,r- Ms, S
mI.T-t ma.7-t)
||v|l + 11(,, o)]l" || 2(k 1)
s -2 |ull + (r,, r)ll" + 1),
which concludes the proof.
C.2 Challenges of extension to heavy-tailed diffusions
When considering heavy-tailed diffusions, that is, when the base distribution v is a Student's t distribution, the stochastic slow-fast systems (4) and (6) could still be used derive sampling algorithms. However, Theorem showwinn convergence to the averaged dynamics does not apply in this setting, as it relies on exponential ergodicity of the frozen process. In the case where the fast process follows overdamped Langevin dynamics targetting a heavy-tailed distribution, convergence is sub-exponential or polynomial rather than exponential [Wang 2006 Chapter 4], and thus the frozen process fails to meet the required condition.
instead of using overdamped Langevin dynamics for the fast process, an alternative is to consider different diffusion processes that target the conditional distribution Pt,z more efficiently. Motivated by the work He al. [2024a], we propose employing natural Itô diffusion that arises in the context of the weighted Poincare iecquality, which has the following expression
As discussed in Section 3.1
dY, = -- VUA. (v)dl + / 20pxx dBr,
1
(25)
39

---

(B) (v++
||= y/1 a(1 Ant)
atd
-(a+a)
X (v) (V) (atd)
Vxxa
where a denotes the tail index of the noising distribution v. This results into the following slow-fast system in the underdamped setting
dX2 M Videt
{
( Vlogu
X VI-
- -M dt V if Aut <
dV.
(26)
Vini VV
TM dt + VVTdB, if Aut
dY (a 1) x, (y)de V m dBt.
However, although the modified Itô diffusion offers improved convergence properties (He et al. 2024a]. it does not guarantee exponential ergodicity of the frozen process under mild conditions. As result, assumption A44 is not satisfied in this setting. We leave the analysis of the convergence of 26 to the averaged dynamics (5) for future work.
C.3 Bias the annealed Langevin dynamics
We formally quantify the bias introduced by the underdamped averaged system diffusion path which is given by ((=)xe(0_1/k) with
relative to the true
1 M log (e(x))
p=( v) exp
(27)
Theorem C.2 (Theorem 4.2 restated). Let QU-ALD = (9.U-ALD) be the path measure of the
diffusion annealed Langevin dynamics 5). and Q (De) that of reference SDE such that the marginals at each time have distribution Pe defined in (27) If go,u-ALD = Po. the KI. divergence between the path measures is given by
KL (Q [IQU-AL))
1 A(p) 16 A(e).
Proof. Let Q be the path measure corresponding the following reference SDE
dZ, (M-L0, (Z,, U.)) dt
fort [0, 1/x].
dU (V log (Z.) TM "M- dz. U.)) dt + v2rdB.,
The vector field i ((01.1, 02 (2, is designed such that (Z,, U.) P. for all E [0, 1/x]. Using the Mikker-Plancee equation, we have that
a... V= (pe (M (M- 01.)) V. (De (V log TM '' " 02.:)) + TT... V (pi01,)) - V. (f102.), LE [0,
This implies that ê; (01.e, (2,1) satisfies the continuity equation and hence generates the curve of probability measures (pr).. Leveraging Lemma A.3 we choose i to be the one that minimises the
(f,))) norm, resulting in 06 (2()) = j being the metric derivative. Using the form of Girsanov's theorem given in Lemma A.4 we have
KL (Q (LQ)
1
Ba
[oa 1/N1
[1,, FO]l [o)(
dt à [ [|,, Folli dt 1 dt
1// 1/111100
2// 1:2 1
1//0 p).
4
x /''
bolt dt 1 A()).
4
where we have used that x |Rl, and the change of variable formula.
To conclude, we note from that the position and the velocity variable independent, and that the marginal distribution of the velocity e remains constant along the diffusion path (ft) [o,1)_As a result, the metric derivative simplifies to
Wa (2(pess.p))
Wa(He48,
|p|,
lim 8-0
lim 6-+0
[
8
40

---

Consequently, the action A(p) satisfies
A(p) A(H),
where μ (14)) is the diffusion path defined in (I.
Finally note that from the data processing inequality, it follows that the KL divergence between the marginals at final time is bounded by KL (f1/x (191/<,U.DALD) I KL (Q I[QU-ALD).
By Cordero-Eneinar Encinar et 2025 Lemmas 3.3 and 4.2], the action A(p) is bounded when the base distribution v is either a Gaussian or Student's 's t distribution, provided that the target distribution æ has finite second order moment and the annealing schedule X satisfies the following assumption.
A5. Let X : R+ + there exists a constant Cx satisfying either of the following conditions
[0,1] be non-decreasing and weakly differentiable, such that if I0)
max (E[0,T]
O,log Aul Ca
or
O.A: VA.(1
max [E[O,T]
Cx2
(28)
A.)
or if v follows Student's I distribution, then condition (28)
holds for some constant Cx-
D Experiment details
The code to reproduce our experiments available at https://// /stps://githhb... gitb.......... com/paulaoak htw.iiinbbbaevrraingg.. git All experiments were implemented using JAX.
D.1 Target distributions
Mixture of Gaussians (MoG) and mixture of Student's t (MoS) distributions The mixture of Gaussians distribution (8-MoG) consists of 8 equally weighted Gaussian distributions with mean my 10 (1 cos + sin(2ri/8)) for shifted the distribution instead of considering the usual benchmark centred at 0 to make sampling more challenging. We note that some of the baselines get stuck in modes close the initial standard Gaussian distribution, unlike our methods.
(9,..., 7} and covariance 0.7 We have
For the mixture of 40 Gaussians in dimensions 2 and 50, the modes are equally weighted, with means sampled uniformly from a hypercube of side length 40, and all the covariance matrices set to the identity matrix J.
The mixture of Student's distributions consists of 10 standard Student's distributions, each with 2 degrees of freedom (t2). Following [Blessing et al.) 2024, Chen et al. 2225), the mean of each component is sampled uniformly from a hypercube with a side length of 10. We evaluate this benchmark in dimensions 10 and 50.
Rings This distribution is defined by the inverse polar reparameterisation of a distribution Ps which has itself a decomposition into two univariate marginals: pr and pe. The radial component p, is a mixture of 4 Gaussian distributions Aii 0.15²) for E (o.....3 3) describing the radial positions. The angular component pe is a uniform distribution over the interval [0, 2x...
Funnel
The density of this distribution is given by ((2) x Nxx10 n7 N(r,0, for 10 with 77 Neal 2003).
(x:)
10
Double well potential (DW) The unnormalised density of the d-dimensional DW is given by π X exp(- (a) 8)2 - Σd Di=+t x²)) with m E and a separation parameter å E oo). This distribution has 2m modes, and larger values of 5 make sampling more challenging due to higher energy barriers. Ground truth samples are obtained using rejection sampling with a Gaussian mixture
proposal distribution Midgley et al. 2023).
41

---

Bayesian logistic regression examples We consider binary classification problems on two bench- mark datasets: lonosphere (dimension 35) and Sonar (dimension 61). Given a training dataset D
where Ij R and Vj, the posterior distribution over the model parameters (w,b), Rd and E R, is defined by
p(w,b(D) x p(w, b) II Bernouli(yy, o( x b)),
where a denotes the sigmoid function. Following Crenioux et al. 2024), we place independent Gaussian priors on the parameters p(w, b) (w;; 1, (b; 2.5").
Statistical physics model: o² distribution The goal is to sample metastable states of the stochastic Allen-Cahn equation 6² model (dimension 100) [Albergo al. 2019
Gabrié et al. 2022 The
o² model is continuous relaxation of the Ising model which is used to study phase transitions in
statistical mechanies. Following [Albergo et al. 2019 Gabrié et al. 2022 Grenioux al. 2024
we consider a version of the model discretised on a -dimensional grd size d 100. Each
configuration of the model is represented by a d-dimensional vector (6. We clip the field to at the boundaries by defining do
0d+1 O. The negative log-density of the distribution is defined as ad 2
In TA (6) = B
(b: Di+1)
2 + ho;
Aad
We fix the parameters a = 0.1 and B 20 to ensure a bimodal regime, and vary the value of We denote by w.y the statistical occurrence of configurations with 0d/2 0 and w- the statistical occurrence configurations with 6a/2 < At = O, the distribution is invariant under the symmetry o -6
expect W+ w.. For h 0. the negative mode becomes dominant.
When the dimension d is large, the relative probabilities of the two modes can be estimated by Laplace approximations at O-th and 2-nd order. Let 24 and ot denote the local maxima of the distribution, these approximations yield respectively
] x det Hi, (d2) (of) x det Hi, (d4) 1/22 where Hμ, is the Hessian of the function In (Ø).
Th (6
(o
For this last benchmark, we only compare our method with SLIPS [Grenioux et al. 2024], as it the only baseline capable of accurately recovering samples from the target distribution while preserving the correct relative mode weights.
D.2 Evaluation metries
To evaluate the quality of the generated samples, we consider the entropy regularised Wasserstein-2 distance (Peyré and Cuturi 2019], with regularisation parameter & 0.05, for distributions where true samples are available. This metric can be efficiently computed in JAX using the OTT library (Cuturi al. 2022).
For the Funnel benchmark, we assess the sample quality using the Kolmogorov-Sminng -iiinov distance. Specifically, we use the sliced version introduced in Grenioux et al. [023 Appendix D.1].
In the Bayesian logistic regression tasks, performance is measured via the mean predictive log-
likelihood, computed as p(w, b(Dies), where Drest is held-out test dataset not used during training.
D.3 Algorithms and hyperparameters
All algorithms are initialised with samples drawn from standard Gaussian distribution to ensure fair comparison. This contrasts with the approach in Grenioux et al. [2024), where the initial samples are drawn using prior knowledge of the target distribution via quantity denoted as R,, which is upper bounded by the scalar variance of the target т. While we observe that initialising from this distribution can improve performance for methods, such information is typically unavailable in real-world scenarios. Nevertheless, for the of model, we consider the same implementation for SLIPS as described in Grenioux et al. [2024], given the complexity of this benchmark.
Besides, to ensure consistency, all algorithms use the same number of energy evaluations, as detailed for each benchmark distribution in Table
42

---

Table 3: Number of energy evaluations for each benchmark
8-MoG
40-MoG
40-MoG
Rings (d= 2)
Funnel
(d=22
(d (d=2) 2)
(d = 50)
(d = 10)
Energy evaluations
3 x 10
3 x 10
5 x 10
3 x 10
5 x 10
5-dim DW
10-dim DW
50-dim DW
Ionosphere
Sonar
(m = 5,8 = 4)
(m = 5,8 =3)
(m = 5,8 2)
(d = 35)
(d = 61)
Energy evaluations
3 x 10
5 x 10
5 x 10
1 x 10
3 x 10
a² model
10-MoS
10-MoS
10-MoS
(d = 100)
(d (d=2] 2)
(d = 10)
(d = 50)
Energy evaluations
1 10
10
10
1 10
Hyperparameter selection Foreach baseline algorithm, we perform grid search over predefined set of hyperparameter values. Selection is based on the corresponding performance metric, computed using 4096 samples. The selected hyperparameters for each algorithm are summarised below.
The SMC and AIS algorithms define a sequence of annealed distributions for k E [0, K]. At each intermediate distribution, we perform n 64 MCMC steps. The parameter RMCMC selected from the interval [20, 100) using a grid step of 4. The total number of distributions K is chosen such that the product
K matches the number of target evaluations specified in Table R
For ULMC algorithm, we tune three hyperparameters: the mass M, the friction coefficient T and the step size h. In all benchmarks, we fix the mass the identity matrix M L. The step size is chosen within the following grid h E {0.000,... 0.009, 0.01..... 0.09, 0.10, 0.11, 0.20). The selected values for the step sizes h are provided Table Additionally, we use non-constant friction coefficient. Specifically, when running the algorithm for L steps, the friction coefficient remains constant at I min mmax = 5 during the remaining steps. The values of I min and I Fmax are chosen within the grid {0.001, 0.005, 0.01, 0.09, 0.1...., 0.9, 1, ..,
for the first L/2 steps and then increases linearly from F min to
For the PT algorithm, the number of chains K is selected from the interval 20) with a grid step of 2. The corresponding temperatures are chosen to be equally spaced on logarithmic scale, with the minimum temperature fixed at 1 and the maximum temperature chosen from the grid {10?, 103, 10º, 10 10 10"3. In our experiments, we use K parallel chains with temperatures [1.0, 5.6, 31.6, 177.8, 1000.0}. Each chain employs Hamiltonian Monte Carlo sampler with n leapfrog steps per sample and step size h specified in Table The grid for the step size is the same as that of ULMC. The number of leapfrog steps n is chosen so that the total number of target evaluations matches the computational budget defined Table 3
Table 4: Selected step sizes for ULMC and PT algorithm across experiments.
8-MoG
40-MoG 2-dim
40-MoG 50-dim
Rings
Funnel
ULMC
0.05
0.05
0.01
0.05
0.03
PT
0.10
0.10
0.03
0.10
0.05
5-dim DW
10-dim DW
50-dim DW
Ionosphere
Sonar
ULMC
0.05
0.05
0.01
0.04
0.04
PT
0.10
0.05
0.05
0.05
0.03
The DiGS algorithm uses K noise levels ranging from ag to 01. The number of noise levels K is selected within the interval [1, 20] with a grid step of 2. The maximum and minimum noise levels, ai and ar, respectively, are chosen from the grid (0.05,... ...,. 0.09, 0.1, 0.2, 1). At each noise level, we perform neaass Gibbs sweeps, each consisting of maala denoising steps using MALA with step size h. The number of Gibbs sweeps nights selected from the interval [50, 500) using a grid step of 50, while the grid for the step size is the same as that of ULMC. The final values of K,
43

---

KKK 01,, Gibbbs and are provided in Table 5 for each experiment. The number of MALA steps, MALAA is determined based on the computational budget.
Table 5: Selected hyperparameters for DiGS algorithm across experiments.
8-MoG
40-MoG 2-dim
40-MeG 50-dim
Rings
Funnel
K
3
1
5
3
aK
0.1
0.1
0.1
0.1
0.1
01
0.5
0.1
0.9
0.5
nobbs
200
300
100
300
200
h
0.05
0.10
0.01
0.10
0.03
5-dím DW
10-dim DW
50-dim DW
lonosphere
Sonar
K
3
5
OKK
0.1
0.1
0.1
0.1
0.1
0.1
0.5
0.8
0.5
0.5
TiGibs
300
200
100
200
200
0.06
0.01
0.005
0.03
0.01
The RDMC algorithm has hyperparameter T, corresponding to the final time of the OU process, its value is provided Table 6 All other parameters follow the implementation detailed in Grenioux et al. 2024).
Table 6: Selected hyperparameters for RDMC algorithm across experiments.
8-MoG
40-MoG 2-dim
40-MoG 50-dim
Rings
Funnel
T
-og((080)
1og(0.75)
log(0.70)
log(0.85)
log (0.90)
5-dim DW
10-dim DW
50-dim DW
Ionosphere
Sonar
T
log(0.70)
log((.70)
log(0.75))
llg((.95)
log(0.95)
For the SLIPS algorithm, we adopt the implementation and hyperparameters provided in the original paper [Grenioux et al. 2024]. To ensure a fair comparison with other algorithms initialised from standard Gaussian distribution, we set the scalar variance parameter R2 Vd. yielding
R2z R 1. While this choice may degrade SLIPS' performance compared to using an optimally tuned R. estimating such parameter in practical scenarios is often non trivial. Nevertheless, for completeness, we report below the performance of SLIPS with its optimal parameters alongside our algorithms. In this setting, the results are comparable. However, key advantage of our method that it does not require estimation of the scalar variance. Moreover, we initialise our algorithms using the same informed choice R, by setting the base distribution (0, I²),, with o R2/v Vd, we observe performance improvement, particularly in the overdamped regime.
For all numerical experiments in the main text, use the underdamped version of our algo- rithms, as it yields improved performance. For MULTALMC. the required hyperparameters include the schedule function A,, the values of As and A, the mass matrix M, E, the fric- tion coefficient T, the step size h, and the number of SROCK steps The grid for the step size is the same as that of ULMC. The scale separation parameter e is chosen within the grid (0.001, 0.005, 0.01, 0.09, 0.10, 0.11..... 0.203. The parameter is chosen from the interval [0.3, 0.7] using grid step of 0.05. In all experiments, we select X5 = 0.01, A = 0.6, M = 1, and Additionally, we consider a time-dependent friction coefficient. Specifically, when running the algorithm for L steps, the friction coefficient remains constant at I miin for the first L/2 steps and then increases linearly from I min I max during the remaining steps. The values of I min and I max are selected within the grid (0.001, 0.005, 0.01...... 0.09, 0.1,.... 0.9. 1. 1.5..... 10). The mass matrix is set to the identity. is worth noting that the number of SROCK steps can be reduced, which may result in small compromise in performance. This effect is further analysed in Appendix F.2 Lastly, the number of iterations determined based on the computational budget specified in Table
44

---

Table Metrics for different benchmarks averaged across 30 seeds. The metric for the mixture of Gaussian (MoG), Rings and the double well potential (DW) is the entropy regularised W distance (with regularisation parameter 0.05), the metric for the Funnel is the sliced Kolmogorov- Smirnov distance and the metric for the Bayesian logistic regression on Tonosphere and Sonar datasets the average predictive posterior log-likelihood on a test dataset. We compare the performance of algorithms (initialised with standard Gaussian distribution) with that of SLIPS with optimal value the parameter Rg refer to as SLIPS (Ra).
8-MoG (L)
40-MoG (L)
40-MoG (L) (d = 50)
Rings ())
Funnel (L) (d = 10)
Algorithm
(d (d=2) 2)
(d = 2)
(d 2)
SLIPS (Rx)
0.66 + 0.11
0.98 + 0.06
20.85 0.63
0.19 0.02
0.029 + 0.007
MULTALMC
0.65 0.07
0.91 0.04
20.58 0.71
0.18 0.02
0.032 + 0.005 0.031 + 0.005
MULICDIFF
0.62 + 0.10
0.93 + 0.06
20.13 0.59
0.19 + 0.03
5-dim DW (L)
10-dim DW (L)
50-dim DW (L)
Tonosphere (1)
Sonar (+) (d = 61)
Algorithm
(m= 5,6 4)
(m = 5,8 = 3)
(m = 5, 2)
(d = 35)
SLIPS (Rtx)
2.05 + 0.30
4.45 + 0.36
14.99 0.59
86.72 0.10
109.38 + 0.12
MULTALMC
2.08 + 0.16
4.48 + 0.40
14.03 0.63
86.85 0.09
109.00 + 0.13 -109.60 0.21
MULTCDIFF
1.95 1 0.24
4.23 1 0.37
13.98 0.56
86.33 0.10
Table 8: Selected hyperparameters for MULTALMC algorithm across experiments.
8-MoG
40-MoG 2-dim
40-MoG 50-dim
Rings
Funnel
Xe
Linear
Lincar
Linear
Linear
Cosine-like
0.10
0.05
0.05
0.10
0.07 0.01, 0.5 0.003
I min, I max
0.07,
0.01, 0.5
0.01,0.5
0.1, 0.005
h
0.005
0.005
0.001
5-dim DW
10-dim DW
50-dim DW
Ionosphere
Sonar
Cosine-like
Cosine-like
Cosine-like
Cosine-like 0.10
Cosine-like
E
0.10
0.07
0.05
0.10
T min, T max
0.01, 0.5
0.01, 0.5
0.01, 0.5 0.001
0.5 0.005
0.1,0.5 0.5
h
0.005
0.001
0.005
8 model
10-MoS 2-dim
10-MoS 10-dim
10-MoS 50-dim
X
Cosine-like 0.05
Linear
Linear
Linear
0.10
0.10
0.05
mins max
0.01,
0.1,0.5 0.5
0.1,0.5
0.1, 0.5
h
0.001
0.005
0.005
0.001
On the other hand, for MULTCDIFF, the following hyperparameters need to be specified: the friction coefficient I, the scale separation parameter E, the final time T of the OU process, the time discretisation 0 To
< Tú. T, and the number of SROCK steps s. The values of T, and s are set to be the same as in the MULTALMC algorithm. The time discretisation is chosen such that the difference Äri Xr, is constant, where denotes the OU schedule. The number of iterations L is determined based on the computational budget for each benchmark.
D.4 Computation time
All experiments were conducted on a GPU server consisting of eight Nvidia GeForce RTX 3090 Ti GPU cards, 896 GB of memory and 14TB of local on-server data storage. Each GPU has 10496 cores as well as 24 GB of memory.
Recall that to ensure a fair comparison, we fixed the number of energy evaluations across all algorithms. As result, the runtime of ULMC, DiGS, RDMC, SLIPS, and our proposed methods, MULTALMC and MULTCDIFF, are broadly similar. In contrast, algorithms such as SMC, AIS, and PT exhibit longer runtimes due their accept/reject steps. This is demonstrated in Figure 2
45

---

which shows average computation times (over 30 random seeds) for each method on the 40-MoG benchmark in 50 dimensions.
800
750
Computation Time (s)
700
650
600
550
500
450
SMC
AIS
ULMC
PT
DIGS
RDMC
SLIPS
MultALMC
Multcoift
Figure 2: Boxplots of computation times for different algorithms on the 40-MoG benchmark in 50
dimensions, averaged over 30 random seeds.
In addition, we analyse how the runtime of our methods scales with dimensionality. To that end, Table(Jincludes performance metrics, number of energy evaluations, and computation times for 40-MoG benchmark across different dimensions, averaged over 30 random seeds. The performance metric used the entropy-regularised Wasserstein-2 - distance, W., with regularisation parameter 0.05.
Table Number of energy evaluations, entropy regularised Wasserstein-2 distance (with regularisa- tion parameter 0.05) and computations times for our sampling methods, MULTALMC (Algorithm) 3 and MULTCDIFF (Algorithm 4 evaluated on the 40-MoG benchmark across different dimensions. The results are averaged over 30 random seeds.
d = 10
d = 20
d = 50
Energy evaluations
x 10
6 x 10
1 x 10
2 x 10
x 10
MULTALMC
0.91
0.04
1.64 1 0.12
4.10 1 0.27
8.76 1 0.34
20.58 1 0.71
W
MULTCDIFF
0.93 + 0.06
1.59 1: 0.20
3.92 1 0.36
8.68 1 0.41
20.13 4: 0.59
MULTALMC
35 + 5
57 9
111 + 13
223 - 28
517 + 28
Time (s)
MULTCDIFF
33
58 77
104 + 15
231 21
502 + 30
E
Additional numerical experiments
E.1
Comparison of overdamped and underdamped dynamics
As noted in the main text, the overdamped version of MULTALMC requires small value of K, which corresponds slowly varying dynamics driven by V log P. perform well in practice. However, this leads to a large number of discretisation steps since we use small step size with respect to this slowly changing dynamics resulting high computational costs. In contrast, the underdamped version achieves better performance without the need for such small step sizes, thanks to the faster convergence properties of underdamped dynamics [Eberle and Lörler 2024). To further analyse this, we compare in Table T0 the performance and number of energy evaluations of the overdamped and underdamped versions of MULTALMC, as specified in Algorithms 2 and respectively, on the 40-MoG benchmark across different dimensions. The results show that the overdamped version
46

---

requires approximately an order of magnitude more energy evaluations (and hence time steps) to achieve performance comparable to the underdamped version.
Table 10: Performance metrics on the 40-MoG benchmark across varying dimensions, averaged over 30 runs with different random seeds, along with the number of energy evaluations for the overdamped and underdamped versions of MULTALMC, as defined in Algorithms 2 and respectively.
MULTALMC
d 2
d 5
d 10
d 20
50
W
1.12 - 0.10
1.88 0.22
4.46 0.39
9.23 0.67
21.33 = 0.98 1 x 10
Overdamped
# evaluations
1 x 10
8 x 10
x 10"
6 x 10"
W # evaluations
0.91 +: 0.04
1.64 + 0.12
4.10 + 0.27
8.76 + 0.34
20.58 + 0.71
Underdamped
3 x 10
x 10
x 10
x 10
x 10
E.2
Analysis of the impact of the number of SROCK steps in the discretisation
The number SROCK steps used our Poooittnss-MULTAAAAN
and MULTCDIFF imple- mented in Algorithms 2 3 and 4 plays a critical role in balancing computational efficiency and numerical performance. Increasing the number inner SROCK steps generally leads to better performance [Abdulle et al., 2018); however, this comes at the cost of greater computational overhead. Ideally, we aim to identify the minimum number of steps required to maintain strong performance and computational efficiency.
In the numerical experiments presented in Section 5 we use five SROCK steps. We further analyse the sensitivity of our accelerated methods, MULTALMC (Algorithm(3 and MULTCDIFF (Algorithm + the number of SROCK steps. To this end, we report the entropy regularised Wasserstein-2 distance for the 40-MoG benchmark in 50 dimensions (Table [1]] when running our methods with varying numbers of SROCK steps while keeping the number of sampling steps for the slow process fixed.
Table 11:
Performance on the 40-MoG benchmark in 50 dimensions using varying numbers of SROCK steps in the discretisation of MULTALMC and MULTCDIFF, as proposed in Algorithms B and respectively. The number of sampling steps for the slow process held fixed. The standard deviation is computed over 30 runs with different random seeds. The regularisation parameter for the entropy regularised Wasserstein-2 distance W is set
0.05.
SROCK steps
3
7
11
MULTALMC
21.11 0.74
20.58 + 0.71
20.29 + 0.68
20.08 + 0.50
19.99 = 0.45
w
MULTCDIFF
20.92 0.63
20.13 = 0.59
20.01 0.57
19.96
0.48
19.89 0.40
E.3 Stability across hyperparameter values
Our algorithms demonstrate robustness and stability across a broad range of hyperparameter values, as shown in Figures and 4 These results suggest that while some tuning is beneficial, our methods do not require highly sensitive or exhaustive hyperparameter optimisation.
F Limitations and future work
We elaborate here on the limitations of our method discussed in Section O
First, the theoretical guarantees presented in Section relyy on stringent assumptions, such as strong convexity of the potential V, and Lipschitz continuity of the score function of the target V log Relaxing these assumptions, e.g., to strong convexity outside compact region or to satisfying weak functional inequalities, could broaden the applicability of our analysis. The former would accommodate target distributions that are multimodal within a compact region and Gaussian-like in the tails, while the latter could capture heavy-tailed distributions.
47

---

24
30
MuitALMC
MuRtALMC MulcCDiff
MultCDiff
28
22
26
21
20
20
19.
0.02 0.04
0.06
0.08 0.10 0.12 0.14
10
Number of steps
(a) Performance against scale separation parameter s.
(b) Performance against number of steps.
Figure 3: Ablation results on the 40-component Mixture of Gaussians benchmark in 50 dimensions.
Results are averaged over 30 random seeds.
Avg Predictive Log-Liken LigLeikelhood
Avg Predictive Log-Likee LikLiikelioodd
-110
-112
-108
116
110
-118
-120
MultCDiff
MultALMC
-111-
MukCDiff
MURALMC
0.02
0.04
0.06 0.08
0.10
0.12
0.14
109 Number of steps
(a) Performance against scale separation parameter E.
(h) Performance against number of steps.
Figure 4: Ablation results the Sonar benchmark. Results averaged over 30 random seeds. The performance metric the average predictive posterior log-likelihood on test dataset.
Additionally, the current method requires manual tuning hyperparameters such as step size 8, scale separation parameters æ and friction coefficient I. Automating this tuming process, similar to the approach Bllsssing al. 20255. valuable direction for improving usability and robustness.
Further research could explore extending the controlled diffusion framework to heavy-tailed target distributions, which pose additional challenges explained in Appendix [..22) or developing more efficient numerical schemes for implementing the proposed multiscale samplers.
G Broader impact
This work introduces MULTALMC and MULTCDIFE, two training-free multiscale diffusion sam- plers that enable efficient, provably accurate sampling from complex distributions. These methods have the potential to significantly reduce the computational and environmental footprint of generative modelling and Bayesian inference. Potential societal benefits include faster scientific discovery, improved uncertainty quantification, and broader participation in high-impact ML research. As with all general-purpon algorithms, misuse is possible, e.g. lowering the cost of harmful synthetic content generation or accelcrating harmful molecule design.
48


Reviews:
[{"review_id": "MS9CvRr7uC", "reviewer": "Reviewer_JLNp", "summary": "This paper considers the problem of sampling using diffusion processes, given oracle access to the gradients of the potential function of the underlying distribution. Their main contribution is a new algorithm based on a multi-scale version of the standard annealed Langevin dynamics. They show that this algorithm provably converges to the true distribution asymptotically, given that the true distribution is strongly log-concave. They complement these theoretical results with some preliminary evaluations both for the Gaussian setting, as well as for the heavy-tailed setting.", "questions": "My main question, as described above, is regarding the motivation of the paper. While the authors briefly justify why this method could be preferred to the other works, these concerns do not really hold for the strongly log-concave setting that their theoretical results hold for. If the authors could justify why the results are interesting in comparison to prior work specifically in this strongly log-concave setting, that would be very helpful.", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Quality: From a technical point of view, the results in this paper are quite nontrivial and above the bar for acceptance.\n\nClarity: The results in the paper are presented clearly, at least within the main text of the paper.\n\nSignificance: This is my main concern for the paper. Namely, it is by now fairly well understood how to sample from strongly log-concave distributions given this sort of oracle access, e.g. the papers that the authors cite of Chen-Chewi-Salim-Wibisono and Lee-Shen-Tian. Moreover, these papers yield non-asymptotic rates, as opposed to the asymptotic convergence demonstrated here. Despite the paper’s claims, I think it is unlikely that this approach will work for complex real-world settings (as it doesn’t learn a score function). Indeed, the evaluations in the paper are on simple, synthetic settings. Thus, I think the main motivation must be theoretical, but the results are then largely superceded by these previous works. \n\nOriginality: The method is original, to my understanding.", "was_revised": true, "final_justification": "After reading the rebuttal and in discussion with the other reviewers, I am okay with increasing my score.\n\nTo clarify, and this was definitely because I was being hasty, I was generally referring to the fact that good, efficient rates are known for sampling from strongly log-concave and smooth distributions under log-grad access, but indeed, the papers I referenced are not quite the right references for this. Regardless, as the authors point out, the novelty of the paper must be for these more complex, empirical settings.\n\nI am not an expert on the state-of-the-art in this regard, but after discussion it seems like these benchmarks are acceptable, albeit not ideal, for this area, which is surprising to me. For this reason, though, I am okay with raising the score.", "rebuttal": "Thank you for your review and insightful comments. We also sincerely appreciate your positive remarks on the quality and clarity of our work.\n\nWe reply to your concerns and questions below, please let us know if anything is unclear.\n\n> This is my main concern for the paper. Namely, it is by now fairly well understood how to sample from strongly log-concave distributions given this sort of oracle access, e.g. the papers that the authors cite of Chen-Chewi-Salim-Wibisono and Lee-Shen-Tian. Moreover, these papers yield non-asymptotic rates, as opposed to the asymptotic convergence demonstrated here.\n\nWe are not entirely sure whether, in your comment, oracle refers to access to the gradients of the target potential, i.e. $\\nabla \\log \\pi = \\nabla V_\\pi$, or instead to what Chen–Chewi–Salim–Wibisono and Lee–Shen–Tian define as a restricted Gaussian oracle (RGO). For completeness, we address both interpretations below.\n\n**1. Oracle = access to $\\nabla\\log\\pi=\\nabla V_\\pi$**: Our algorithms are designed for sampling tasks, where the goal is to obtain samples from a complex target distribution $\\pi$. In this setting, it is assumed in the literature that we have access to $\\log \\pi(x)$ (up to an additive constant) and its gradient $\\nabla \\log \\pi(x) = \\nabla V_\\pi(x)$. This contrasts with typical generative modelling tasks where the objective is to generate new samples from an unknown data distribution, given only access to a dataset of existing samples. To make this distinction clearer, we provide the following table highlighting the key differences between the two tasks.\n\n|                       | **Sampling**                                                                                  | **Generative Modelling**                                                                                               |\n|-----------------------|------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------|\n| **Available information**     | $\\log \\pi(x)$ (up to a constant) and $\\nabla \\log \\pi(x) = \\nabla V_\\pi(x)$                   | Samples/dataset from underlying data distribution $\\pi_{\\text{data}}(x)$           |\n| **Some typical algorithms**| Markov chain Monte Carlo, sequential Monte Carlo | VAEs, score matching, GANs\n\nHaving said this, when $\\pi$ is strongly log-concave, standard Langevin Monte Carlo (LMC) algorithms can be readily applied, with known convergence guarantees in this setting [1], without requiring the use of the more sophisticated frameworks of Chen–Chewi–Salim–Wibisono or Lee–Shen–Tian, which rely on access to an RGO. However, we would like to emphasise that, when developing our algorithms, our main focus is on complex multimodal distributions, where plain LMC and related methods typically struggle. To address this challenge, we employ methods based on annealed Langevin diffusions, which are better suited for these harder sampling problems.\n\n**2. Oracle as in Chen-Chewi-Salim-Wibisono and Lee-Shen-Tian**: We would like to clarify that our setting does not assume access to an oracle in this sense. The works of Chen-Chewi-Salim-Wibisono and Lee-Shen-Tian rely on access to what they refer to as a restricted Gaussian oracle (RGO), which corresponds in our notation to exact sampling from the conditional distribution $Y\\sim\\rho_{t,X}(y)$ (or $\\pi^{X|Y}$ in their notation). Some of their theoretical results explicitly assume this oracle access and therefore do not account for the complexity of sampling from $Y\\sim\\rho_{t,X}(y)$, which in our case corresponds to the fast process. Our theoretical analysis does include the complexity of this sampling step. \n\nBesides, while both aforementioned works also provide theoretical results when implementing the oracle approximately, via a rejection sampling scheme (Corollary 6 in Chen-Chewi-Salim-Wibisono and Corollary 1 in Lee-Shen-Tian), these approaches assume access to the minimiser of the target potential—an assumption that can be impractical in higher-dimensional settings.\n\n> Despite the paper’s claims, I think it is unlikely that this approach will work for complex real-world settings (as it doesn’t learn a score function). Indeed, the evaluations in the paper are on simple, synthetic settings. Thus, I think the main motivation must be theoretical, but the results are then largely superceded by these previous works.\n\nWe would like to highlight that we do estimate the drift term of the averaged process $\\bar X_t$, Equation (3), (which corresponds to the score function of the true moving target distribution $\\mu_t$). This estimation is made possible via the auxiliary fast process $Y_t$, which enables us to approximate the drift of the marginal $X_t$ without direct score learning. \n\nBesides, our algorithms demonstrate strong performance across a variety of established benchmarks in the sampling literature that go beyond the strongly log-concave setting (some of them arising from real-world applications such as Bayesian inference). In addition, we include a real-world example given by a statistical physics model which has a dimension of 100 (see Figure 1a). Importantly, this example is not only high-dimensional but also inherently challenging due to the free-energy barrier between the modes [2].\n\n> While the authors briefly justify why this method could be preferred to the other works, these concerns do not really hold for the strongly log-concave setting that their theoretical results hold for. If the authors could justify why the results are interesting in comparison to prior work specifically in this strongly log-concave setting, that would be very helpful.\n\nAs we acknowlegde in the discussion section, the strongly-convex assumption required for deriving theoretical guarantees is indeed restrictive. However, the motivation for providing partial theoretical guarantees was to demonstrate that our method is principled and to provide an initial template for more general theory. This does not mean that the proposed methods are restricted to the assumption used in the theoretical analysis. Developing a full theoretical framework for multiscale systems under weaker assumptions remains technically challenging and constitutes a substantial line of research in its own right. For this reason, it falls outside the scope of the present work. Nevertheless, we hope that the connections we establish with the averaging literature will encourage applied mathematicians working on multiscale dynamics to engage with this problem. We are genuinely excited and eager to see what theoretical developments our initial approach may inspire in future work.\n\n**Conclusion**\n\nWe sincerely appreciate your recognition that our work is “quite nontrivial and above the bar for acceptance” (as noted in your Quality comment). We also fully acknowledge your concerns regarding the scope of our theoretical results, which are currently confined to the strongly log-concave and asymptotic settings. \n\nThat said, we would like to emphasise that while our proof-of-concept theoretical guarantees focus on the strongly log-concave case, our method is explicitly designed for — and evaluated on — real-world scenarios that are non log-concave and high-dimensional. To this end, we **have benchmarked our algorithms against a broad set of state-of-the-art methods on challenging, high-dimensional sampling problems** demonstrating both practicality and competitive performance.\n\nIn light of this, we kindly ask you to consider reevaluating your score. We believe that current results in our manuscript provide a comprehensive evaluation of our samplers, including a practical (and highly challenging) statistical physics model. \n\nWe look forward to engaging with you further during the rebuttal period to address any remaining concerns and welcome your feedback on our responses. Thanks again!\n\n[1] Durmus, A. et al. (2019) Analysis of Langevin Monte Carlo via Convex Optimization. Journal of Machine Learning Research, 20 (73).\n\n[2] Gabrié, M. et al. (2022) Proc. Natl. Acad. Sci. U.S.A. 119 (10)."}, {"review_id": "SrAtfoaAKb", "reviewer": "Reviewer_HLYz", "summary": "This paper proposes a multi scale dynamics-based approach to improve the sampling efficiency of score-based samplers.  Existing method typically rely on nested iterative MCMC frameworks to sample from an unnormalized distribution without access to samples.  To alleviate the complexity of these approaches this paper proposed a slow-fast dynamics based on different SDEs to estimate the score function and different points along a diffusion path without resorting to averaging through MCMC-generated samples as in existing approaches.  Results on several simple examples shows the superior performance of the proposed method.", "questions": "1. How does the empirical performance change with the dimensionality of the underlying signal? Do the number of steps and scale separation parameter have to be tuned very carefully as the dimensionality changes?\n2. Is there any intuition on how to choose the hyper parameters of the proposed algorithm carefully without too much experimental tuning?  What happens when the conditional distribution mixing time is set too slow or too fast?", "limitations": "Yes", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n\n1. The goal of the paper and the proposed method is clear.  The background and related work is thorough and well written.\n2. The proposed method with parallel slow-fast dynamics is fairly novel and addresses the computational burden of existing approaches.\n3. The drawbacks of each proposed scheme are explained and it helps motivate additional improvements with subsequently proposed algorithms\n4. The experimental evaluation is strong and show the competitive performance of the proposed scheme.\n\n**Weaknesses**\n\n1. The practical viability of the approach to some real-world settings is not immediately clear given the constraints and assumptions placed within the theoretical analysis.", "was_revised": true, "final_justification": "I believe that the paper is well written and the analysis is thorough.  Despite not being evaluate on very high dimensional datasets, the experiments are fairly extensive and provide evidence for the viability of the proposed method.", "rebuttal": "Thanks a lot for your thoughtful review and helpful comments. We also greatly appreciate your positive remarks on our work.\n\n> The practical viability of the approach to some real-world settings is not immediately clear given the constraints and assumptions placed within the theoretical analysis.\n\nOur algorithm is specifically designed to sample from *diffusion paths* with the aim of targeting highly complex probability distributions in practice. In our experiments, we demonstrate that the proposed algorithms exhibit strong empirical performance in general settings beyond strong convexity. This includes a real-world application in statistical physics involving a 100-dimensional system (see Figure 1a), as well as benchmarks from Bayesian inference tasks such as Ionosphere and Sonar.\n\nWe acknowledge that the assumptions made in the theoretical section are arguably strong. Our objective in providing this theoretical guarantee was to demonstrate that the method is principled under suitable conditions. Moreover, the analysis serves as a template that can be extended to more general assumptions in future work.\n\nAnalysing multiscale systems under weaker assumptions presents significant technical challenges and constitutes a substantial line of research in its own right. For this reason, it falls outside the scope of the present work. Nevertheless, we hope that the connections we establish with the averaging literature will encourage applied mathematicians working on multiscale dynamics to engage with this problem. We are genuinely excited and eager to see what theoretical developments our initial approach may inspire in future work.\n\n> How does the empirical performance change with the dimensionality of the underlying signal? Do the number of steps and scale separation parameter have to be tuned very carefully as the dimensionality changes?\n\nAs the dimension of the underlying signal (i.e., the $(Y_t)\\_{t\\geq 0}$ process), and consequently the dimension of the $(X_t)_{t\\geq 0}$ process, increase, convergence to the target distribution becomes slower, requiring a larger number of time steps. However, convergence still occurs. This is illustrated in Table 9 (Appendix D.4), where we report performance metrics (specifically, the entropy-regularised Wasserstein-2 distance) alongside the number of energy evaluations (which is in one-to-one correspondence with the number of steps) and computations times for our sampling methods on the 40-MoG benchmark across different dimensions.\n\nOn the other hand, we observed that **our algorithms demonstrate robustness and stability across a broad range of hyperparameter values**, as shown in the table below. We plan to include a plot (extending the results of this table) on a finer hyperparameter grid in an updated version of the paper. These results suggest that while some tuning is beneficial, our methods do not require highly sensitive or exhaustive hyperparameter optimisation.\n\nTable 1. Ablation results on the 40-component Mixture of Gaussians benchmark in 50 dimensions for the scale separation parameter $\\varepsilon$ and the number of steps. Results are averaged over 30 random seeds. Notably, the performance remains stable across a wide range of hyperparameter values.\n|        |    $\\mathbf{\\varepsilon}$   | 0.01   | 0.03 | 0.05 | 0.07   | 0.10      | 0.12|0.15|\n|----------------|---------|-------------|---------------|---------|-------------|-------------|-------------|-------------|\n| **MultALMC**          | $W_2^\\delta$  |   25.48±1.32    |    20.99±0.83    |   20.58±0.71|  21.68±0.92|  24.50±1.17| 26.89±1.80| 28.43±2.57|\n| **MultCDiff**             | $W_2^\\delta$    |    24.33±1.41     |     20.72±0.58     |    20.13±0.59 |    21.02±0.70     |   23.81±1.23| 26.21±1.66| 28.01±2.42|\n\n|        |    Number of steps  | $\\mathbf{5\\times 10^5}$   | $\\mathbf{7.5\\times 10^5}$ | $\\mathbf{1\\times 10^6}$ | $\\mathbf{2.5\\times 10^6}$   | $\\mathbf{5\\times 10^6}$ |\n|----------------|---------|-------------|---------------|---------|-------------|-------------|\n| **MultALMC**          | $W_2^\\delta$  |   22.45±1.22    |    21.83±0.95    |   20.47±0.72| 20.15±0.50|  19.40±0.36| \n| **MultCDiff**             | $W_2^\\delta$    |    22.77±1.04    |     21.16±0.76     |   20.08±0.53 |    19.69±0.41     |   19.12±0.37| \n\n\nTable 2. Ablation results on the Sonar benchmark for the scale separation parameter $\\varepsilon$ and the number of steps. Results are averaged over 30 random seeds. The performance metric is the average predictive posterior log-likelihood on a test dataset. Notably, the performance remains stable across a wide range of hyperparameter values.\n|        |    $\\mathbf{\\varepsilon}$   | 0.01|0.03 | 0.05 | 0.07   | 0.10      | 0.12| 0.15|\n|----------------|---------|-------------|---------------|---------|-------------|-------------|-------------|-------------|\n| **MultALMC**          | Performance  |   -119.84±1.41    |     -116.16±1.07    |   -112.91±0.83|  -110.52±0.49|  -109.05±0.13| -109.84±0.52| -113.14±1.06|\n| **MultCDiff**             | Performance    |    -118.73±1.67     |     -115.83±1.39     |    -112.20±0.97 |    -110.11±0.51     |   -109.60±0.21| -109.78±0.55| -112.96±1.23|\n\n|        |   Number of steps   |  $\\mathbf{5\\times 10^5}$   | $\\mathbf{7.5\\times 10^5}$ | $\\mathbf{1\\times 10^6}$ | $\\mathbf{2.5\\times 10^6}$   | $\\mathbf{5\\times 10^6}$ | \n|----------------|---------|-------------|---------------|---------|-------------|-------------|\n| **MultALMC**          | Performance  |  -111.14±0.52    |     -109.06±0.15    |   -108.73±0.20|  -107.22±0.13|  -106.49±0.14| \n| **MultCDiff**             | Performance    |    -110.86±0.58     |     -109.55±0.21     |   -109.02±0.19 |    -107.61±0.16     |   -106.37±0.09| \n\n\n> Is there any intuition on how to choose the hyper parameters of the proposed algorithm carefully without too much experimental tuning? What happens when the conditional distribution mixing time is set too slow or too fast?\n\nThe main hyperparameters in our algorithms are the step size and scale separation parameter ($\\varepsilon$). The intuition for choosing the step size is similar to that used in classical Langevin-based algorithms (see [1,2]).In particular, this parameter typically scales inversely with the dimension of the problem. Regarding the scale separation parameter, we empirically observe that values in the range 0.01-0.11 yield good performance across a wide range of problem settings (see Table 8, Appendix D.3, where we report the optimal value of $\\varepsilon$ for each benchmark distribution).\n\nIf the conditional distribution's mixing time is set too slow, then the obtained samples from the process $Y_t$ will not correspond to that of the correct conditional distribution. This results in the drift term of the slow process ($X_t$) being poorly estimated. That is, the drift term will not accurately reflect the drift term corresponding to the true moving target, which in turn fails to capture the correct dynamics of the underlying moving target. Intuitively, this leads to a lag in tracking the moving target distribution, which may cause some modes of the final target distribution to be missed. \n\nOn the other hand, if the conditional distribution mixing time is too fast, the particles from the slow process $(X_t)$ may exhibit overly noisy or \"wiggly\" trajectories, hindering convergence over time.\n\nWe hope that the additional discussions and ablation studies have addressed the reviewer’s questions, and we look forward to receiving further feedback!\n\n[1] Chewi, S. et al. (2024) Analysis of Langevin Monte Carlo from Poincare to Log-Sobolev. Foundations of Computational Mathematics.\n\n[2] Milstein, G. N. (1994) Numerical integration of stochastic differential equations. Springer."}, {"review_id": "1TxJGZk5pV", "reviewer": "Reviewer_MkbF", "summary": "This paper proposed a novel multiscale sampling method for score estimation and develop two algorithms based on annealed Langevin dynamics and controlled diffusions, respectively. This paper validate the proposed method on multiple synthetic data with various dimensions. This paper also provide convergence guarantees.", "questions": "All of the questions are above.", "limitations": "Can this multiscale score estimation method scale up to real-world data.", "rating": 4, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n- Score estimation with multiscale sampling approach is interesting and novel. The proposed method use the multiscale sampling approach to improve the sampleing efficiency and quality.\n\n- This paper provide theoretical guarantees for proposed method.\n- The proposed method achieve promising results on low-dimension data.\n\n\n\nWeaknesses:\n\n- I'd like to know the results on real-world data, at leaset like mnist or cifar10 because I find the performance on high-dimensional and complex data (table 2) doesn't has consistent improvement.\n- Lack of comparison of computational cost.", "was_revised": true, "final_justification": "I'll keep my score.", "rebuttal": "Thank you for your thoughtful review and insightful comments. We also greatly appreciate your positive feedback on our work.\n\nWe reply to your concerns below, please let us know if anything is unclear.\n\n> I'd like to know the results on real-world data, at leaset like mnist or cifar10 because I find the performance on high-dimensional and complex data (table 2) doesn't has consistent improvement.\n\nWe would like to clarify that our algorithms are designed for **sampling tasks** where we have direct access to a complex unnormalised target distribution and its gradient, and the goal is to generate samples from it. This contrasts with **typical generative modelling tasks** where the objective is to generate new samples from an unknown data distribution, given only access to a dataset of existing samples. To make this distinction clearer, we provide the following table highlighting the key differences between the two tasks.\n\n|                       | **Sampling**                                                                                  | **Generative Modelling**                                                                                               |\n|-----------------------|------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------|\n| **Available information**     | $\\log \\pi(x)$ (up to a constant) and $\\nabla \\log \\pi(x) = \\nabla V_\\pi(x)$                   | Samples/dataset from underlying data distribution $\\pi_{\\text{data}}(x)$           |\n| **Some typical algorithms**| Markov chain Monte Carlo, sequential Monte Carlo | VAEs, score matching, GANs                      |\n\nTherefore, direct sampling from MNIST or CIFAR10, given access to the dataset of images, does not fit naturally into our framework. To apply our methods to these datasets, one would first need to train an energy-based model [1,2] to learn the unnormalised target distribution corresponding to the dataset and then use our proposed method to draw new samples from the learnt distribution. However, in such a setup, the quality of the generated samples would depend not only on the performance of our sampler but also on how well the energy-based model captures the true data distribution. This would make it difficult to isolate and evaluate the effectiveness of the sampling algorithm itself.\n\nThat said, we have evaluated our methods on a wide range of established benchmarks for sampling algorithms [3,4,5] (some of them arising from real-world applications such as Bayesian inference) and included a real-world example involving a statistical physics model which has a dimension of 100 (see Figure 1a).\n\n> Lack of comparison of computational cost.\n\nTo **ensure a fair comparison across different algorithms we use the same number of energy evaluations**, as detailed in Table 3 (Appendix D.3). Thus, the computational cost is kept the same for all compared methods. \n\nAdditionally, in Appendix D.4, we report the computational resources used, along with the computation times (averaged over 30 random seeds) for each algorithm on the 40-component mixture of Gaussians benchmark in 50 dimensions (see Figure 2).\n\nTo further analyse the computational cost, we provide below a comparison in terms of energy evaluations and computation times required by each algorithm on the 40-MoG benchmark (dimension 50) to achieve comparable performance, as measured by the entropy-regularised Wasserstein-2 distance. Results are averaged over 30 random seeds.\n\nTable 1. Comparison of the computational cost across methods on the 40-MoG benchmark in 50-dimensions. Results are averaged over 30 random seeds.\n| Algorithm | $\\mathbf{W_2^\\delta}$       |   Time (s)   |  Energy evaluations |\n|-----------|------------------|------------------|----------|\n| **SMC** |  25.80±2.03   |  1402±96    |   $1\\times10^7$    | \n| **AIS**|  27.14±2.26     | 1327±72     |     $1\\times10^7$   | \n| **PT**   | 23.36±0.39      | 593±56    |  $4\\times10^6$  |\n| **DiGS**  | 23.68±0.91     | 438±31    |    $4\\times10^6$    |\n| **RDMC**  | 24.19±0.82     | 856±40     |  $8\\times10^6$  |\n| **SLIPS**  | 23.71±0.65      | 509±27    |  $5\\times10^6$  |\n| **MultALMC**  | 23.40±0.45      | 379±35     |  $3\\times10^6$  |\n| **MultCDiff**  | **23.31±0.42**     |**363±21**     |  $3\\times10^6$  |\n\n\nWe note that our proposed algorithms consistently require fewer energy evaluations to achieve comparable performance than alternative methods. We plan to include this discussion in an updated version of the paper.\n\nWe hope that this addresses your concerns, and we sincerely request any further feedback that you may be able to provide us during the rebuttal period. Thank you for your continued feedback!\n\n[1] Du, Y. et al. (2021) Improved Contrastive Divergence Training of Energy-Based Model. In Proceedings of the 38th International Conference on Machine Learning.\n\n[2] Song, Y. and Kingma, D. (2021) How to Train Your Energy-Based Models. arXiv:2101.03288.\n\n[3] Grenioux, L. et al. (2024) Stochastic Localization via Iterative Posterior Sampling. In Proceedings of the 41st International Conference on Machine Learning.\n\n[4] Vargas, F. et al. (2024) Transport meets Variational Inference: Controlled Monte Carlo Diffusions. The Twelfth International Conference on Learning Representations.\n\n[5] Richter, L. and Berner, J. (2024) Improved sampling via learned diffusions. The Twelfth International Conference on Learning Representations."}, {"review_id": "mfxpLZLR3S", "reviewer": "Reviewer_RNkH", "summary": "This work introduces learning-free methods for the sampling of complex unnormalized distributions. The methods aim at estimating the score of the target distribution and are based on the insight that the path of distributions associated to a reverse diffusion process can be expressed as a convolution of the target distribution with a noising distribution. It follows that the score of these distributions is an expectation over a distribution that depends on the noising distribution and the target distribution. The methods rely on the concept of introducing SDEs that represent stochastic processes at different time scales. Simulating these coupled SDEs corresponds to carrying out the aforementioned score estimation (fast process) and the corresponding Langevin dynamics (slow process) which generates samples from the target distribution.", "questions": "1. As pointed out in the paper, the multiscale approach avoids the need for a nested algorithm. However, a nested algorithm that does not introduce SDE with different time scales is a natural competitor for the introduced multiscale methods. In such a nested approach, samples ~ $\\rho$ in (2) would be drawn in the inner-loop by simulating the corresponding SDE. While some of the methods in the experimental section rely on such a nested approach it is not clear whether any of those corresponds to a variant of the introduced methods where the multiscale SDE simulation is replaced with a nested SDE simulation.\n\n2. How is it ensured that hyperparameter optimization was carried out with comparable computational resources for the compared sampling methods? What were the hyperparameter grids used for each method and benchmark?", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. Employing multiscale SDEs for sampling a diffusion path appears to be novel and promising.\n2. Well written, clear notation.\n3. Experimental evaluation is extensive and the results indicate that the new methods perform well.\n\n\nWeaknesses:\n\n1. Missing ablations on the heavy-tail distribution in for MULTALMC: the heavy-tail and Gaussian versions of MULTALMC are only compared on the MoS data set but not on any other benchmarks.\n2. The theoretical result relies on log-concavity of the energy function which makes it practically hardly relevant. This problem is, however, openly pointed out by the authors.", "was_revised": true, "final_justification": "The paper represents a valuable contribution to the community. While the presented theoretical contribution may be of limited practical relevance, the strong methodological section and thorough experimental evaluation elevate this work above the bar for publication. It matches all requirements for a publication, and the authors' response to my questions clarified all of my potential concerns.\nConcretely, the authors addressed a question concerning a nested competitor method. Their theoretical argumentation and the newly added experimental results provide in total a satisfying response.\nQuestions about the methodology underlying the hyperparameter tuning and the allocation of computational resources in the experimental evaluation were also answered appropriately. With the additional details provided on the hyperparameter tuning, this work should be well reproducible and it appears to be based on a sound experimental evaluation.", "rebuttal": "Thank you for your thorough review and helpful comments. We also greatly appreciate your positive feedback.\n\nWe reply to your questions below, please let us know if anything is unclear.\n\n> As pointed out in the paper, the multiscale approach avoids the need for a nested algorithm. However, a nested algorithm that does not introduce SDE with different time scales is a natural competitor for the introduced multiscale methods. In such a nested approach, samples ~ $\\rho$ in (2) would be drawn in the inner-loop by simulating the corresponding SDE. While some of the methods in the experimental section rely on such a nested approach it is not clear whether any of those corresponds to a variant of the introduced methods where the multiscale SDE simulation is replaced with a nested SDE simulation.\n\nReverse Diffusion Monte Carlo (RDMC) proposed by Huang et al. (2024) [1] corresponds to the framework that you describe when the averaged process is given by the reverse diffusion process (that is, the setting described in Section 3.2 of our work). We include comparisons of our algorithms with RDMC in Tables 1 and 2. \n\nHowever, RDMC does not explore the cases where the averaged process is given by the time inhomogeneous Langevin diffusion (i.e., Equation (3) of our work) either for Gaussian diffusion paths, where the base distribution $\\nu$ is a Normal distribution, or for heavy-tailed diffusions with a Student’s t base distribution. \n\nThe nested inner-outer loop algorithm that you suggest based on annealed (time-inhomogeneous) Langevin diffusion with a Gaussian base distribution will generally lead to a worse performance compared to Reverse Diffusion Monte Carlo. This is due to the bias inherent in annealed Langevin dynamics, as discussed in Section 4.2 of our work. \nMoreover, for completeness we show in the table below, that the nested inner-outer loop algorithm based on annealed Langevin diffusion with a heavy-tailed base distribution underperforms our multiscale approach for the mixture of Student's t distributions benchmarks, when the number of energy evaluations is set to be the same for all methods.\n\nTable 1. Comparison of our MultDALMC method with a nested inner–outer loop algorithm, both using a heavy-tailed base distribution. Results are averaged over 30 random seeds. We evaluate two variants: one based on the modified Itô diffusion process (Equation (7) in our work) and another using an overdamped Langevin diffusion for the fast process and inner loop. Performance is measured by the entropy-regularised Wasserstein-2 distance.\n|        |    $\\mathbf{d=2}$   |$\\mathbf{d=10}$ | $\\mathbf{d=50}$ |\n|----------------|---------|-------------|-------------|\n| **MultALMC w/ modified Itô diffusion**          |   **0.85±0.06** |   **3.42±0.23**   |      **17.21±0.52**    | \n| **Inner-outer   w/ modified Itô diffusion**           |   1.12±0.29   |    5.01±0.57     |      20.39±1.14    |\n| **MultALMC   w/ overdamped Langevin**           |   0.91±0.11   |    3.97±0.28     |      18.56±0.69    |\n| **Inner-outer   w/ overdamped Langevin**          |   1.34±0.20   |    5.62±0.81     |      21.88±1.03    |\n\n> How is it ensured that hyperparameter optimization was carried out with comparable computational resources for the compared sampling methods? What were the hyperparameter grids used for each method and benchmark?\n\nBelow, we provide the hyperparameter grids used for our proposed algorithms as well as for all baseline methods. The design of these grids was carefully guided by the original papers introducing each method. We allocated sufficient computational resources to all baselines to ensure that the selected combination of hyperparameters was optimal, some methods required more extensive tuning than others. More specifically, each hyperparameter configuration was evaluated over 30 independent random seeds to account for variability and ensure robustness. When a method exhibited high sensitivity to certain hyperparameters, we expanded the grid around high-performing regions for finer exploration.\n\nBesides, it is important to remark that to ensure a fair comparison we used the same initialisation based on a standard Gaussian distribution for our methods and all baselines in Tables 1 and 2. However, this initialisation is not optimal for the SLIPS algorithm [2], which benefits from a more tailored initialisation. For completeness, we also report the performance of SLIPS with optimal initialisation in Table 7, Appendix D. This optimal initialisation assumes knowledge of the distribution’s scalar variance which is unavailable in practical scenarios.\n\n**Hyperparameter grids for the different algorithms**\n- **SMC and AIS**: The hyperparameter we tune is the number of MCMC steps $n_{\\text{MCMC}}$, its value is selected from the interval $[20, 100]$ using a grid step of 4.\n- **ULMC**: The step size is chosen within the following grid $h\\in\\(0.001,\\dots, 0.009, 0.01, \\dots, 0.09, 0.10, 0.11, \\dots, 0.20\\)$. The friction coefficients $\\Gamma_{\\min}$ and $\\Gamma_{\\max}$ are chosen within the grid $\\(0.001, 0.005, 0.01, , \\dots, 0.09, 0.1,\\dots, 0.9, 1, 1.5, \\dots, 10\\)$. The mass matrix is set to the identity.\n- **PT**: The grid for the step size is the same as that of ULMC. The number of chains is selected from the interval $[1, 20]$ with a grid step of 2. The corresponding temperatures are chosen to be equally spaced on a logarithmic scale, with the minimum temperature fixed at 1 and the maximum temperature chosen from the grid $\\(10^2, 10^3, 10^4, 10^5, 10^6, 10^7\\)$.\n- **DiGS**: The grid for the step size is the same as that of ULMC. The number of noise levels $K$ is selected within the interval $[1, 20]$ with a grid step of 2. The maximum and minimum noise levels, $\\alpha_1$ and $\\alpha_K$, respectively, are chosen from the grid $\\(0.05, \\dots, 0.09, 0.1, 0.2, \\dots, 1\\)$. Finally, the number of Gibbs sweeps $n_{\\text{Gibbs}}$ is selected from the interval $[50, 500]$ using a grid step of 50.\n- **RDMC**: It has a single hyperparameter $T$ which we search within the grid $T \\in \\(0.001, 0.005, 0.01, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.5\\)$ as suggested in [1].\n- **SLIPS**: Following [2], the number of MCMC steps is selected from the interval $[20, 100]$ with a grid step of 2. The value for $t_0$ is chosen by approximate equal log-SNR spacing in $[-5.0, 0.0]$. The value of $\\eta$ is selected from the interval $[4,6]$ using a grid step of 0.2.\n- **Ours**: The grid for the step size is the same as that of ULMC. The scale separation parameter $\\varepsilon$ is chosen within the grid $\\(0.001, 0.005, 0.01, \\dots, 0.09, 0.10, 0.11, \\dots ,0.20\\)$. The friction coefficients $\\Gamma_{\\min}$ and $\\Gamma_{\\max}$ are selected within the grid $\\(0.001, 0.005, 0.01, , \\dots, 0.09, 0.1,\\dots, 0.9, 1, 1.5, \\dots, 10\\)$. The mass matrix is set to the identity. The parameter $\\tilde\\lambda$ is chosen from the interval $[0.3, 0.7]$ using a grid step of 0.05.\n\nWe look forward to engaging further during the rebuttal period to address any remaining concerns you may have. Thanks again!\n\n[1] Huang, X. et al. (2024) Reverse Diffusion Monte Carlo. The Twelfth International Conference on Learning Representations.\n\n[2] Grenioux, L. et al. (2024) Stochastic Localization via Iterative Posterior Sampling. In Proceedings of the 41st International Conference on Machine Learning."}]

Author rebuttal:

Meta-review: This paper proposes a novel sampling method from a given unnormalized target distribution that is high-dimensional and has many modes through a diffusion path. In specific, it applies a fast stochastic differential equations (SDE) to efficiently estimate the drift term of the averaged intermediate score (target slow SDE) without any training. Experimental results on a number of toy examples and real-world datasets show that the proposed multi-scale dynamics-based sampling (coupled slow and fast SDEs) outperforms baseline samplers.

Overall, the use of multi-scale dynamics for diffusion-based sampling seems to be technically sound and novel. Most concerns raised by the reviewers are well addressed by the authors, including a clear description of the target sampling problem different from generative modeling, a computation cost, hyperparameters’ sensitivity, a limited theoretical analysis on the strict constraint and assumption. Based on the consensus of all positive ratings between the reviewers, I would recommend the paper to be accepted.

However, I think that one main concern regarding the practicality on real-world datasets still needs to be further resolved. For example, given a learned energy function (even though it is not perfectly correct) on mnist or cifar, the sampling performance comparison between the proposed fast inner process and other baselines or variants would be necessary. Moreover, it is necessary to analyze whether the proposed method yields greater benefits as the dimensionality increases or the number of modes grows.

In addition, regarding the paper’s organization, rather than allocating limited theoretical analysis to the main body, it would be more appropriate to present the details of the proposed sampling methods (Algorithms 1 and 2) there and to further strengthen the experimental section.
Decision: Accept (poster)