Alternating Diffusion for Proximal Sampling with Zeroth Order Queries
Authors: Hirohane Takagi
Atsushi Nitanda
Venue: ICLR 2026 Poster (2026)
Keywords: Sampling
Diffusion-based Monte Carlo
Zeroth-order methods
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We introduce a zeroth-order proximal sampler with Gaussian-mixture score estimation that converges rapidly without rejection sampling.

Abstract: This work introduces a new approximate proximal sampler that operates solely with zeroth-order information of the potential function. Prior theoretical analyses have revealed that proximal sampling corresponds to alternating forward and backward iterations of the heat flow. The backward step was originally implemented by rejection sampling, whereas we directly simulate the dynamics. Unlike diffusion-based sampling methods that estimate scores via learned models or by invoking auxiliary samplers, our method treats the intermediate particle distribution as a Gaussian mixture, thereby yielding a Monte Carlo score estimator from directly samplable distributions. Theoretically, when the score estimation error is sufficiently controlled, our method inherits the exponential convergence of proximal sampling under isoperimetric conditions on the target distribution. In practice, the algorithm avoids rejection sampling, permits flexible step sizes, and runs with a deterministic runtime budget. Numerical experiments demonstrate that our approach converges rapidly to the target distribution, driven by interactions among multiple particles and by exploiting parallel computation.

Published as a conference paper at ICLR 2026
ALTERNATING DIFFUSION FOR PROXIMAL SAMPLING WITH ZEROTH ORDER QUERIES
Hirohane Takagi - Atsushi Nitanda2.5.4.
1 Graduate School of Information Science and Technology. The University of Tokyo, Japan
Institute High Performance Computing, Agency for Science, Technology and Research, Singapore
3 Centre for Frontier AI Research, Agency for Science, Technology and Research, Singapore
4 College of Computing and Data Science, Nanyang Technological University, Singapore
i ht atakagiliss s. u u-tokyo. ac. ip
8 at atsushi ni www.iiiinnnaaaaaaaaaaaaaaaa edu.sq
ABSTRACT
This work introduces a new approximate proximal sampler that operates solely with zeroth-order information of the potential function. Prior theoretical analy- ses have revealed that proximal sampling corresponds to alternating forward and backward iterations of the heat flow. The backward step was originally imple- mented by rejection sampling, whereas we directly simulate the dynamics. Unlike diffusion-based sampling methods that estimate scores via learned models or by invoking auxiliary samplers, our method treats the intermediate particle distribu- tion Gaussian mixture, thereby yielding a Monte Carlo score estimator from directly samplable distributions. Theoretically, when the score estimation error is sufficiently controlled, our method inherits the exponential convergence of proxi- mal sampling under isoperimetric conditions on the target distribution. In practice, the algorithm avoids rejection sampling, permits flexible step sizes, and runs with a deterministic runtime budget. Numerical experiments demonstrate that our ap- proach converges rapidly to the target distribution, driven by interactions among multiple particles and by exploiting parallel computation.
INTRODUCTION
Sampling from probability distributions (x) x e fx(()) is fundamental task statistics and ma- chine learning, with applications in Bayesian posterior inference and score-based generative model- ing. Methods like the Unadjusted Langevin Algorithm (ULA) (Roberts & Tweedie 1996 Durmus & Moulines 2017) and the Metropolis-Adjusted Langevin Algorithm (MALA) (Roberts & Rosen- thal] 1998] Roberts & Stramer (2002) are widely used. Their convergence under strong convexity
assumptions has been established in sharp nonasymptotic terms (Dalalyan) 2017 Dwivedi et al.,
2018). and it has further been shown that ULA enjoys exponential convergence in KL divergence under functional inequalities (Cheng Bartlett 2018 Vempala & Wibisono, 2019).
Beyond Langevin-type approaches, there has been growing interest in alternative sampling schemes with nonasymptotic guarantees. Among them, proximal sampling (Liang & Chen 2023b) intro- duces auxiliary distribution close to the sunget-typicall the Gaussian convolution of π (Lee) al., (021) and alternates conditional updates between the target and auxiliary variables. From theoretical perspective, each update can be interpreted as alternating forward heat-flow step and a reverse denoising step, which yields exponential convergence under suitable functional functional-iequuality
assumptions on the target distributions (Chen et al. 2022 Kook al.
2024
Wibisono 2025).
Despite this line of analysis, scalable implementations of proximal samplers remain challenging. Existing implementations rely local optimization of combined with rejection sampling (Liang) & Chen 2023alb) by Fan et al. 2023: Liang & Chen 2024), which necessitates small step sizes (i.e., weak convolution) to maintain acceptance and thus incurs many iterations and high overall cost.
"This work was conducted as part of research internship at Agency for Science, Technology and Research (AxSTAR), Singapore.

---

Published as a conference paper at ICLR 2026
Ideal proximal sampling
Approximate proximal sampling via alternating diffusion
Pk
Pk
Heat flow
Reverse dynamics
P
Surrogate dynamics
k+4
One step (discretization score estimation)
Pk+1
X Pk+1
π
π
π
π
Gaussian limit (not required)
Figure 1: Illustration of the ideal proximal sampling (left) and our approximation (right). Heat flow and reverse dynamics are defined between x and π but applied to intermediate Although these do reach their targets in one step, the ideal version attains exponential convergence. Compared to the rejection sampling-based implementation of proximal samplers, our approach allows for larger step sizes (i.e., stronger convolution), which reduce the iterations to reach the target distribution.
These bottlenecks have spurred diffusion-based Monte Carlo, which simulates denoising stochastic differential equations (SDEs) (Huang et al. 2024a He et al.) [024). This raises natural question: can proximal sampling be implemented in its theoretical form and in scalable way, directly through Gaussian convolutions and diffusion processes, without relying on rejection sampling?
Another practical consideration is the efficient sampling of many particles in parallel. In such set- tings, particle-bar based algorithms that introduce additional interaction terms or gradient-flor struc- tures promote faster mixing while maintaining diversity (Liu & Wang 2016/
al. [2200] Boft & Vanden-Eljjden 2023 Lu et al. 2024 [lin et al., (0055). However, proximal sampling was originally formulated appeal and the practical demands of multi-particle sampling.
single-particle iterative framework, leaving gap between its theoretical
1.1 CONTRIBUTIONS
In this work, propose and analyze diffusion-based approximation of proximal sampling. The key idea is to interpret the auxiliary samples as forming a Gaussian-mixture approximation and to exploit this structure for approximate time-dependent estimation in the denoising dynamics. Our main contributions are as follows:
We introduce a new algorithm that serves approximate proximal sampler. It is learning-free, gradient-free with respect to f, and rejection-free with a fixed runtime.
We extend the theory of proximal sampling to show that our diffusion- besttusioo--beeddd approximation inherits comparable convergence rates under suitable assumptions. This implies that our algorithm, which alternates perturbation and denoising, converges to the target distribution without requiring initialization from the Gaussian limit, unlike standard diffusion models.
We provide empirical evidence that our method converges rapidly to representative tar- gets compared to existing implementations of proximal sampler, and that marti-particle extension improves both wall-clock efficiency and sample diversity.
2 PRELIMINARIES
In this section, we provide a brief overview of proximal sampling, which underlies our proposed method, along with its convergence properties.
2

---

Published as a conference paper at ICLR 2026
2.1 DEFINITIONS
Let and y be two probability measures on R with μ & u. We define the Kullback-Leibler (KL) divergence, the Rényi divergence of order q 1, and the relative Fisher information
H1( = -f dj du
9 1-1 log /
1
d dv
du
du, Ruj((1) =
du,
Ju(1) =
V
log
djc
dv
Note that setting = 1 yields R1, = Hp),). We say that v satisfies a log-Soboley inequality (LSI) with constant CLSI > 0 of (p) - CCsssJ(u).
2.2 PROXIMAL SAMPLING
Let :R4 I be potential function, and consider the target distribution x (x) exp(-
For a given step size h O, following Lee et al. (2021). we define joint distribution on IR x R as
x -- 2/12
X.X (I, x exp
(f((x)
(1)
2h
Marginalizing over y yields π (x) x I xxxXY(2,, dy. defines new distribution TT" as n (1) fexp (-f(x) 2 y/1²) da. This corresponds to the Gaussian convolution x² = (l,,
On the other hand, marginalizing over
The proximal sampler with step size h iteratively updates a current particle æ4 E Rd as follows:
Dx+4 (N(,,
Ik+1 x
(3)
where x (x) exp(-f(x) | The restricted Gaussian oracle (RGO) (Lee
24
al. 2021 assumes access an oracle which enables sampling from
(1r).
Several works have implemented the RGO (Liang & Chen 2023a Fan et al. 2023) through rejec-
tion sampling, achieving an expected complexity of O(1) under regularity conditions on f. How- ever, this requires choosing sufficiently small depending on f and d, which in practice leads to many iterations, and the computational cost further fluctuates due to randomness.
2.3 CONNECTIONS TO DIFFUSION PROCESSES
Chen et al. (0022) established an improved convergence analysis of proximal sampling by inter- preting each update through dynamics interpolated by SDE, illustrated in the left panel of Figure I] The forward step (2) amounts to Gaussian convolution, which corresponds to evolving the heat equation or its associated SDE, for E [0, h).
dz, dB,,, 0,44 = A
(4)
where B, is a standard Brownian motion. Writing P. for the Gaussian convolution kernel (O0 tla), π have (4 MrP and in particular x² ==== PA. The backward step (3) can be viewed as the reverse operation of the forward step conditioned on Zn- It governed by the following SDE starting from 7 74, and the corresponding Fokker-Planck equation; for 1. E h),
dz V log(( ph_p))) dt + dB, dep div log( x PPr-)) + AApE. (5)
where B denotes the backward Brownian motion. Letting Q, denote its transition kernel, we obtain 14p 1 Qa π xx)= (x) 15 (y) dy. In particular, we have x. =
Hence the two steps in proximal sampling can be viewed as SDEs between π and Y with P: and Q as their transition kernels. The proximal sampler applies the dynamies 4 and (5) with the start distributions p law(zz) and P := law(y1-1/22 then evolve through Gaussian convolution conditional integration. Although single step does not reach π Y or π iterating the forward-backwwrr backward scheme leads to contraction towards the target distribution, made precise in the following theorem.
the k-th iteration, respectively, which
3

---

Published as a conference paper at ICLR 2026
Theorem (Chen et al. (2022). Theorem 3). Assume that x. satisfies LSI with constant CLsI- For any >> 0 and any initial distribution po the k-th iterate Pk of the proximal sampler with step size h satisfies, for q 2
R R. (p (1 + h/ (1s1))ak1
gg )
Theorem provides the exponential convergence guarantee of proximal sampling. Its proof follows by combining the forward and backward contraction properties each step stated in the next lemma.
Lemma
2022). Appendix A.4). Assume that π satisfies LSI with constant CLsI and p < n < For q 2 1, we have
(Chen et al.
xg
Ri (pPr)) R Ph (pPr) S
(1 h h
R (p')
Rg (p'Qn) Rg Qn (d(Q))
h1+1/C1s))2
In this paper, we aim to realize the proximal sampler by directly simulating the associated SDEs. In the forward direction, convolution with xx can be simulated exactly by injecting Gaussian noise, so the bound (7) is directly applicable. In contrast, the backward dynamics (5) cannot be simulated exactly, and our surrogate version of (5) will provide an approximation in place of (8).
APPROXIMATE MULTI-PARTICLE PROXIMAL SAMPLING
Building on the SDE interpretation in Section propose approximate the backward dynamics by replacing π with surrogate distribution constructed from the current particles. This enables Monte Carlo estimation of the score function using only evaluations of f, without requiring gradient or rejection sampling. At a high level, each iteration of our diffusion-based proximal sampler con- sists of (i) evolving particles by the forward heat flow (4) and (ii) applying the surrogate dynamics, which time-discretizes and approximates the reverse flow
3.1 SCORE ESTIMATION FOR SURROGATE DYNAMICS
We derive the update rule in the backward step, given the particles at the k-th iteration: Xe (xr))N and 1 {リ5 J. We approximate the backward dynamics (5) by replacing x with surrogate distribution (- 1 X₁) constructed from these particles, as follows:
N
(y1)
@x+1( Yx41/2, Xx) x
Σ
π
(x)
N
Ōk 1/2 (1) Xx)
1=1
where 9 Xx) N E.N(y; x, hla). The conditional distribution with the reweighting
term ensures that, N + x., sampling y 4x+1/2 recovers the target distribution x² By substituting the explicit form of म x and π we obtain
exp (-f(x) 24 13
(x+1(0 YK41/2, x N
dk 1111/2(31 Xx)
M( hla) 1x+1/2 (1) Xx)
exp((f((a) =: IN 1tt (2)xx(())) (10)
k+1/2
A+1/2
where yx x is the unnormalized density an N-component weighted Gaussian mixture (see Ap- pendix for the full derivation). ic+1 involves an inverse reweighting with respect to a Gaussian mixture 1/22 This reduces weights when yj is located in areas where X₃ is overly concentrated, while amplifying their importance in sparse regions. The resulting term, derived from the empirical particle system, may help prevent particle collapse and promote exploration. We later verify this effect in our experiments.

---

Published as a conference paper at ICLR 2026
Algorithm Zeroth-Order Diffusive Proximal Sampler
Input: potential function f: f: R, initial samples (a) J. step size h, iterations K, diffusion steps T, number of interim samples M, noise schedule {or fo217 Jt=D with o7 and o o min
b All operations for 1,.... N. = ,...,,,,,,, .,d 1, 1...., M are evaluated in parallel.
2: for == 0....... K -
#Step 1: Forward heat flow 4
yk+ 2k + vhsl ( with e itd La).
(3) (i)
3:
4: Initialize 2T Ik V with &
(0)
(
(()
(0) Lid N(O,
# #Step 2: Surrogate version of reverse cynamics 15
5:
for I T, T 1, 1,...: do
6:
Set At + of -7-11
7:
Compute o, maj = (2 ) Wis = wy ( by (10 and set wis t Whs/ Es was
Draw M samples (i.)) lid E, Wis.Vmin 0²).
8:
20 ~
Compute C₁ exp( () and set Ci f- ct,1/ E. E.
9:
# Update each particle using the Euler Maruyama step (i.) //o VAre ((i) with m() iid (1) + c,
10:
2/-1
(1)
2
N((0, (a).
11: end for
12:
Set (1)
k+1
20
(6)
13:
end for
1N
14:
return [ si=l
The surrogate score is defined as (2, V Vlo log( :) (Z, ). which serves the drift term
the surrogate reverse dynamics, replacing log (+* Pa .) This expression can be rewritten in expectation form using Bayes' rule:
in the backward dynamics (5)
Fn E 1(20) /C/
(z)2) = E k+1/2 /41/2 IN
(11)
oo
where b- C-= IN
kx11/2 (20 -1(o) dro-
(12)
k+ 1+1 (xx) (x01)((00) ) is a Gaussian mixture and the conditional law z z0 under the perturbation kernel (Nro of Ia), Bayes' rule implies that the posterior k++12
See Section A.1 for the derivation. Since 9N
distribution UN (wo =) is itself a Gaussian mixture:
k+1/2 (zolz)
IN
(13)
N( Mis (h o
where σ² = (h-1 + 2- m²(2) = d y + o, 22) wz(2) =
Q Xx))
(14)
The detailed derivation can be found in Section A.2 This formulation enables practical score esti- mation via Monte Carlo sampling from the Gaussian mixture without requiring model training.
3.2 ALGORITHM AND COMPLEXITY
The surrogate score [I involves an expectation with respect to the
posterior (3) In practice, we approximate this expectation using Monte Carlo estimation by drawing M samples for each particle. Combining the forward step (Gaussian perturbation) and the discretized surrogate reverse dynamics with 7 denoising steps, we obtain the iterative sampling procedure that approx- imates the proximal sampling. is worth emphasizing that the proposed method relies only on a zeroth-order oracle of and does not require gradient information, unlike cradient-based sampling methods such as Langevin Monte Carlo. A complete description of the algorithm is provided in Algorithm
5

---

Published as a conference paper at ICLR 2026
Regarding computational cost, straightforward implementation would require KTMN evalua- tions of the potential f for outer iterations. This complexity can be further reduced to KT under the assumption of parallel oracle that evaluate f simultancously multiple samples. Such a design naturally exploits parallel computation, making the method efficient in modern computing environments.
4 THEORETICAL ANALYSIS
In this section, provide theoretical analysis that essentially quantifies the deviation from the ideal proximal sampling caused by time and space discretization. Our method splits the macro step size into 7 substeps of length 17 h/T for the backward dynamies. For each l E ( let t = (E 1)m7 be the substep start time, and let Pt, denote the current law at time ty. The following lemma establishes the contraction of the KL divergence
a single outer iteration of 1 up to an error due to time discretization. We assume that л satisfies an LSI with (5) constant CLSI- Focusing on iteration k, write p law(x/")) (()) and Pk++122 1aw()) A complete description the assumptions and proofs deferred to (opppenddxx A
Algorithm
Lemma (Discretization-nlly one-step bound). Assume N, M - cc, so the backward update uses the exact score. Suppose further that within each substep E, the temporal variation of the score between le and le + L (1. E [0, ri) integrating over pactt is bounded by L2, the score of the reference measure Vt, is Lipschitz in space with constant Lu,s and there exists a uniform entropy bound (R) > 0 along the substep path. Then for u 2 1,
Har (2x H1x (px(Pp+)) (1+ (1+h/Crss))--/(eu)) au Aa A/4)) (1)
where
0u == =u -22 h1+h/C1s))
(K) = L,, (CLsI + h)Hl + (k))(k)
((k) 2nL L3 nL /d with Lu == sup Lu,e and L, sup Lea
The additional discretization error term in Lemma scales as Oh/TT) with respect the number of substeps T. Therefore, if T is chosen appropriately, the backward step essentially inherits the convergence rate of the ideal proximal sampler in (8).
When N, M co, we assume that at each substep start te the score estimator admits the bound A 2 supp:=.....? E,,,, [11 (:, - (-), where St, is the true at and (k) is the Monte Carlo estimator. The following result incorporates this condition.
N.M te)
Proposition (Main one-step bound with split errors). Let 2 Under the assumptions of Lemma with finite N, M, and the score estimation bound, the k-th iterate of Algorithm 7 satisfies
HaY (P% (Dx+1/2) H (pxt (1 h1+h/Casi))-//ee 1/(202) (A(( A (1
(K)
combined with the KL contraction the forward step (Lemma 1 7 shows that H2 (pt)) decreases geometrically in k. up to the accumulated approximation errors. Indeed, com- bining the forward contraction with the approximate backward bound yields
Proposition
Ha (Px..) & r H₂ (px) + au (A 1(4)),
where = (1 h/CLsI) 11((2)))) > Unrolling this inequality gives, for any i 2 1,
1 +
Hx (px
11x (,x((x))) + au
1-1
1=0

---

Published as a conference paper at ICLR 2026
In particular, if the per-step emorr satisfies the uniform bound AT A A for all j 2 then
using Ei=0 s - yields
H (pk)
pit H (xx
+ Out
r A.
Consequently, ensure H (px) suffices to choose
12H (/e)
log (p
c(r - 1)
2
and A S
logr
2(1 1/(2x²)) log(1 + h/CLsI)
20%
Equivalently, this yields the iteration complexity k = (llgg( H (00 )/8\/ logr) under the admis- sible error level A = O(e(r - 1)/()r)))
As for the score estimation error (k) reflects the Monte Carlo error arising from the finite sample sizes N and M. As shown in Section A.4 for instance, under bounded conditional variances of π obtain (11/N) error, and under bounded expectations of et and finite fourth moments under the Gaussian mixture proposals we obtain (1/MM) error, yielding A2 = (1/NN 1/M).
From Theorem the ideal proximal sampler requires about (2 log(1 h/ h/CLLs)) outer iterations times a logarithmic factor in the initial divergence. When CLSI >> h, this factor is approximated by O(Crsii))) suggesting that larger is favorable. While rejection-sampe sampling implementations of RGO suffer from upper bound on h, our method can take h large long T = (h) control the discretization error. This does not change the overall computational cost (outer iterations x T steps), but allows h reflect the global structure of T such inter-mode distances rather than only the local smoothness of f, which we confirm to be practically advantageous in experiments.
EXPERIMENTS
We design two experiments to showcase the strengths of our algorithm as an approximation to proximal sampling. First, we revisit the Gaussian Lasso mixture (Liang & Chen 2023b) to test whether our method accelerates convergence beyond RGO-based proximal sampling, especially via parallel particle updates. Second, we study uniform distributions over bounded, nonconvex, and disjoint domains, comparing our method with In-and-Our (Kook et al., 2024), the proximal sampler for uniform distributions originally analyzed for convex bodies. These experiments respectively demonstrate the benefits larger step sizes and the applicability of our approach when gradients unavailable.
5.1 GAUSSIAN LASSO MIXTURE
Setup. Following the experimental setting in Liang Chen (023b), we the target distribution as Gaussian Lasso mixture:
Vdel exp (-5 -0)0x(1)))))))))1))e -1) + 2d-1 exp expp-1140 (-114x(111
(x) xxxx
22/22T)
where Q = USU d 5, 11, 15, 16, 18), and U is an arbitrary orthogonal matrix.
We compare against the RGO baseline with 100 chains and step size h=1/135, exactly matching the experimental setting of [iang & Chen (002b)), where it outperformed ULA and MALA. Our method uses step size 1=1/10 under two settings: N=100 interacting particles, or N=1 run with 100 chains. For comparability, into one iteration, so that the parallel fevaluattion cost is comparable to ours.
thin the RGO bascline by grouping every 10 single-step updates
Convergence is measured by KL divergence to the target distribution, estimated using Büth et al. 2025) with a fixed reference 1000 particles from a long RGO run. At each evaluation, 1000 particles pooled across 10 successive iterations. The histogram in the right panel of FFiuree 3 shows this reference distribution. All experiments are repeated with 10 random seeds, reporting the mean and variance of the KI. estimate. Detailed settings given in Table 1 and Section B

---

Published as a conference paper at ICLR 2026
divergence particies aggregated every steps)
RGD 1th thimiing=11)
(diffusion, 20)
Wia Interaction
Edtiratedd diveggeze
200
400
600
800 1000
Iteration
Figure 2: Convergence of estimated KL divergence, averaged over 10 random seeds with shaded areas indicating variances, Our method (orange) outperforms both the proximal sampler with RGO (blue) and an ablated variant of our algorithm without particle interactions (green). achieves the same accuracy as RGO in about 10x fewer iterations x100 faster when accounting for thinning).
ITCO 3S92N!
RAG Infoernee
1.0
08-
28
0.--
02
30
Figure 3: Onc-dimensional marginals of x along the third coordinate. The red curve is the ground- truth. Our method around 100 iterations (left, orange) already matches the RGO-based sampler at ~~1000 iterations (left, blue), and with only 200-300 iterations (right, orange) it closely aligns with the reference obtained from sufficiently long run following Liang & Chen (2023b) (right, blue).
Our method converges substantially faster than the RGO baseline. With step size 13.5 times larger, our method reaches comparable KL divergence in 1100 iterations, while the proximal sampler needs ~950 iterations (ss 9500 RGO updates). Our method also surpasses the ablated variant with ==1 independent parallel chains. indicating that particle interactions are essential to accelerate mixing. The marginals in Figure indicate that our method is already approaching the target by 100 iterations (vs. ~50 for RGO), and it closely matches the long-run reference by 250 iterations. more detailed observation deferred to Section B (Figure 5), which shows that increasing the step size h leads to faster convergence.
Results. Figure 2) shows the convergence
Discussion. These results demonstrate that relaxing the stringent step size restriction of proximal sampling yields a practical gain of nearly an order of magnitude convergence speed. The benefit amplified when leveraging multiple interacting particles, which facilitate more efficient exploration of the mixture components. This effect can be explained by the surrogate target distribution in (9) which re- re-weights via (#x+1/2(y Xx)) to down-weight overpopulated regions and encourage exploration of sparser ones, thereby accelerating convergence heyond the step-size effect alone.
5.2 UNIFORM DISTRIBUTIONS ON BOUNDED DOMAINS
Setup. We evaluate method uniform sampling domain K C R3, and compare it with the In-and-Out (Kook et al. 2024). In-and-Out proposes U141 N(x,, hTa) and repeatedly resamples ø441 [( hTa) until x++1 E K or retry threshold reached. It converges under warm starts, but practical efficiency requires a convex
bounded, nonconvex, and disjoint
We take K the union of two disjoint solid tori in R Ti is torus centered (10,0, with major radius 10 and minor radius 1. and T₃ is a torus centered at (-13, 0, 0) with major radius and minor radius 1. We denote K T, U T2. Particles start from (0, Iz). Our sampler uses
8

---

Published as a conference paper at ICLR 2026
10 10
30
D
Out mod-3 2-201
Figure 4: Empirical distributions a 3, 10, 200 iterations with seed 0 the two-tori domain. In-and-Out (top row) finds Ti which overlaps with the initial standard Gaussian distribution, but fails to reach T2. Our method (bottom row) generates some particles outside the domain but gradu- ally drives particles toward 72, demonstrating its ability to explore both components.
if æ K and 100 otherwise, inducing an approximate uniform law. We run 1,000 particles and monitor their occupancies, visualized from the third coordinate direction. Details of the experimental setting are given in Table 2
the potential f())
Results. Figure 4 shows that In-and-OU Out converges uniformity inside the near torus T1, consis- tent with its guarantee under warm starts, but fails to reach the distant torus T2. Our method first fills the near torus, then transitions to the remote one, where many particles eventually accumulate.
Discussion. This experiment confirms that our algorithm works without gradients and can explore disconnected modes where exact proximal steps stagnate. The noisy score approximation facili- tates such transitions, akin to effects observed in diffusion-basee black-box optimization (Lyu al,, samplers include Projected Langevin Monte Carlo (Bubeck et al. 2018 MYULA (Brosse et al., 2017), and penalized Langevin dynamics based on distance to the constraint set (Gurbuzbalaban et al." 2024), as well proximal samplers for convex bodies ing projection or separation oracles Dang & Liang, 2025). All these methods critically rely on projection- or separation-type operations and are therefore limited to simple convex bodies where such mappings are computable. By contrast, our sampler only requires a membership oracle or simple outside penalty, making the zeroth-order oracle framework applicable broader range of domains.
2025 Other constrained-doman
CONNECTIONS TO OTHER SAMPLING METHODS
Diffusion-basse based Monte Carlo methods transport samples from a Gaussian initialization the target distribution. Some rely on training models for score estimation (Vargas et al. 2023 Richter &
Berner 2024), while others drive auxiliary samplers (Huang al., 2024a] He al. (024)) Several
works a so develop acceleration and correction techniques within this paradigm (Lu et al. 2022 Kim
Ye 2023 Li al., 2024). Beyond diffusion biffusion-bbsed pushforwards, related approaches construct
explicit density paths (Fan et al. 2024 Guo et al. 2025).
Our approach also specifies and employs SDE dynamics between two fixed distributions, but un- like the one-way pushforward paradigm, it repeatedly applies these dynamics. This removes the restriction of restarting from a Gaussian base by using a variance-expanding (VE) diffusion. In the context of variance-preservingg (VP) diffusions, sampling error partly arises from the discrepancy between the Gaussian equilibrium and the distribution obtained after finite-time mixing (Li et al., though the VE process does not
2023 Pierret & Galerne, 2025). Our analysis suggests that, reach the Gaussian equilibrium in any bounded horizon, convergence can still be ensured by alter- nating finite-time noise addition and denoising dynamics through repeated cycles.

---

Published as a conference paper at ICLR 2026
A particularly relevant comparison is with ZOD-MC (He et al. 2024). iiffusion-based method that performs denoising from Gaussian initialization. while inner proximal sampler is used for score estimation. Our method inverts this structure: instead of nesting a sampler inside a pushforward loop, we simulate proximal-style SDE dynamics directly. Moreover, unlike Z0D-MC, which requires access the minimizer of the potential- or, in practice, gradient information to locate local solutions for proximal updates-our algorithm op- erates solely with zeroth-order oracle information. Another related line SLIPS (Grenioux al.) 2024), which alternates between denoising noisy observations and updating auxiliary variables. This alternatig scheme resembles ours, but the overall iteration follows one-way dynamics that converge to the target distribution the time horizon grows.
ZOD-MC, the outer loop transports samples
Finally, our method is also close in spirit to Diffusive Gibbs Sampling (DiGS) CChnn et al. 2024). which performs Gibbs updates by alternating perturbation and denoising, updating both the state and an auxiliary variable each step. However, the designs differ: DiGS employs VP diffusion with auxiliary samplers such as MALA and provides convergence guarantees, whereas method uses VE diffusion consistent with proximal sampling, yields provable contraction guarantees, and exploits parallel multi multi-particle computation for lightweight score estimation.
CONCLUSION
We have introduced ceffusion-baseee approximation of proximal sampling, which simulates the backward SDE using only zeroth-order information of f. The key idea is to approximate the in- termediate particle distribution by a Gaussian mixture, enabling score estimation without auxiliary samplers or additional model training.
Our method remains within the proximal sampling framework, while also being closely connected to diffusion-based Monte Carlo methods. Our theoretical analysis shows that it achieves comparable when discretization and score estimation errors properly controlled. Unlike rejection-sampling implementations of RGO, which are constrained to small step sizes by local properties of f. our algorithm may accommodate larger step sizes that reflect global features such as inter-modd mode distances.
convergence
Finally, incorporate interaction terms among multiple particles, which empirically accelerate convergence. A complete theoretical characterization of these interactions, together with efficient parameter tuning across iterations, remains an important direction for future work. More generally, understanding how such empirical gains scale to higher-dimensional settings remains open chal- lenge, as the curse of dimensionality well-known issue in broader context of zeroth-order methods and rejection sampling.
ACKNOWLEDGMENTS
This research is supported by the National Research Foundation, Singapore and the Ministry of Dig- ital Development and Information under the Visiting Professorship Programme (award number AIVP-2024-004), Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not reflect the views National Research Foundation, Singapore and the Ministry of Digital Development and Information.
10

---

Published as a conference paper at ICLR 2026
REPRODUCIBILITY STATEMENT
We provide complete assumptions and full proofs of all theoretical results in Appendix A The de- tailed settings of our experimental evaluations are described in Section 5| Section In addition. the source code for reproducing the empirical results is available at https: gi hub. com/ ht 2od-ps/
REFERENCES
Nicholas M Boffi and Eric Vanden-Eijnden Probability flow solution of the fokker-planck equation. Machine Learning: Science and Technology, 4(3):035012, jul 2023. doi: 10. (0.1088/2632-2155)
ace2aa. URL http::ttp///.c....m. //dx::://dd.. /x doi www..org//000 112.0..1.8.000
Nicolas Brosse, Alain Durmus, Éric Moulines, and Marcelo Pereyra. Sampling from a log-
concave distribution with compact support with proximal langevin monte carlo. In Satyen Kale and Ohad Shamir (eds.), Proceedings of the 2017 Conference on Learning Theory. volume 65 of Proceedings of Machine Learning Research, pp. 319-342. PMLR, 07-10 Jul 2017. URL https: /ttpss//rroceeeee ieg.eediin.n..... ngs. m www.ppesss/////pe wt
htmml
Sébastien Bubeck, Ronen Eldan, and Joseph Lehec. Sampling from a log-concave distribution with projected langevin monte carlo. Discrete & Computational Geometry, 5944)7777-78, Jun 2018.
ISSN 1432-0444. doi: 10. 10.10007/000454--18899992 101454-018-9990
URL
Carlson Moses Büth, Kishor Acharya, and Massimiliano Zanin. infomeasure: a comprehensive python package for information theory measures and estimators. Scientific Reports, 15(1):29323, Aug 2025. ISSN 2045-2322. doi: (0.1008/41198-002 14053-5. BtLhttts: //do: /
10.1038/ 9115/8-19988-8 025--4
Djalil Chafai. Entropies, convexity, and functional inequalities, On t-entropies and --Sobolev inequalities. Journal of Mathematics of Kyoto University, 44(2):325 363, 2004. doi: 10.1215/ kjm/1250283556. URL Attt::://bl..orgII0..21122
Wenlin Chen, Mingtian Zhang, Brooks Paige, José Miguel Hernández- and David Barber. Diffusive Gibbs sampling. In Ruslan Salakhutdinov, Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp (eds.), Proceedings of the 41st Inter- national Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pp. 7731-7747. PMLR, 21-27 Jul 2024. URL actps: pprcee://prcceedigg... mllr. press/v235 che32//che22255 html
Yongxin Chen, Sinho Chewi, Adil Salim, and Andre Wibisono. Improved analysis for proximal algorithm for sampling. Po-Ling Loh and Maxim Raginsky (eds.), Proceedings of Thirty Fifth Conference Learning Theory, volume 178 of Proceedings of Machine Learning Research, pp. 2984 3014. PMLR, 02-05 Jul 2022. URL https: http:///oooeeee...c.cce chen22c. html
/ mlr ptwss//11ss/ wtt.s////r........s////////
Xiang Cheng and Peter Bartlett. Convergence of langevin meme in kl-divergence. In Firdaus Janoos, Mehryar Mohri, and Karthik Sridharan (eds.). Proceedings of Algorithmic Learning Theory, vol- ume 83 of Proceedings Machine Learning Research, pp. 186-211. PMLR, 07-09 Apr 2018. URL hittps:/pproceed: wtt.g/asmrr....../.......m.. /vteesssvv... v cheng 8a httm1
Arnak Dalalyan. Further and stronger analogy between sampling and optimization: Langevin monte carlo and gradient descent. In Satyen Kale and Ohad Shamir (eds.). Proceedings of the 2017 Conference on Learning Theory, volume 65 of Proceedings of Machine Learning Research, pp. 678 689. PMLR, 07-10 Jul 2017. URL https: http://rooeediingsssmmmrr.rreeee da a lyanllnllll 1aan17a.htm html
Thanh Dang and Jiaming Liang. Oracle- based uniform sampling from convex bodies. arXiv preprint arXiv:2510.00 0209293, 2025.
11

---

Published as a conference paper at ICLR 2026
Alain Durmus and Éric Moulines. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. The Annals of Applied Probability, 27(3):1551 1587. 2017. 16- 16-A2P. URthttpss://
doi: 10.1214/
org.
Raaz Dwivedi, Yuansi Chen, Martin J Wainwright, and Bin Yu. Log-concave sampling: Metropolis- hastings algorithmss are fast! In Sébastien Bubeck, Vianney Perchet, and Philippe Rigollet (eds.), Proceedings of the 31st Conference On Learning Theory, volume 75 of Proceedings of Machine Learning Research, pp. 793-797. PMLR, 06-09 Jul 2018. URL https: https::/proceeddiii /ttp://eoocec....../....
mlr. wwwpppeessss... httess/v15/wwivddd dwivedill Ba. html
Jiaojiao Fan, Bo Yuan, and Yongxin Chen. Improved dimension dependence of proximal algo- rithm for sampling. In Gergely Neu and Lorenzo Rosasco (eds.). Proceedings of Thirty Sixth Conference on Learning Theory. volume 195 of Proceedings of Machine Learning Research, pp. 1473-1521. PMLR, 12-15 Jul 2023. URL htttps://rroceddnggs.. wttppp..t.....////////
Far htmll
Mingzhou Fan, Ruida Zhou, Chao Tian, and Xiaoning Qian. Path-guided particle- based sampling. In Proceedings of the 41 International Conference on Machine Learning, ICML'24. JMLR.org, 2024.
Futoshi Futami, Issei Sato, and Masashi Sugiyama. Accelerating the diffusion ensemble sampling by non-reversible dynamics. In Hal Daumé III and Aarti Singh (eds.). Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Reseurch, pp. 3337-3347. PMLR. 13-18 Jul 2020. URL https: / /ttp://pprcece....//oep
mlr. www.ppees//e. 222as//111//Ettaa222aaa.a.a
am1
Louis Grenioux, Maxence Noble, Marylou Gabrié, and Alain Oliviero Durmus, Stochastic localiza- tion via iterative posterior sampling. In Forty-first International Conference on Machine Learning,
2024. URL http://///ppeeeeeeeeeeeee///aana
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. haattp/e//peneeeeeeeeeeettttttoo
Mert Gurbuzbalaban, Yuanhan Hu, and Lingjiong Zhu. Penalized overdamped and underdamped langevin monte carlo algorithms for constrained sampling. Journal of Machine Learning Re- search, 25(263) 2024. URL http://:: tttp://apppe/////2222111122 1443.html
Ye He, Kevin Rojas, and Molei Tao. Zeroth-order sampling methods for non-log-concave dis- tributions: Alleviating metastability by denoising diffusion. In A. Globerson, L. Mackey, D. Belgrave, A. Fan, U. Paquet, J. Tomczak, and C. Zhang (eds.), Advances in Neu- ral Information Processing Systems, volume 37, pp. 71122-71161. Curran Associates, Inc., 2024. URL https://prrcceedinggs http://ew.e....ao/////////////// http://pppeeeeepppppppeeeeeep f1les/paperr/2247 co cc/paper
fi
Xunpeng Huang, Hanze Dong. Yifan HAO, Yian Ma, and Tong Zhang. Reverse diffusion monte carlo. In The Twelfth International Conference on Learning Representations, 2024a. URI.
https: / /ttp///ppp//eeeeeeee/ee/ee/eeeee net. / forumm (o-.rumnii-IIIEMMMPF.
Xunpeng Huang, Difan Zou, Hanze Dong, Yi-An Ma, and Tong Zhang. Faster sampling with- out isoperimetry via diffusion-based monte carlo. In Shipra Agrawal and Aaron Roth (eds.), Proceedings of Thirty Seventh Conference on Learning Theory, volume 247 of Proceedings of Machine Learning Research, pp. //proceedings.. mlr. prwss / 47 hwtp44//hhang22aaaaaaam
PMLR, 30 Jun-03 Jul 2024b. URL https:
Vasily Ilin, Peter Sushko, and Jingwei Hu. Score-based deterministic density sampling. arXiv preprint (eXXxx22004181.00 2025.
Beomsu Kim and Jong Chul Ye. Denoising MCMC for accelerating reffussionbaseee generative mod- els. In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt. Sivan Sabato, and Jonathan Scarlett (eds.). Proceedings of the 40th International Conference Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 16955-169977 PMLR, 23-29 Jul 2023. URL www.httpss... htttps://procee::1 ngs. ml wwwppsssv20000 k m2. w html
12

---

Published as a conference paper at ICLR 2026
Yunbum Kook, Santosh S. Vempala, and Matthew S. Zhang. In-and-out: Algorith- mic diffusion for sampling convex bodies. In A. Globerson, L. Mackey, D. Bel- grave. A. Fan. U. Paquet. J. Tomczak, and C. Zhang (eds.), Advances in Neural In- formation Processing Systems, volume 37, pp- 108354-108388, Curran Associates, Inc., 2024. URL httpp::/pproceee/// ngs neurips http://parrrflllssppppp////00000 Ei
Paper Conference pdff
LF Kozachenko and Nikolai Leonenko. Sample estimate of the entropy of random vector. Problemy Peredachi Informatsii, 23(2):9-16, 1987.
Alexander Kraskov, Harald Stögbauer. and Peter Grassberger. Estimating mutual information. Phys. Rev. E, 69:0661 38, 38, Jun 2004. doi: 10. 20.1100/hyyseee......0 69.06669......... URL hollhttts://nni aps. 69. 066138
Yin Tat Lee, Ruoqi Shen, and Kevin Tian. Structured logconcave sampling with a restricted gaussian oracle. In Mikhail Belkin and Samory Kpotufe (eds.), Proceedings of Thirty Fourth Conference on Learming Theory, volume 134 of Proceedings of Machine Learning Research, pp. 2993-3050. PMLR, 15-19 Aug 2021. URL
.m presssvv333 3/tp hetp///331//eee....m
html
Gen Li, Yu Huang, Timofey Efimov, Yuting Wei, Yuejie Chi, and Yuxin Chen. Accelerating convergence of score-based diffusion models, provably. In Ruslan Salakhutdinov. Zico Kolter, Katherine Heller, Adrian Weller, Nuria Oliver, Jonathan Scarlett, and Felix Berkenkamp (eds.), Proceedings of the 41st International Conference Machine Learning, volume 235 of Pro- ceedings of Machine Learning Research, pp. 27942-27954. PMLR, 21-27 Jul 2024. URL
https: petps://prre http:///oooceeegg..iooe www.:mr....rp.mommmmmmmmm....eee wvw.sv2235//2 wtw.:/.ppress//2////52
Puheng Li, Zhong Li, Huishuai Zhang, and Jiang Bian. On the generalization properties of diffusion models. In A. Oh, T. Naumann, Globerson, K. Saenko, M. Hardt, and S. Levine (eds.), Ad- vances in Neural Information Processing Systems, volume 36, pp. 2097-2121 Curran Associates, Inc., 2023. URL /ttpsss/prooceednggg... neuripg.....u......ep..m........ cc ctt paper_ httpapprr__llseppprrraraa/ 2023/fil 2200:1LLe00000000000000000222112 Core pdf
Conference
Jiaming Liang and Yongxin Chen. A proximal algorithm for sampling from non-smooth potentials. In Proceedings of the Winter Simulation Conference, WSC '22, pp. 3229-3240. IEEE Press, 2023a.
Jiaming Liang and Yongxin Chen. A proximal algorithm for sampling. Transactions Machine Learning Research, 2023b. ISSN 2835-8856. URL https: http:://///eeeeeee.eeeeee wwt.or/wwwweeeeeentttt. forum?
Id=CkXOw1hf27
Jaming Liang and Yongxin Chen. Proximal oracles for optimization and sampling. arXiv preprint arXiv:2404.0223. 02239, 2024.
Qiang Liu and Dilin Wang. Stein variational gradient descent: A general purpose bayesian in- ference algorithm. In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016. URL hatps://proeediins..eei... cc/ppprrrrppppprrrrppppppppppe.. filess/aper/2016// b3a.b/bbaBbbcm hlle://aeeeieeee2332233000003333
pdf
Cheng Lu, Yuhao Zhou, Fan Bao, Jianfei Chen, Chongxuan LI, and Jun Zhu. Dpm-solver: A fast ode solver for diffusion probabilistic model sampling in around 10 steps. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Neu- ral Information Processing Systems, volume 35, pp. 5775-5787. Curran Associates, Inc., 2022. URL httpp///t////ooooooooooorr/ssssp
File /ele//
Janneng Lu, Yue Wu, and Yang Xiang. Score-based transport modeling for mean-field fokker- planck equations. Journal of Computational Physics, 503:112859, 2024. ISSN 0021-9991. doi: http://ciianeeeeer...ee/ess////2
URL https: wwww. scw.sceeeeiecd.rtt.
13

---

Published as a conference paper at ICLR 2026
Yueming Lyu, Atsushi Nitanda, and Ivor Tsang. Nonparametric distributional black-box optimiza- tion via diffusion process. In ICLR 2025 Workshop on Deep Generative Model in Machine Learn- ing: Theory, Principle and Efficacy, 2025. URL https:/ /otpneeveee/eeeeeeeeeee net forum? iid DOXsyHprRN
Emile Pierret and Bruno Galerne. Diffusion models for gaussian distributions: Exact solutions and wasserstein errors. In Forty-second International Conference on Machine Learning, 2025. URL
attps://openreet
net http:/ooommmmddxxxeeeeeee
Lorenz Richter and Julius Berner. Improved sampling via learned diffusions. In The Twelfth In- lernational Conference on Learning Representations, 2024. URL
G. O. Roberts and O. Stramer. Langevin diffusions and metropolii- algorithms. Method- ology And Computing In Applied Probability, 4(4))337-337, Dec 2002. ISSN 1573-7713, doi: 10. 00.103//A1003356111133
Gareth Roberts and Jeffrey S. Rosenthal. Optimal scaling discrete approximations langevin diffusions. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 60(1): 255-268, 1998. 2985986
ISSN 13697412, 14679868. URL http:://wwt /www. jater. orsttoboor//ttaale/
Gareth O. Roberts and Richard L. Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 2(4):341 - 363,, 1996.
Francisco Vargas, Will Sussman Grathwohl, and Arnaud Doucet. Denoising diffusion samplers. In The Eleventh International Conference Learning Representations, 2023. URL https:
//openreview.o net/forunm opterveeiee...eeefooooreooooeeae wtt.orommm..ssaaaa.aaaa..
Santosh Vempala and Andre Wibisono. Rapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 32. Cur-
wttp:/tppss/ppppceee............ ran Associates, 2019. /
files/
Andre Wibisono. Mixing time of the proximal sampler in relative fisher information via strong data processing inequality. arXiv preprint (tti::202200622533 2025.
14

---

Published as a conference paper at ICLR 2026
A PROOFS AND TECHNICAL DETAILS
A.1 DERIVATION OF THE SURROGATE ScoRE ESTIMATOR
We begin by substituting xx with
1 N
(3g
+ (or 1 4 X1.)
X 2s (2r)
9h+ 1 (15 Xx)
where (3 Xx) = N, hla). This follows from the fact that the forward heat
flow transports the empirical Dirac mixture t Ei oR, to the Gaussian mixture 9k+1/2
Using
2% 1 リo
xXIY=b( Sexp( (f(()- z 1/u k5 12) du
x (1;) x / eep((sS(u) # llu - Mo 12 du,
we obtain
exp(- f(x))
2h 1 - 131111))
(111 Yo+g, Xx)
N
gx+ (y/ (xx)
Nz: yj. hld) 1 (x)
k+1/2 (x) exp(-f(x))) exp(-f(xx) =: NN
1+1122 where ON is the unnormalized density of N-component weighted Gaussian mixture. As N 00,
lim 1+1/2 IN (x) = N(z,v,haa)
9x+1/2(y Xx) 2 Xx)dy const,
N
and therefore ( Yk+1/2, Xx) - n (.).
Since 0c+1 can be regarded as the product of (i) weighted Gaussian mixture term and (ii) the exponential factor exp(-f), the surrogate score function & 5,(2,) == Ph-)))) can be computed using only the particles and evaluations of f. Following the derivation in Lyu et al. 2025). we obtain
(
8x(2) V log
IN (xx) p(z xo) e ((0o))
1+1/2
J J k+1/2 (2O) Vplz xo) f(xo) dæo
f AN
k+1/2 (Io) p(z
xo) e (aa)
Here
(2To77)a/2 exp (- (_1|2 -12--200 20l1
p(2 x0)
of
20?
so that
I0 2
Vp(z wo) = p(z 200)
o²
Substituting into the above expression yields
f 9N (xo) p(2 xo) 2og e f dro
k+1/2
(2) =
f IN
1/2 (Io) p(= ro) e (m d
f 9N (wo (2
k+1/2
2) " e (
f IN
1/2 (No 2) e (oo
which is exactly the score function of the surrogate dynamics
15

---

Published as a conference paper at ICLR 2026
Remark 1. Lyu al. (2025) consider distributional black-box optimization for maximizing f. When (xo) Gaussian mixture, they show that setting ( x p* (20) exp(-f(au) (xx)) / as the denoised distribution at time 0 in VP diffusion, the score each time t can be computed directly from samples of the posterior Gaussian mixture, without any auxiliary sampler.
(zo)
A structurally similar mechanism appears in particle-based proximal sampling. To sample from T xxx (x) x cxp(-f(x)) (in practice, by moving samples from 0x+1/2(B) y Xx) toward
00+1( (x Yx+1/2, Xx )), we place the denoised distribution gN / (x) exp(-f(x)), where gN k+1/2 is the Gaussian mixture obtained by smoothing the empirical particles Yx+1/2. The score identity yields a backward update that uses only zeroth-order evaluations of f.
A.2 DERIVATION OF THE (UNNORMALIZED) GAUSSIAN MIXTURE POSTERIOR
When sampling from an N-component Gaussian mixture, first draw component index accord- ing to the relative weights and then sample from the corresponding Gaussian component. For this reason, we often ignore constant multiplicative factors of the mixture distribution for simplicity.
k+ 1+1/2 (x0) x Ej Oky ( My, hLa) with aj = 1/ Xx). By Bayes' rule,
We write 9N
*+ (xo) p(= wo)
k- 2 (v0 12) = IN
9N
S AN (u) p(2 u) du
Σ aj ay.NaNNNNN (x ys h1d N of
1 aj NN, Uj h(a) (2 u du
N
Σ N Uj hla) N(); Ta). j=1
(15)
We now complete the square in the exponent of the product N NZz: yz. hla) N (2; 2o, of (a) viewed as function of zo. Using the identity
1 - -00 1oll? +
1 20?
|| - zoll²=
20 -00 m2(2)11 + 2(h of) 112-1111 15l2.
2h
with
σ² = -1-1++022
hot m,(x) == 02 02 -0²²)
of
h
h + of
+ o
V/ +
+ t 2,
we obtain the product-of-
factorization
NZzo yo hla) (2; (2;20.0f(a) AN (h
o2),)
Substituting this into 15 yields
N (a
k+ IN (no 2) &
k+ 1/2 (ro 2) is again Gaussian mixture with common covariance σ² and updated Therefore IN
1/20
means my (2). Writing the relative weights as
(x)x)) = a, (2z; yo (h + oa) = 0K+1/2(15 Xx)
(h (h+
obtain the desired posterior decomposition
N
UN (wo
k+ 3 x) α
Σ
wy
j=1
Finally, note that using W(z) E3 wj(z), the normalized posterior distribution be written
k+ 9N
(zo 22)=
Σ
wy
w(=)
16

---

Published as a conference paper at ICLR 2026
A.3 ANALYSIS FOR DIME-DISCRETIZATION ERROR
Our algorithm simulates the surrogate version of diffusion process (5) by dividing the dynamics over horizon h into 7 steps. We first analyze the discretization error under the assumption that the score function is perfectly estimated (i.e., N, M - oo). Each step simulates the time evolution over an interval of length 7j h/T. corresponding to the segment (e 1)m + [[e- 1)7, Erj] with t E [0. for f= 1, T. In the ideal dynamics (5). the drift term is time-dependent, whereas in the discretized scheme we approximate by fixing the drift at the beginning of each step. Following
an argument similar to that in Vempala & Wibisono (2019), we obtain the following result.
Lemma 3. Fix step index 8 E (1,...,T) T) and let 1) == h/T. Let (ut))eto.n) denote the ideal backward marginals in (5) at elapsed time (e - 1)n + t for E [0, m/I. and (He)eelo.n) denote the frozen-drift approximation within this step. Then for 2 and all E [].
d di Hv. (pe) S -(2-1t
1
4w² 4
(u + u² E (155(()))))) 5e(2))11].
(16)
Equivalently, integrating over t [0, yields
Hv, (pr)) - (Ho) - - 2 4u2
1
S." Ju (He) dt + u² [- ('Emocl1BB(2a) - {x(x))]? dr.
Proof. For convenience to analyze the behavior at (E - 1)n + with t 10 which corresponds to
the time interval [[E- - 1)1,, Enj in original backward process defined in 5). we restate the SDE and the associated Fokker-Planck equation as: for & [0. 7).
dzz = (z ) dt + dB,,, law (Z = V/-
=1,
0,y, div (v) +
1 Au = V (vr(-5e (V-
log
2
where we define the score function s, == V log (T Ps Ph-((e-1aot).
In contrast, the time-discretized approximation corresponds to the process
dæ, (0)) di dB,, law(zz) = Mt:,
(17)
where the drift is frozen at the beginning of the step.
The associated Fokker-Planck equation is
1
aspl(E | 20) div((u(ze I 20) So(2o))
Ape(ze I 20)
2
= V - (palEe (-)0-5) +
v Ploguelz 20))).
Averaging over 20 ~ No. this becomes
1
log pe) (e))
0.14 V
(50(20) | 2x] + 2
The time derivative of the KL divergence is then d
d Hv, (14)
dt / 14 log Me dz V1
dt
= / [((h)t log 1t + (Bi14 14 Oru )] dz
/|v (-Eo-Ipp [5((0) zel + V log 2Vlogue)
log Ht
ve
Me v. (vr( (v((-3+- v Viog dz
us
/l (Mel (-Eyoje (30(20)
1 1
log m1), v log Hi
2x] +
+ (ve
+
1 log vt).
H4 v
log #) dz u1
2
= { 1 V log
1t + Byog (So(20) 2x] 3t, Vlog
vi uh dz.
ii
17

---

Published as a conference paper at ICLR 2026
Here, the third equality substitutes the Lokker-Plannee equations for depi and drrn and also uses that the integral of 0.j4e vanishes due to mass conservation. The fourth equality uses integration by parts (assuming sufficiently fast decay infinity) and the identity.
Simplifying, we obtain
d
1 -- (me) + Eus (EHon (50(20) |
|xo] St, V log >
Hi V4.
Hu, (w))
dt
Applying (a, b) S u² u²|1a|12 1/2+ + 1. 1011, we conclude
1
1 H|, Me, (me) 5 4u²
Jv. (pr) + u² B (1050)))))) - {=(=))l?)
This establishes the discretization error bound of one diffusion step in terms of the deviation between the frozen score ig and the ideal time-dependent score &t-
Assumption (Smoothness of the interim distribution). The interim distribution π² Pn- 1)n is Ly-smmothh, i.e,, the potential gradient V log(T< (h-(e-1)m) is L--Lipschitz
Assumption (Lipschitz condition along the time direction). The expected score satisfies Lipschitz-rype condition along the time direction:
Eu, (1)))(2) (2+))1] S C2,,t².
Corollary 1. Assume Vu == x P(e-177 1ym satisfies LSI with constant Cesi(Vo) and under Assumption and Assmpliinn 2
Hv, (Pi) S (
1
Jv. (He) +2² (4n2Lt CaS((Po) Heo (Ho) + ndC),
(18)
dt where C sup 2TT
2
4u2
12 + 10//e/d.
(19)
Proof. We bound the second term of the right hand side in the inequality in Lemma using sumption andd Assumption 2]
As-
E (20.=)) (20) - 3. (z+))l < Eyg, B =.) [1]50(=0) - =0 + E,.. []10((z) =x(2))1
< L (o.0 [1 =x11?] + C2,e 12.
Under the discretization update
2t = 20 + 30 (2o) 1 VLg, 5~ S~N(0,
we have
(20. [1l20 - zel1
VeEl
= 2 Ep (ooolll5((z)) [=)=(=)))11]1 td
< & (LSi(vo) Hup (Mo) + 2d + td.
The last incquality comes from Lemma 12 in Vempala & Wibisono (2019) with AathhAssumption I and Pg satisfying Talagrand's inequality with constant (Ls) Putting them all together, obtain 188 where is independent
by taking the supremum as (19).
Proposition 2 (One-step bound for the diffusionnapproxxmmaed proximal sampler without score timation error), Assume Assumption 1 and (Assumption 2) hold, and that xx satisfies an LSI with constant CLSI(T Let C be (19). Let overall step size h 0 be split into T steps with == h/T. In the regime N,M + oo where the score estimation error vanishes, the distribution at the 1i-th iteration P satisfies
H& (px)
Hgx (px)) <
+ 2u -2222 (T ((L+ (x - CL
+ ndC)
for u 2 1, where is the supremum of KL divergence between the interim distribution of the backward denoised distribution at timestep En and n² Qin and T. = sup, Tue
18

---

Published as a conference paper at ICLR 2026
Proof. As shown in Corollary 13 of Chaafa (2004). the logarithmic Sobolev constant satisfies
CLSI( P.) s CLSI(T*)) + t.
For E [0, n), the law corresponds to the ideal denoising path at time (e - 1)n + I, namely Vt X (h- Hence its logarithmic Sobolev constant can be bounded as
CLS)(v)) CLSI( + (e 1)m t =: t.
From Corollary combined with LSI, we obtain
dt Hv, (4n) S - 1 121/(2n²) ae - Hv, (1a) + In (L)))(P))), (1) + ndC].
Applying Gronwall's inequality, where multiply by the integrating factor: (ar - 1)- (0-+11/(22))
di { (a t) -))--22))2(22)) He, (wi)} (ar t) -1+1/(222) u² 4n/ Lm) (mo) + mdc
and integrating [0, m/,, obtain
(ave - 1) -111+//(2u)) Hv, (1n) Otg 1+11(22"²) H (po)
2u 1/(200 1/(2n)) (at - 11) //222)))))))) 4n7²
L LÇCLSI
Hro (HO) + ndc).
Recall that Me, i/ defined in the 8-th step of the time-discretized diffusion. Using the fact that ph,, My are equivalent to po: vo in the (E + 1) step and og - = 01gt1, we iterate this inequality over L...., obtain
Cis((Tx)-111/22) HAN (pktt) - (Cisi( h)) H.Y (Pk+1/2) CIsi( ) L + h) + ndc)
(Cis((C) h)) 1(22))
where H == sup, Hn, (Ho), the supremum of the KL divergence between the updated distribution in each diffusion step and the corresponding ideal distribution without discretization error and L = sup, Lt.
Equivalently,
HaY (06+2/2)
H (ph+1) <
(1 (T=n/
+ 2u 22tCCis( ((1 + h/Cr.srTTT 1) - 1) -1)
L. ) +h)H + sdd)
This corresponds to the time-discretized version of the KL contraction in S with 1. We complete the proof by applying inequality (7) in Lemma 1 with
19

---

Published as a conference paper at ICLR 2026
A.4 ANALYSIS FOR SCORE ESTIMATION ERROR
In this section, we evaluate the error of score estimation under a set of assumptions. In particular, prove that the estimation error scales as O11/N + 1/M), which converges zero in the limit of N, M + 0o.
Recall that (1 Xx)) 3/4 E; Nix 24,, hla) and we further define (unnormalized) Gaussian mixture density
IN x + (r U = 1 Na hla)
UK+ (1/ Xx)
which is conditioned on U = 1/22 X1) s.t. X ~N p Yo41/2 ~N (x+1/2(-Xxx)
At time f, the true score in the surrogate dynamics is expressed as
(x( v log(( + lo E ()))))))(())))
We approximate this by
M EM
exp( ()) -2)/ 21)((h-1) 5(m) ~ IN lz zz). (20)
kt
(2,,)
M Dm exp( (-f((m)))
Compared to the setting with only time discretization in Section A.3 the drift term 80 in one de- noising step the discretized dynamics 17) is replaced by the estimator (x (8 1)y) given
in (20). Using this substitution, the statement of Lemma 3, originally expressed as (66). is modified
accordingly as follows:
Hv, (Ma) S - -G2 ---11 --1=2) dl
1
Ju. (mL) = u² Exo.il [S (6 1)m) - 5=(2)]]? U, 2|]]
Auz
where ~ 4k+1/2.x and Alu ~M IN
x+11// 20)-
The expectation in the second term on the right-hand side is evaluated as:
Epo. [lo..[||...a (zo: (E - 1)m) $2 U, 2|u]
$ 21E (20) 3 + pm || (20, (e 1)n) 21 U, 2|u]
$ 21E us.. [1] [1]0)(20] - ==(x,)]1] + A [|| (e - 1)m) - 50(2))) [0]
+ [1 ( (e 1)n) N ( (e 1)n) 12 |U, U1]
(21)
which means we can independently evaluate the errors from time discretization (the first term), from finite N (the second term) and from finite M (the third term).
Score error from finite N. Let == h (E - 1)7, so that 2, coincides with 2p, the initial sample of the ê-th denoising step. For 2n E IRd, define the conditional law, where denoised 2 is conditioned on Z7,
x xx 27) x exp (xP(-)(3)-
|| 1
21
2T
with normalizer
7(z+) = / eexp (xxp(-1(3) 1|2- 2
dě.
2r
Let
G( (0) = Z(z+) /
1
IN (8 UU) exp (f(())
kt
11
2,
1/
2r
/
ØN
k+t
(# ( x x(3 27) dē.
Here assume:
20

---

Published as a conference paper at ICLR 2026
Assumption 3 (Bounded Covariance of the conditioned target distribution). For all 8 1, ,T and 27, there exists Cv > s.t. sup (aI,, Σ Cv.
Assumption (Lower bound of density ratio between the perturbed distributions).
π x 2-) * NN(0, (1))
N R RR0
@k+1/2(u Xx))
for a.e. u, hence G(ar |U) K.
Lemma 4 (finite-N eπror is O(1/N)). Under Assumption 3 and Assumption 4 the expectation of the evror term induced by the finite base sample size N satisfies
Exo du 1/a [1 -0 7) (2-r(2))11
2x 2xv
Pp ((122)) (N Il (xtzlxu) + x²(xx * *N Il a [xx)]
Nr²
In particular, under bounded diviergges, the second term on the right-hand side of the inequal- ity 21 is O(1/N).
Proof. By Tweedie's formula time T = h- (l-
1)17 we have
1
[ Egwnx x (-12-)18 27
Eng [ 7) -2++(2+))
Byo
E
72
(1((12>)
|1/
9N
k+
(2
{
-1 x xx($ { 2イ dš 2
72 Eyp
1/2
G(2r
(5)
Here used that IN (o =+) x
*+
IN A+ #(2 U) e) (-f(2)
18-2+11²) and substituted the 27
definition of G(z+ U) to obtain the normalized form.
Let E(z+) ʃ xx 2-) dé. Since J( (9N (EI0)
G(=,10) 1) xx(2 (+) dē 0, center the
integrand and apply Cauchy- Schwarz to obtain
Epp [[,[[11 [| 7) - sh 3h-r(2>))1 u]]
72 Ens (3-
IN (2
k+ A
- (2
U)
π xxx (r)) d3 ||]
G( UU)
k+ 1 (8 U)
IN
Byso Var V
Var (112=))
73
G(Zr U)
Cv Eyo / G(z,
(IN (2) G(zr U)) x xx(2 2+ de
kt
72
Uy2 15)
Cv
[ [ooNN (aN+ -(z U)) x (8 2. da
Epo
22
The second inequality is from Assumption 3 and Assumption 4 Ayields the last inequality. Using the unbiasedness identities with uy ~ 9k+1/2( Xx),
N(2;uj,,Ta)) 1, B xx 2+) * NO, hTa)
Eu Xx) (u+11/2(15 [xx]
1x+1/2(U, Xx)
and the decomposition ((g ((g-1) + (1 -G))² 2 1)2 + 2 1)2, we get
E Epo000441/4
[| h- -) - sh-r(2+)11
2x -Cv E Epo.gu+1/2
( ++
1 -)) UU) 1)²) xx(3) 2+) da
(2|U)
72
21

---

Published as a conference paper at ICLR 2026
Standard variance calculations for importance- weighted kernel mixtures yield
I
f + (oN (21 (2) 1)2 1 (2 2+) d
Bo x (+12,)
Var
N hua)
N
0041/21yy
E (0(2z) 1)²x* 2+) dē
= po x (#|=+)
Var
x iu 1/2
(241/2("g\xx)
Each variance upper bounded by the corresponding x² diverggee,,, yielding the claim.
Score estimation error from finite M. Here assume
Assumption 5. The following fourth moments under IN
k+ 11- 2r) are finite:
E [exp(A[(2))] < <c, [ ž, exp(e(((()))) < oo.
This can be satisfied, for instance, when the variance of each component in the Gaussian mixture 1+1/2 IN lu is sufficiently small compared to the growth rate of f. Lemma 5 (finite-M M error is O(1MM), Fix 2, and U ~NN Gk+ 4 (- Xx). Under Assumption S there exists a constant CM such that
CM
E Pa
[10
( h 7) 8N(2,, h 1eah---)11 1] 0 1
M
Proof. Write the estimator
20
at time T as a ratio of empirical means as follows:
M
1 M
Am
(2lm) - --(((m))) Bu = e [(3(m)).
CN, =
Am = M
-BM*
M
m=1
m=
Let A B[AM] and B E[BM] under g(m) Lid k+
9N
11 2r). Then & (ār, h - = A/(TB)
and
1
AM
A
2
AM - ABM
|| - 121
72
BM
B
72
By
Use the decomposition AM B
ABM
(AM
A)B
+ A(B
ABBr) obtain the crude bound
Am
1
A
2 211
All2
A| |Bas - B1
BM
B
B3
B2 B
Since e -((3) > O, Jensen implies B S M
En c Hence
||
1
M ||AM - A
2)(((m))
||
1
2)(2(m))
|l )
B1I
72
B
M
Taking conditional expectation with respect to 2/v ~ IN with exponents (, to each term, we have
1/2
l and applying Hölder's inequality
E |1 11 U]
< 7 (E(X))) 1 (BIIIAM - All")* +
|IA1
( - B>)) $}.
(B(BB)-
B²
where Xas = M EME e ())(2(m))
We compute the second moment of XM explicitly:
E
E|
1
E
Σ
(2((2()) (((cemm))
M
mamm'
M -
(BIE2 2
E
+
M
M
22

---

Published as a conference paper at ICLR 2026
To analyze the fourth moment of Σm= Um where E[um] 0, note that
2
M
Σ
Σ (Uss 1z))
(y) (Ck, vt).
=
Taking expectations, only the terms which each index appears an even number of times contribute, and hence
M
MB M(M (B 112) + 2M U2)
Um
=
m=
Thus, obtain
E AI - A
[I((2 2+)e -((a) El(E 2+)e (ee--(3)))1 11*)
3
Var e -1(3) (1)3)))
E
Ma
M2
E
[BM ( BY S M3
1
E (c- Ele--Eee----- Ele- ((3)))))))
3
Var (--(()) f(8)
M2
Substituting these bounds back gives
[[ - 1 72
2
(E(X21))
{
177 B[11 E [|
+
3
Var (2 - )e 5(3) )-))=-(3)) 2(3)) 2)
M
||||
1
E
(e (e-))) Ele (eee(())))))
1(3)
3
Var le 5(3) 13/
M3
M
If the fourth moments finite under Assumption 5 O(M-1). Therefore,
E O(1 and each rooted bracket is
CM M 7
E [I (2 r)1 U] <
for a constant CM depending only on the above moments.
Combining Lemma 4 and Lemma 5
conclude the overall score estimation error is 0(1/N 1/M) under several appropriate assumptions.
23

---

Published as a conference paper at ICLR 2026
B EXPERIMENT DETAILS
We provide parameters of each method for our experiments in Table and Table All experiments were conducted on an Intel Xcon CPU Max 9480 without GPU acceleration.
KL divergence estimation for the Gaussian Lasso experiment. Following the setting of Liang & Chen (2023b), we first ran the proximal sampler for 100,000 burn-in iterations, and then continued for 400,000 iterations. From this long trajectory we randomly collected 1000 samples to serve as reference particles. At each evaluation the experiment, the KL divergence was estimated using 1000 current particles and these 1000 reference particles. For estimation we used the k-nearest- neighbor estimator (Kozachenko Leonenko, 1987) implemented in Büth et al. (2025), with = recommended by Kraskov et al. (2004). Although this estimation involves sampling error, repeated the entire experiment with 10 different random seeds and, for each seed, computed the KL of particles. We then reported the mean and variance of
divergence based on the corresponding these estimates across the 10 runs to provide a more robust evaluation.
Table 1: Parameter setting for the experiment in Section 5.
Method
Parameters
Proximal Sampler with RGO
Initial distribution: (1) ~ (O)
v1/2
Step size: 7) 135
Number of independent chains: 100
Thinning: 10
Ours
Initial distribution: zo
~ ~N(00, Ta)
Step size: 1/10
Diffusion steps: = 10
Noise schedule: linear interpolation between 0 and h
Number of particles: N = 100
Number of interim samples: M 4000
Ours without interaction
Same Ours, except:
Number of particles: N = 1
Number of independent chains: 100
Table Parameter setting for the experiment in Section 5.2
Method
Parameters
In-and-Out
Initial distribution: TD _/ Ia) Step size: h =
Number of proposals for rejection sampling: 10000
(Particles are discarded if not accepted)
Number of independent chains: 1000
Ours
Initial distribution: xi NV
Step size: h =
Diffusion steps: T 10
Noise schedule: linear interpolation from 0.01 to 1
Number of particles: N = 1000
Number of interim samples: M 300
24

---

Published as a conference paper at ICLR 2026
h=1/20
h-10/135
1.2
h-1/10
h=2/10
1.0
KL divergence
0.8
0.4
0.2
0.0
150
200
300
350
lteration
Figure 5: Convergence KI. divergence for different step sizes h, with all other parameters set as Each curve shows the mean 10 random seeds, with shaded areas indicating
Ours in Table [ variances. Larger step sizes lead to faster convergence.
0.5
N=100, 4=0000
Exiiattd divergence
N00 M- 1000
0.4
N=1000, M:400
0.2
0.1
50
100
150
200
250
300
350
400
Iteration
Figure 6: Comparison of hyperparameter tuning while keeping the number of particles used during the algorithm's 's execution (i.e.. M N). For readability, the first iteration has been omitted. While increasing the number of particles N used for approximating the target distribution stabilizes the KI. divergence, changing both and results in very similar convergence patterns.
Effeet of step size. In addition to the default choice h = 1/10, we conducted experiments with other values of h (Figure 5}. We observe that larger step sizes accelerate convergence, which can be attributed to the heat flow more rapidly bridging the two modes. This phenomenon resembles diffusion-basee Monte Carlo methods, where pushforward dynamics from a Gaussian initialization cover the target distribution.
Hyperparameter sensitivity under fixed computational cost. Furthermore, Figure 6 shows that when the total number of particles used during the algorithm's 's execution (i.e., .M N) fixed, the algorithm exhibits comparable behavior across different choices of N and M. Increasing N im- proves the approximation of the target distribution, thereby stabilizing the estimated KI. divergence. Conversely, increasing M enhances the Monte Carlo estimation of the diffusion scores. Thus, under fixed computational budget, the method remains robust these hyperparameters.
25


Reviews:
[{"review_id": "Yl5nz1d0fm", "reviewer": "Reviewer_bjDQ", "summary": "This paper studies proximal sampling for a distribution satisfying the log-Sobolev inequality without implementing the restricted Gaussian oracle (RGO) by rejection sampling. Instead, it uses the idea of diffusion-based Monte Carlo to simulate the backward SDE (i.e., eq (5)) for the backward step in proximal sampling, that is, the RGO. By tracking $N$ particles, the proposed method approximates $\\pi^X$ by (9), and hence approximates the score function in (5) by (11). The key idea behind the approximation (9) is importance sampling. Thus, different from the exact implementation of RGO using first-order oracle of potential $f$ (for solving an optimization problem) and rejection sampling, Algorithm 1 approximately implements RGO by simulating (5) over $T$ internal steps within the stepsize $h$ of the (outer) proximal sampling.\n\nThe paper presents a complete analysis of Algorithm 1 with both discretization error and score estimation/Monte Carlo simulation error. It also provides two numerical experiments: 1) Gaussian-LASSO mixture; and 2) uniform sampling over non-convex sets. Numerical results in Section 5.1 demonstrate the advantages of the proposed method: it allows the use of a larger stepsize $h$ compared with proximal sampling based on rejection sampling, leading to faster convergence.", "questions": "1. The key idea behind the approximation (9) is importance sampling. It would be nice if this insight could be made explicit in the paper.\n\n2. Typo: line 710, a \"$y$\" is missing in the Gaussian.\n\n3. Line 130, Fisher information should be relative Fisher information.\n\n4. Some related papers might be missing:\nMixing Time of the Proximal Sampler in Relative Fisher Information via Strong Data Processing Inequality, Andre Wibisono\nProximal Oracles for Optimization and Sampling, Jiaming Liang and Yongxin Chen\nOracle-based Uniform Sampling from Convex Bodies, Thanh Dang and Jiaming Liang\n\n5. Similar to (Kook, Vempala, and Zhang, 2024), the last paper above by Dang and Liang studies proximal sampling for uniform sampling on convex bodies. It presents two implementations of RGO using projection and separation oracles. It would be interesting to compare with this paper as well in Section 5.2.", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Unlike the proximal sampling based on rejection sampling, this paper approximately implements RGO via simulation backward SDE and score function using a particle system (i.e., Monte Carlo simulation). The advantage (as reflected in the experiments) is that the stepsize $h$ of the proximal sampling could be taken large (as long as $T={\\cal O}(h)$). In particular, $h$ does not depend on dimension $d$ and properties of $f$ such as smoothness $L$. Moreover, the method does not require the first-order oracle of $f$, such as a subgradient; instead, it only assumes the zeroth-order oracle.", "weaknesses": "A major shortcoming is that the technical challenge addressed by the work is not very clear. Simulating SDEs with particle systems is a well-studied idea that has appeared frequently in the literature, so the methodological novelty seems limited, in view of the comparison with related works in Section 6. Moreover, the technical depth of the analysis is uncertain, as many of the proof techniques appear to be adapted from existing works, such as Vempala and Wibisono (2019).\n\nAnother notable shortcoming lies in the experiment in Section 5.2. In-and-Out finds T1 but fails to reach T2. It is claimed that the proposed method successfully explores both T1 and T2; however, based on Figure 4, it is difficult to conclude that the proposed method clearly identifies either region. The yellow points appear scattered all over, so there is no specific pattern that can be observed. The uniform sampling results are not convincing.\n\nFinally, the comparison only with RGO in Section 5.1 might be limited.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ZEVDVpx8TO", "reviewer": "Reviewer_zSYe", "summary": "The method provides an alternate implementation of the proximal sampler/RGO, based on an estimate for the scores using a particle ensemble of samples. This method yields a score estimator which appears as a mixture of Gaussians; in particular, the score estimator is implementable using zeroth order oracle queries. Theoretical analysis for this algorithm is provided, showing rates polynomial in the problem parameters. Experimental evidence then shows that this algorithm outperforms standard implementations of the RGO several benchmarks of interest.", "questions": "For this algorithm, when running the particle ensemble, do we take the entire ensemble of N particles as samples? This has the risk of inducing extra errors from correlation, but is unlikely to matter much in practice.\n\nI am also surprised that we only see $1/N$ in the error, compared to the usual curse of dimensionality for non-parametric estimation. What is the intuition here?\n\nLine 373 has some space formatting issues.\n\n\nNotation for R\\’enyi divergence switches between $\\mathcal R$ and $R$.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "The method provides a means for implementing the RGO in the proximal sampler using only \\emph{zeroth} order queries, which is a more general computational model than the standard gradient oracle model. Furthermore, it does not need any convoluted tricks (MALA + underdamped) to implement the proximal sampler, compared to prior theoretical proposals.\n\nThe method appears to work extremely well in practice when compared to the standard implementation of the RGO. Generally, I suppose it is not too surprising that these particle methods can perform well in practice, although theoretical guarantees may be harder to establish.\n\nThe method seems to be easily parallelizable.", "weaknesses": "The theoretical guarantees are not particularly strong; I would be surprised if in the LSI setting this could improve upon the guarantees for the usual implementation of the RGO. Indeed, as the error scales as 1/N in the number of particles, so we should expect $N \\asymp \\varepsilon^{-2}$ or polynomial in the accuracy (compared to the standard implementation). Of course, there are other drawbacks to the theory.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "VmmZjI8SKW", "reviewer": "Reviewer_rzzM", "summary": "This paper addresses a canonical isoperimetric sampling problem via a localized reverse heat flow and develops theoretical guarantees by leveraging the proximal-sampling analysis paradigm, reflecting the close connection between proximal operators and localized reverse diffusion. Unlike standard proximal samplers, the proposed method does not require first-order (gradient) information of the target distribution. This is enabled by a Tweedie-style identity that reformulates the denoising score, thereby removing the need for explicit gradients. While related identities have been used for importance weighting in [huang2024a], the probability path considered here is distinct—combining a global proximal path with a local reverse diffusion path—and the paper analyzes a Sequential Monte Carlo (SMC) scheme, which was not studied in [huang2024a].", "questions": "Refer to the above parts.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Although many components (proximal samplers, reverse diffusion Monte Carlo, and SMC) are drawn from prior work, the paper integrates them coherently to deliver a zero-order sampling scheme with provable convergence. The design is guided by a clear insight: Langevin-based methods inherently require gradients, whereas reverse-diffusion samplers can operate using only zero-order information.\n2. The paper explicitly elucidates the close connection between proximal sampling and reverse diffusion, an important and conceptually insightful observation that may inform the design of gradient-free samplers.\n3. The exposition is clear and well-structured: the intuition, assumptions, and main theorems are presented transparently, and the experiments are comprehensive, effectively demonstrating the proposed method’s efficacy.", "weaknesses": "1. The complexity analysis appears incomplete. The paper provides a one-step KL contraction bound with an additive term, but lacks an end-to-end error and cost analysis that composes these bounds over the full trajectory. A global, non-asymptotic complexity guarantee would substantially strengthen the contribution.\n2. The theoretical results are largely asymptotic. In practice, performance with finite $M$ (particles) and $N$ (time steps) is crucial. The assumptions used to control the one-step error via $M$ and $N$ seem idealized; we recommend validating their plausibility empirically on synthetic data and reporting sensitivity to $M$ and $N$.\n3. The comparison with RGO is unclear. What constitutes an ``iteration''? Does one pass that updates the entire sequence $\\{x_i\\}_{i=1}^N$ count as a single iteration or as $N$ iterations? Please clarify the accounting, and provide convergence comparisons against wall-clock time to enable fair, hardware-agnostic evaluation.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Rw2U2wTBPl", "reviewer": "Reviewer_RPkJ", "summary": "Method for non-parametric diffusion based approximation of proximal sampling exploiting zeroth-order information\nwith particle interaction. The authors use a Gaussian mixture surrogate for denoising the particles in the alternating proximal sampling process. They avoid the computation of gradient of $f$ for guidance in the SDEs and use of wasteful rejection sampling for denoising step. This leads to a multi particle sampling algorithm with deterministic runtime and exponential convergence guarantees with accurate score function esitmation", "questions": "The Monte Carlo estimation of the score function is increasing the overall computation with $M$ - But in the experiment section I didn't understand how the authors used large $M$ and compare to algorithms that don't have that many parallel threads in a fair setup. Can the authors provide a description of the fairness?\n\nMode coverage in the target is attributed to both the ability to use larger noise scale $h$ and the inverse weight of $x$ in high concentration regions. But the two seem to be in contention in $\\hat q_{k+1/2}(\\cdot|X_k)$, i.e close to $X_k$ and with small $h$ the component weights in Eq 10 will be smaller. Doesn't increasing $h$ at the same time make the weights more uniform?\n\nHave the authors studied the convergece time vs $M$ and $N$ when $N\\times M$ is fixed?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The authors present a very clear motivation and methodology for their work. The background is covered in a structure manner and the theoretical analysis and the experiments support the claims in the paper. The paper is very well written and the contributions are noteworthy. \n\nThe proposed sampling method uses simple surrogate to bypass inefficient use of rejection sampling and approximation of the score function. The method relies on noisy approximation of the score function of the denoising SDE, which derives it from zeroth order information and multi-particle interactions, improving the coverage of the disconnected modes of the distributions, and alleviating the need for projection and convex support.", "weaknesses": "The introduction stresses that the proximal algorithm should be scalable. However the experiments are limited to small $d$ scenarios. It would be informative to have additional results showing convergence time vs $d$ and compare the algorithm with the baselines.\n\nTypo: 146 $B_t^\\leftarrow$ is backward Brownian motion.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)