Chain of Log-Concave Markov Chains
Authors: Saeed Saremi
Ji Won Park
Francis Bach
Venue: ICLR 2024 poster (2024)
Keywords: Markov chain Monte Carlo
Langevin MCMC
log-concave distributions
kernel smoothing
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

Abstract: We introduce a theoretical framework for sampling from unnormalized densities based on a smoothing scheme that uses an isotropic Gaussian kernel with a single fixed noise scale. We prove one can decompose sampling from a density (minimal assumptions made on the density) into a sequence of sampling from log-concave conditional densities via accumulation of noisy measurements with equal noise levels. Our construction is unique in that it keeps track of a history of samples, making it non-Markovian as a whole, but it is lightweight algorithmically as the history only shows up in the form of a running empirical mean of samples. Our sampling algorithm generalizes walk-jump sampling (Saremi & Hyvärinen, 2019). The "walk" phase becomes a (non-Markovian) chain of (log-concave) Markov chains. The "jump" from the accumulated measurements is obtained by empirical Bayes. We study our sampling algorithm quantitatively using the 2-Wasserstein metric and compare it with various Langevin MCMC algorithms. We also report a remarkable capacity of our algorithm to "tunnel" between modes of a distribution.

Published as a conference paper at ICLR 2024
CHAIN OF LOG-CONCAVE MARKOV CHAINS
Saeed Saremi', Ji Won Park', Francis Bach²
I Frontier Research, Prescient Design, Genentech, South San Francisco, CA Inria, Ecole Normale Supérieure, Université PSL. Paris, France
ABSTRACT
We introduce a theoretical framework for sampling from unnormalized densities based on a smoothing scheme that uses an isotropic Gaussian kernel with single fixed noise scale. We prove one can decompose sampling from density (minimal assumptions made on the density) into sequence of sampling from log-concave conditional densities via accumulation of noisy measurements with equal noise levels. Our construction is unique in that it keeps track of history of samples, making it Mank-vann as whole, but it is lightweight algorithmically history only shows up in the form of a running empirical mean of samples. Our sampling algorithm generalizes walk-jump sampling (Saremi & Hyvärinen, 2019). The "walk" phase becomes a Markvvar)) chain of (log-concave) Markov chains. The "jump" from the accumulated measurements is obtained by empirical Bayes. We study our sampling algorithm quantitatively using the 2- 2-Wasserstein metric and compare it with various Langevin MCMC algorithms. We also report a remarkable capacity of our algorithm to "tunnel" between modes of a distribution.
the
INTRODUCTION
Markov chain Monte Carlo (MCMC) is an important class of general-purpose algorithms for sampling from an unnormalized probability density of the form p(z) = e f(=)) /Z IR.. This a fundamental problem and appears in variety of fields, e.g., statistical physics going back 1953 (Metropolis et al., 1953), Bayesian inference (Neal, 1995), and molecular dynamics simulations (Leimkuhler & Matthews, 2015). The biggest challenge facing MCMC is that the distributions of interest lie in very high dimensions and are far from being log- log-concave, therefore the probability mass is concentrated in small pockets separated by vast empty spaces. These large regions with small probability mass make navigating the space using Markov chains very slow. The second important challenge facing MCMC that the log-concave pockets themselves are typically ill-conditioned highly elongated, spanning different directions for different pockets-which only adds to the complexity of sampling.
The framework we develop this paper aims at addressing these problems. The general philosophy here is that of smoothing, by which we expand the space from Rdd to Rm for some integer m and "fill up" the empty space iteratively with probability mass in an approximately isotropic manner, the degree of which we can control using single smoothing (noise) hyperparameter The map from (noisy samples in) 1nd back to (clean samples in) Rd is based on the empirical Bayes formalism (Robbins, 1956; Miyasawa, 1961; Saremi & Hyvärinen, 2019; Saremi & Srivastava, 2022). In essence, single "jump" using the empirical Bayes estimator removes the masses that were created during sampling. We prove general result that, for any large m. the problem of sampling in Rd can be reduced to sampling from a sequence of log-concave densities: once log-concave, always log-concave. The trade-off here is the linear time cost of accumulating noisy measurements over m iterations.
More formally, instead of sampling from p(x). we sample from the density p(y1m) that is associated with Y1.m (Yi...., Y,), where Y. X N₁, t E [m]. N, X, and N, ~ N(O, o² I) all independent. As we show in the paper, there is a duality between sampling from p(x) and sampling from p(y1 ) in the regime where m 1 is small, irrespective of how large is. This is related to the notion of universality class underlying the smoothed densities. Crucial our formalism is keeping track of the history of all the noisy samples generated along the way using the factorization
p(y = p(v) II P(yelvie-1).
(1.1)

---

Published as a conference paper at ICLR 2024
eee
vi mm)
1/m + (mhe)
91
97
Figure 1: Chain of log-concave Markov chains. Here, (yi ( (e(n.]] are samples from a Markov chain, which is used generate independent draws from p(ye|vi1t-11 for E [m-]. The blue arrows indicate the non-Markov aspect of our sampling scheme: the accumulation of noisy measurements. The wiggly arrows indicate the denoising "jumps". In this example, p(ye (1-r-1) is log-concave for all t, but the jumps asymptotically sample the target density (a mixture of two Gaussians) as increases.
An important element of this sampling scheme is therefore non-Markovian. However, related to our universality results, this history only needs be tracked in the form of an empirical mean, so the memory footprint is minimal from an algorithmic perspective. See Fig. for schematic.
more technical summary our contributions and the outline of the paper are follows:
In Sec. 2, we prove universality results underlying the smoothed densities p(y1.m)-
We study anisotropic Gaussians in Sec. 3, proving a negative resuli regarding the condition number of p(y1.m) in comparison p(y1) in the becomes a segue to our factorization (1.1), where in remarkable contrast we show that the condition number monotonically improves upon accumulation of measurements.
universality class. This analysis
Sec. 4 is at the heart of the paper, where we prove several results culminating in Theorem 1, which shows that broad class of sampling problems be transformed into sequence of sampling from strongly log-concave distributions using our measurement accumulation scheme. (This is feasibility result; particular, we do not prove here that the log-concave sampling strategy is optimal.) We examine the theorem by algebraically studying an example of a mixture of Gaussians in detail. In Sec, 4, we also outline our general sampling algorithm.
We validate our algorithm on carefully designed test densities in Sec. 5. In particular, our algorithm results lower 2-Wasserstein metric compared sampling from p(x) using Langevin MCMC (without any smoothing). We also qualitatively report the capacity of our log- log-coneave sampling scheme to tunnel to mode distribution with small probability mass in small number of steps when it is initialized at mode with much higher mass.
1.1 RELATED WORK
Our solution, sketched above, has roots in walk-jump sampling (Saremi Hyvärinen, 2019) and its recent generalization (Saremi & Srivastava, 2022). Both papers were framed within the context of generative modeling, i.c., sampling from an unknown distribution from which one has access to independent samples. In contrast, this work lays the theoretical foundation for the fundamental problem of sampling from an unnormalized density when there are no samples available. addition, regarding the recent development, Srivastava (2022) regarding the distribution p(y).)) being well-conditioned not correct. This nontrivial negative result motivates our analysis of the non-Markovian scheme for sampling p(y)..)-
show analytically that the intuition expressed by Saremi &
2

---

Published as a conference paper at ICLR 2024
Our methodology is agnostic to the algorithm used for sampling from p(gelvra-1) in (1.1). However, we have been particularly motivated by the research on Langevin MCMC which is a class of gradient- based sampling algorithms obtained by discretizing the Langevin diffusion (Parisi, 1981). There a growing body of work on the analysis of Langevin MCMC algorithms of various complexity (overdampled, Metropolis-adjuste distributions (Dalalyan, 2017; Durmus & Moulines, 2017; Cheng et al., 2018; Dwivedi et al., 2018; Shen & Lee, 2019; Cao et al., 2021; Mou al., 2021; Li et al., 2022).
underdamped, higher-order) for sampling from log-concave
There is a significant body of work on sequential methods for sampling, rooted in annealing methods in optimization (Kirkpatrick al., 1983), which became popular in the MCMC literature due Neal's seminal paper on annealed importance sampling (Neal, 2001). Diffusion models -ick-tic al., 2015; Ho al., 2020) are a related class of sequential methods for generative modeling. Our sequential scheme is distinguished from earlier methods from separate angles:
(I) Although we sample the conditional densities with Markov chains, we condition on all the previous samples that were generated. As whole, our scheme is strictly non-Markovian. (II) In our sequential scheme, we are able to guarantee that we sample from (progressively more) log-concave densities. To our knowledge, no other sampling frameworks can make such guarantees. (III) Compared to diffusion models, the noise level in our framework is held fixed. This is an important feature of our sampling algorithm and it underlies many of its theoretical properties. (IV) All prior sequential schemes rely on a noising/annealing schedule which is hard to tune, and their performance is sensitive to the choice of the schedule (Karras al., 2022; Syed et al.. 2022). In contrast, our sequential scheme is free of scheduling and relies on only two parameters: the noise level a and the number of measurements m.
Notation. We use p to denote probability density functions and adopt the convention where we drop the random variable subscript to when the arguments are present, e.g., p(z) = Px (æ), p(y21 = P (32). We reserve f to be the energy function associated with p(z) x e f(x), We use to denote the spectrum of a matrix, e.g., Amax (C') is the largest eigenvalue of C. We use the shorthand notations [m]] (.....,,) Vl.m (y),,... Mm), and Ii m -1 y/.
2 UNIVERSAL (o, m) -DENSITIES
Consider the multimeasurement (factorial kernel) generalization of the kernel density by Saremi & Srivastava (2022) for m isotropic Gaussian kernels with equal noise level (kernel bandwidth) a:
p(yim) x [ e- (1(()) exp
x well )ell")dz.. dzz
(2.1)
gd
20
We refer to p(y1:m) the (o, densiyy Equivalently, Yi/l lid N(I,021) t [m]. Clearly, p(y1m) is permutation invariant p(y1..... Fnm = (((y(())..... Mr(m)), where π [m] + [m] is a permutation of the measurements. We set the stage for the remainder of the paper with calculation that shows the permutation invariance takes the following form (see Appendix A):
logi()) (TTimii tami m 11)
(17
+ cst,
(2.2)
202
m
where
p(y;o) = log / e-5(m) exp
(
1-22211
1
wll d
(2.3)
The calculation is straightforward by grouping the sums of squares in (2.1):
Σ
I I
voll m (-Ilz Vi-ml 11? + |7
leall?) + cst.
In addition, the Bayes estimator of X given Yi.m = Fi:m simplifies as follows (see Appendix A): E[X|y1=m] F1.m + i V m. m
(2.4)
These calculations bring out a notion of universality class that is associated with p(y):m) formalized by the following definition and proposition.
Our method can be viewed a discretization scheme in the stochastic localization method, which we plan formalize in future research.
3

---

Published as a conference paper at ICLR 2024
Definition (Universality Class). We define the universality class [ō] the set of densities p(y1:m)), in the family of (a, -endsiiiessss such that for all (a.m) [ō] the following holds: m-1/2g Ø. Proposition If fm ~ p(ymm)) let Pa.m be the distribution of ( and define pa = Prr1-- Then Pa.m- = Pm 1/2g. P In other words, Pa.m is identical for densities in the same universality class.
Proof. We are given X e-((x) Y, X + Et E N(0, o20, independently for 1 E [m]] follows Iim- X E, where & ~ N(0,02T), where o² mm- Using (2.4), B[X [y:m is
distributed as
x x + E;
which is identical for all densities p (().m) [ō].
2.1 DISTRIBUTION OF E[X VS. px: UPPER BOUND ON THE 2- WASSERSTEIN DISTANCE
Our goal is obtain samples from px. but in walk-jump sampling the samples are given by E[X 11:mm] where 1/1m ~ p(ymm)) (Saremi & Hyvärinen, 2019; Saremi & Srivastava, 2022). Next, address how far Po.m is from the density of interest px.
Proposition 2. The squared 2- Wasserstein distance between px and pa.m is bounded by
W2(xx.Pa.N
d.
The proof given in Appendix As expected, the upper bound is expressed in terms of o² o² /2/m. m.. A close inspection of the proof shows that the bound above is loose as it is obtained from the rate resulting from replacing the empirical Bayes estimator [[X]y with the empirical mean Y 1am- Note, however, that when the prior p(z) is "strong" (e.g., low entropy), the dependence on o² m can be significantly improved.
THE GEOMETRY OF (a, m) DEENSTIES
In this section we analyze at the problem of sampling from p(Vi.)) where consider p(x) to be an anisotropic Gaussian, ~ N(0, C)2, with a diagonal covariance matrix:
C' diag(?f...... -2).
(3.1)
The density px is strongly log-concave with the property T I f2f(x) < T therefore its condition number is --2 2
"max- Log-coneave densities with > 1 considered ill-conditioned. Since Y1 ~ N(0,, C + o² 1). the condition number for (single-measuremmen) smoothed density, which we denote by Ka,1 is given by:
Kig, = 1011(11- -1+00010-22 maxx )/
(3.2)
Next, we give the full spectrum of the precision matrix associated with (a, m) mm)dennities.
Proposition 3. Consider anisotropic Gaussian density X ~ N(O, C) in R.d, where Cij 78grj- Then the (o, me-diiiity is centered Gaussian in Rmnd Y1m N(0, F.Jh)- For m > 2, the precision matrix Fo is block diagonal with d blocks (indexed by i) of size m X m, each with the following spectrum: (i) There are m - 1 degenerate eigenvalues equal to o -2 (ii) The remaining eigenvalue equals to (o² + mm?)) The condition number Ka,m associated with the (a, m--density is given by:
Amaxx (Fam
Ker.m Xmni (Fa.m)
mmax
Remark (The curse of sampling all measurements at once). The above proposition is a negative result regarding sampling from p(y1:m) if this is important "if"-al measurements y1om are sampled in parallel (at the same time). This is because mo -2 = 5 remains constant for m > 1 for (a, m) [ō [ve-even worse, the condition number Ka,m is strictly greater than K5, for m > 1.
This negative result regarding the sampling scheme by Saremi & Srivastava (2022), we call joint multimneasurement sampling (JMS), leads to investigation below into sampling from p(y1:m) sequentially using the factorization (1.1) that we call sequential multimeasurement sampling (SMS). Now, we perform the analysis in Proposition 3 for the spectrum of the conditional densities in (1.1).
4

---

Published as a conference paper at ICLR 2024
Proposition 4. Assune X ~ N(0, C) is the anisotropic Gaussian in Proposition 3. Given the factorization of ply1m) ) in (1.1). t 1, the conditional density p(y (1:t-1) is a Gaussian with a shifted mean, and with diagonal covariance matrix:
Asi
[-1
2 +
-2 pog
Σ
Ati)
Uti
Σ yki
est,
Ati
where Aui is short for Asi (1+ T²T The precision matrix associated with p(ye ((3 1). denoted by Fij1:t-1, has the following spectrum
(F =1 o';
with the following condition number
1 (t + 0 2
"min
1 1- (t + o22TTixx
Lastly, the condition number Kajia-1 1 is monotonically decreasing as t increases (for any m >> 1):
1 < Km[1:m-1 < .... < K3/1:2 < 1221 <
(3.3)
where = Ko,1 is given by (3.2).
The proofs for Proposition 3 and Proposition 4 are given in Appendix C. These two propositions stand in clear contrast to each other: in the SMS setting of Proposition 4, sampling becomes easier by increasing t as one goes through accumulating measurements Vi-4 sequentially, where in addition K1 can itself be decreased by increasing Next, analyze the SMS scheme in more general settings.
CHAIN OF LOG-CONCAVE MARKOV CHAINS
Can we devise sampling scheme where we are guaranteed always sample log-concave densities? This section is devoted to several results in that direction. We start with the following two lemmas.
Lemma 1. Assume Vx E Idd, 4²f(x) x LI and [[)]l > A|| zoll
A for some zo. Then,
Vy E Rdd
(loog)) - (-1+ n202
3Ld
3A"
+3 | m113) I
p2a²
1206
0²
The proof is given in Appendix D.
Lemma' 2. Consider the density p(x) associated with the random variable X in Rd and the m)) density given by (2.1). Then in expectation, for any m 2 the conditional densities become more log-concave upon accumulation of measurements."
Eyg 2 logp(yn) > Ema V2 log logp(y22u) (()2))))) > ... b>= Eyy.. V2 Mmm
plym //1mmm1)) 1)
Proof. The full proof of the lemma is given in Appendix D, where we derive the following:
2 Nm log p(ym 11) = - -T + o cov X [y1:m)
The proof follows through since due to the law of total covariance the mean the posterior covariance Ey.... cov (X (1:m) can only go down upon accumulation of measurements.
These two lemmas paint intuitive picture that expand in the remainder of this section: (i) by increasing a we can transform density to be strongly log-concave (Lemma 1) which we can expect sampling
sample first measurement from, (ii) and by accumulation of measurements to become casier, where in Lemma this is formalized by showing that on average the conditional densities become more log-concave by conditioning on previous measurements. Next, we generalize these results with our main theorem, followed by an example on mixture of Gaussians.
Note that here assumption is made the smoothness of p(z).
5

---

Published as a conference paper at ICLR 2024
logo(stoon
0.6
322
23
Pim
Vim
10
(a) 0 =
(b) 0 2/2
(c) (m() vs. m
Figure 2: (a.b) The negative conditional Hessian for two values of are plotted a function of 7 and m assuming X is distributed according to (4.4) in ID, where we set = 3, T (see (4.5)). (c) The upper bound in (4.1) sharp for this example; o² ((m) is plotted vs. m for different
Theorem 1 (Once log-concave, always log-concave). Consider to be a random variable Rd4 with a compact support, i.e., almost surely || 2 & R², and take X Z No, No Z. No N(00721) Then, for any m 2 1, the conditional Hessian is upper bounded
V2 log (mmmm//mmmmm))) (1.m-1) $ {(m)I,
(4.1)
Um
where:
R²
((m)
- 1) + (m7 02ح²
(4.2)
T2 m1? + 02
is a decreasing function of m. in particular:
<'(m) = 72(2R202 + 0272 m
E 0.
0² + mT²)3
As a corollary, p(y1) associated with Y1 x N2, Ni N(O, I² is strongly log-concave if
2²> -
(4.3)
and stays strongly log-concave upon accumulation of measurements.
Proof. The full proof is given in Appendix D and it is direct consequence of the following identity:
1 0² mr² + 0 11) I (mr² + ozح² cov Z (1)ms
V2 Vm log (omll(mmm)
which we derive, combined with cov(Z (1mmm) & R2 due to our compactness assumption.
Remark Theorem I spans a broad class of sampling problems, especially since can in principle be set to zero. The only property we loose in the setting of 0 is that the upper bound ((m) does not monotonically go down as measurements are accumulated.
4.1 EXAMPLEE: MIXTURE OF TWO GAUSSIANS
In this section we examine Theorem by studying the following mixture of Gaussians for 1/2:
p(x) N(x; p1, 7 (1 a1) -1,
(4.4)
This is an instance of the setup in Theorem 1, where p(z) = j=((z - 11) + 3 6(z + n). and R² = p1
By differentiating (2.2) twice we arrive at the following expression for V 9m log pl Um V1.mm
Vm log p 11) 2 ym log 1ogp(y)m) = 2m 1)1 + m
(Timin m 1/2 1/2g), (4.5)
where H(y:))= = (P see (2.3) for the definition of In Appendix D we show that for the mixture of Gaussian here, (4.4) with a 1/2, we have
2p4
(1 (1+
+1cchhh cosh 24 22TTy y @ 1 72
))
H(y;o) = +22
- -T
(4.6)
² + 72

---

Published as a conference paper at ICLR 2024
which takes its maximum y 0. By using (4.5), it is then straightforward to show that (4.1), (4.2), and (4.3) all hold in this example, with the additional result that the upper bound is now tight. In Fig. 2, these calculations
visualized in ID for μ = 3, T 1, and for different values of a; in panel (c) we also plot 1 1/m which is the large m behavior of o² o²{(m). This can be seen from two different routes: (4.2) and (4.5).
Remark 3 (Monotonicity). The monotonic decrease of the upper bound in Theorem 1, together with the monotonicity result in Lemma may lead one to investigate whether the stronger result
Vy logp(s)) & V log pogg((22111 (2) ..... b V log (1l could hold, e.g., for the mixture of Gaussians is sharp for this example. For (4.7) to hold, cov (Z (Z(ymm) x cov (Z(y)mm-1) almost surely. However, we can imagine a scenario where y/ + ym- is very large, so that cov (/11mm)) is such that y1 +... ym-1 + ym is close to mIE where cov( (Z)V1m) will be large.
(4.7)
studied here, especially since the upper bound (4.1)
while you
4.2 ALGORITHM: NON- MARKOVIAN CHAIN OF (LOG-CONCAVE) MARKOV CHAINS
Below, we give the pseudo-code for our sampling algorithm. the inner loop, MCMC, is any MCMC method, but our focus in this paper is on Langevin MCMC algorithms' that use Vy, p(g (11-1)) to sample the new measurement Y, conditioned on the previously sampled ones V1:--1
Algorithm 1: Sequential multimeasurement walk-jump sampling referred to by SMS. See Fig. 1 for the schematic. A version of MCMC, is given in Appendix E.
1: Parameter noise level a
2: Input number of measurements m, number of steps for each measurement n4
3: Output X
4: Initialize Y = 0
5: for fort= prrt=[1..... [I.....,. m| do
6: Initialize y-
7: for = [1.....,............ n4] do
y = MCMC, (Y, (i-1) Y 11-11)
8:
9:
end for
10:
Y = --Y((m)
11:
Y1.- = -Y Y + (Yr - 11=-1)/t
end for
13: return
X -- V according to (2.4)
4.2.1 ESTIMATING V log p(y)
So far we have assumed we know the smoothed score function g(y; == ((log)) - o). and in experiments below we consider cases where we know g o) in closed form. In general settings, we would like to estimate g in terms of the unnormalized p(x) x e f(x), Given (2.3) and (2.4), we write a an expression for g(y; o) which can be turned into various estimators:
g(y: o) = 22 [E[ - v), X ~ e f-iiio where f(x;u, f(x) 20 | w112. (4.8)
See Appendix F, where we also give two estimators for depending on how E[X] in (4.8) is estimated: Iplugin (F.1) is plug-in estimator obtained by importance sampling, Ilangevin (F.2) is obtained by
sampling X using Langevin MCMC. Finally, use (2.2) to write V log p(yo 11) in terms of g:
Vn log ponlpllaaaaaa V, logr(u):) = 7
) -
a
(Ui+ - -y)-
(4.9)
The above expression is used in the inner loop of Algorithm I sample from p(y, W1.t-1) sequentially, We conduct experiments to investigate this aspect of the problem (replacing g with g) in the appendix.
'we experimented with variety of Langevin MCMC algorithms to sample from p(ye (y:=-1) in the inner loop of Algorithm 1. The results are reported in the appendix due to space constraints. After extensive tuning, we found the algorithm by Sachs al. (2017) be the hest performing for the test densities considered.
7

---

Published as a conference paper at ICLR 2024
Elliptical Gaussian
Mixture Gaussioos
Eliptical Gaustian
Mixtre Ganesians
10.
10
600
Logg
(a)
(b)
(c)
(d)
Figure 3: (a, b) Sliced 2- Wasseessteinn distance vs. and (c, d) Sliced 2- W-sassesteein distance vs. a
for varying d for the elliptical Gaussian and Gaussian mixture target densities using SMS.
EXPERIMENTS
We evaluate the performance of Algorithm alongside related sampling schemes on carefully designed test densities. We compare the following sampling schemes:
Sequential multimeasuremen walk-jump sampling (SMS), Algorithm with m = 1000,
Joint multimeasurement walk-jump sampling (JMS),
Single-measurement walk-jump sampling ("m = 1*),
Langevin MCMC by Sachs al. (2017) without any smoothing ("o 0"),
Exact samples from pa.m ("Oracle").*
Metric. We use the sliced 2- Wasserstein metric (Bonneel et al., 2015; Peyré & Cuturi, 2019) to quantify the consistency of the obtained samples with the target density px- We use 1,000 projection directions drawn from the Gaussian distribution (Nadjahi et al., 2021, Eq.9).
MCMC algorithms. For all the results in this section, implement MCMC sampling based on underdamped Langevin diffusion (ULD). The particular algorithm used for the results shown this section extends the BAOAB integration scheme using multiple time steps for the O-part (Sachs et al., 2017). In Appendix H. we present the full comparison across other MCMC algorithms, including other recent ULD variants (Cheng al., 2018: Shen & Lee, 2019) well as the Metropolis-adjusted Langevin algorithm (MALA) (Roberts & Tweedie, 1996; Dwivedi et al., 2018).
Score estimation. In Appendix we compare sampling with the analytic score function, the plug-in estimator of the score function given in (F.1) with varying numbers of MC samples n, and the Langevin estimator of the score function given in (F.2).
Hyperparameter search. The hyperparameters were tuned for each sampling scheme and the total number of iterations was kept fixed in our comparisons reported here. See Appendix G for details.
5.1 ELLIPTICAL GAUSSIAN
The elliptical Gaussian features a poorly conditioned covariance: X ~ N(O, C),, where we set 7 in (3.1). For each d, the noise level a and other hyperparameters of the sampling algorithm, such as step size and friction, were tuned. Fig. 3(a) plots the 2-Wasserstein distance with varying d. Our main observation here that SMS outperforms JMS, which is expected from our theoretical analysis in Sec. 3. In addition, as Fig. 3(c) shows. SMS was robust to the choice of a particularly for large d in our experiments. Finally, the underdamped Langevin MCMC by Sachs et al. (2017) does quite well when friction is carefully tuned; we would like to investigate this for larger condition numbers in future research.
7f 0.1,
'Oracle is short for the SMS oracle. For the test densities considered here we can draw exact samples from fo.m m. the distribution of V..]. This haseline is used to separate the issue of the closeness of fa.m Px from the problem of sampling f..m itself.
8

---

Published as a conference paper at ICLR 2024
Intiibaation
hitidbat
6
15442100
(a) SMS, 3 trajectories
Langevin MCMC
(c) SMS, 100 samples
(d) Langevin. 100 samples
Figure 4: Tunneling phenomenon. Trajectories of three walkers (a) under our SMS sampling scheme using Langevin MCMC by Sachs al. (2017) the inner loop of Algorithm (b) the same Langevin MCMC without any smoothing. Purple and green indicate the beginning and end of trajectories, respectively. (c, d) The final samples for 100 random trajectories (using identical seeds). Fewer samples reach the smaller mode in Langevin MCMC compared to SMS.
5.2 MIXTURE OF GAUSSIANS
To evaluate mixing of multiple modes, we consider the test density (4.4) with = 1/5, 7 1, and μ 3. La. where la is the d-dimensional vector (I....,1) 1) As Fig. 3(a) shows, SMS achieves consistently low (sliced) 2-Wasserstein distance with increasing d, whereas other sampling schemes deteriorate in performance. We observe, in Fig. 3(b), that SMS outperforms the best-per rest-performing underdamped Langevin MCMC (a 0) in experiments for at least one value for all d. Higher requires larger a. In addition, we would like to highlight the following:
(I) Vanilla walk-jump sampling (m 1) highly ineffective the dimensions increase. This is in contrast to the sampling from anisotropic Gaussian (already log-concave) in Sec. 5.1. (II) The optimal a here is in fact larger than the noise level needed to make p(y)) log-concave. This related the benefits of sampling from better-conde cetter-oonditioed log-con log-concave distributions, which is well-known in the literature. In Appendix J, we include results for mixture of correlated Gaussians supporting the same conclusion.
5.3 TUNNELING PHENOMENA
Fig. 4 illustrates the trajectories of three walkers (a) under our SMS sampling scheme and (b) using Langevin MCMC. Each walker has the same random seed between (a) and (b) and was initialized at (3, 3), the dominant mode with 80% of the mass. With SMS, walker is able tunnel to the smaller mode fairly quickly, whereas for Langevin MCMC (without smoothing) all three walkers are stuck around the dominant mode. In panels (c) and (d) also show the histogram of final samples in the same setup with initialization at (3, 3) for 100 walkers after 100 K steps.
In summary, we have consistently observed that the same (Langevin) MCMC algorithm, when used in the inner loop of Algorithm 1, is more effective than when used without Gaussian smoothing.
6 CONCLUSION
In this paper, we established a theoretical framework that reduces the general problem of sampling from unnormalized distribution that of log-concave sampling defined by single noise parameter. We conclude with two main limitations of this work at the present time: (I) Our results do not make clear if the log-concave sampling strategy is optimal. The issue of "optimality" is challenging as it is inherently problem-dependent and additionally depends on the MCMC algorithm used in the inner loop of Algorithm 1. (II) Related to the issue of optimality is the fact that for general sampling problems, the smoothed score functions required to sample from p(ya (1tt111:1))) f E [m] using gradient-based methods need to be estimated. This is a complex problem and should be investigated in future research. Finally, immediate application of the machinery we developed here is the problem of generative modeling, as the m smoothed score functions needed in running Algorithm can indeed be learned (approximated) using empirical Bayes least-squares denoising objectives. This approach similar to training diffusion models; however, our sampling scheme fundamentally different, as it relies on the accumulation of measurements, controlled by single noise parameter,

---

Published as a conference paper at ICLR 2024
REFERENCES
Nicolas Bonneel, Julien Rabin, Gabriel Peyré, and Hanspeter Pfister. Sliced and radon Wasserstein barycenters of measures. Journal of Mathematical Imaging and Vision, 51:22-45, 2015.
Yu Cao, Jianfeng Lu, and Lihan Wang. Complexity of randomized algorithms for underdamped Langevin dynamics. Communications in Mathemutical Sciences, 19(7): 2021.
Sourav Chatterjee and Persi Diaconis. The sample size required in importance sampling. The Annals of Applied Probubility, 18951109 2018.
Pratik Chaudhari, Anna Choromanska, Stefano Soatto, Yann LeCun, Carlo Baldassi, Christian Borgs, Jennifer Chayes, Levent Sagun, and Riecardo Zecchina. Entropy-SGD. Biasing gradient descent into wide valleys. In International Conference on Learning Representations, 2017.
Xiang Cheng, Niladri S. Chatterji, Peter L. Bartlett, and Michael I. Jordan. Underdamped Langevin MCMC: non-asymptotic analysis. In Conference on Learning Theory, pp. 300-323, 2018.
Arnak Dalalyan. Theoretical guarantees for approximate sampling from smooth and log-concave densities. Journal of the Royal Statistical Society. Series B (Statistical Methodology), pp. 651-676, 2017.
Alain Durmus and Éric Moulines. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. The Annals of Applied Probability, 27(3):1551 1587, 2017.
Raaz Dwivedi, Yuansi Chen, Martin J. Wainwright, and Bin Yu. Log-concave sampling: Metropolis- Hastings algorithms are fast! In Conference on Learning Theory, pp. 793-797, 2018.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems, 33:6840-6851, 2020.
Aapo Hyvärinen. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(Apr):695-709, 2005.
Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine. Elucidating the design space of diffusion- based generative models. Advances in Neural Information Processing Systems, 35:26565-26577, 2022.
Scott Kirkpatrick, C. Daniel Gelatt Jr., and Mario P. Vecchi. Optimization by simulated annealing. Science, 220(4598):671-688, 1983.
Ben Leimkuhler and Charles Matthews. Molecular dynamics: With deterministic and stochastic
numerical methods. Interdisciplinary Applied Mathematics, 39:443, 2015.
Ruilin Li, Hongyuan Zha, and Molei Tao. Sqrt(d) dimension dependence of Langevin Monte Carlo. In International Conference on Learning Representations, 2022.
Nicholas Metropolis, Arianna W. Rosenbluth. Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller. Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6))1087-1092, 1953.
Koichi Miyasawa. An empirical Bayes estimator of the mean of normal population. Bulletin the International Statistical Institute, 38(4): 1961.
Wenlong Mou, Yi-An Ma, Martin J. Wainwright, Peter L. Bartlett, and Michael I. Jordan. High- order Langevin diffusion yields an accelerated MCMC algorithm. Journal of Machine Learning Research, 1911:1911-1959 2021.
Kimia Nadjahi, Alain Durmus, Pierre Jacob, Roland Badeau, and Umut Simsekli. Fast approxima- tion of the wasserteess distance using concentration random projections. Advances in Neural Information Processing Systems, 14411-22444 2021.
Radford M. Neal. Bayesian learning for neural networks. 1995. PhD thesis, University of Toronto.
Radford M. Neal. Annealed importance sampling. Statistics and Computing, 11:125-139, 2001.
10

---

Published as a conference paper at ICLR 2024
Giorgio Parisi. Correlation functions and computer simulations, Nuclear Physics B, 180(3)3378-38) 1981.
Gabriel Peyré and Marco Cuturi. Computational optimal transport: With applications to data science. Foundations and Trends Machine Learning, 1165-6)355-607, 2019.
Herbert Robbins. An empirical Bayes approach to statistics. In Proc. Third Berkeley Symp., volume 1, pp. 157-163, 1956.
Gareth O. Roberts and Richard L. Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, pp. 341-363, 1996.
Matthias Sachs, Benedict Leimkuhler, and Vincent Danos. Langevin dynamics with variable coeffi- cients and nonconservative forces: from stationary states to numerical methods. Entropy, 19(12): 2017.
Saeed Saremi and Aapo Hyvärinen. Neural empirical Bayes. Journal of Machine Learning Research, 20(181):1122) 2019.
Saeed Saremi and Rupesh Kumar Srivastava. Multimeasurement generative models. International Conference Learning Representations, 2022.
Ruoqi Shen and Yin Tat Lee. The randomized midpoint method for log-concave sampling. Advances in Neural Information Processing Systems, 32, 2019.
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. International Conference on Machine Learning, pp. 2052222666 2015.
Saifuddin Syed, Alexandre Bouchard-Côde George Deligiannidis, and Arnaud Doucet. Non- reversible parallel tempering. scalable highly parallel MCMC scheme. Journal of the Royal Statistical Society Series B: Statistical Methodology, 84(2):321-350 2022.
Cédric Villani. Topics in Optimal Transportation, volume 58. American Mathematical Society, 2021.
Martin J. Wainwright and Michael I. Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning, 1(1-2):1-305 2008.

---

Published as a conference paper at ICLR 2024
A DERIVATIONS FOR SEC. 2
A.1 DERIVATION FOR (2.2)
In this section we derive a general expression for the (a, m)--ssiities p(3imm), short for po (3/mm)
p(y):m) = / p(x)
dr 1=1
jud
|pc 2
20 || - voll²) dx
exp
(2 yd/2
/
e-((x)
exp
x yall )d.
de.
Rd
2
The (a, m)ddenittyy p()..m) is permutation invariant under the permutation of measurement indices:
[m]
[m]]
In the calculation below derive general form for p()1.m) where this permutation invariance becomes apparent in terms of the empirical mean of the m measurements
1am
m
and the empirical mean of {llb || We start with rewriting of log p(yim//z):
m
-2 log p(y)mla))
11ye x=|2 + est
=m||2:|2 2:
Yc, x)) +
||
est
1=1
-m(1I= Fi.mll 11
|
11?
+
llval())
+ est,
where est is constant that does not depend yiam- Using the above expression, arrive at:
log p(y1:m) log / (=1()) (x) dr + est
log
e)(()) exp(
m
Fimm 1²)de
17:m 1 m
> TM
1=
||
+ cst.
2o
2m 102
The equation above reduces to (2.2) with the following definition
()) log / E (f0)) exp
20 |1 y/122) wll2)dtr. dx.
A.2 DERIVATION FOR (2.4)
Next, derive the expression for (y1..) E[ [1:mm given in (2.4):
E E[
= yt + V²V, log p(y1:m)
1/1 + o² m V -1/2 + o-----2000 -2 Dh))
(A.1)
V -1/2gr)
Vim + m
The first equation above comes from the generalization of the Bayes estimator to factorial ker- nels (Saremi & Srivastava, 2022), and for the second equation used (2.2).
12

---

Published as a conference paper at ICLR 2024
B SMOOTHED SCORE FUNCTIONS
In this section, we give several different expressions related to smoothed densities used the paper. consider Y N(X, 1) where X ew ((() Rd, thus with the density
p(y) 7 f(x) 1
zo 11l do x f (-f(a)e g zoz w- da.
e
(2na2)d/22
We have the following expressions for log p(y):
p(g) log
e --(() (0021)=221 10=1==vll dx est,
logp(y) log / (-f(ottlee 2 |x|l do + cst,
leading to
v(log))()) e f(0) 2 1|=- 2 dr
(E|X y),
(B.1)
e (e-((a)e केत v- dx
j (e- V (xr u)e 25z 1x11 duc
v(logr)(y)
e Mo+tve z ||0|| dx
e [(x) Vf मेン 11e EE (X))(()))))
(B.2)
J (e-S(()e के u1l dx
This in turn leads to three expressions for the Hessian:
1
v (logr)(u)
[x(x y) 1w] El ylyl
1
E xx T Tly] - EEE - - cov
1
(B.3)
0²
(logg)(()) E[v²f(X]
E E f(X) cov (xx
(B.4)
(log))(y) 1
[[( - [[ -V
σ²
σ² [B(xVf(x)Tlv)
Eb E σ²
cov(X, Vf(X))y)
B. GAUSSIAN EXAMPIEE
If X N(p, C.), then Y N(H, C. + d2)), and then
[[X[Y] μ + = p),
= -C
Therefore, X| is Gaussian with mean μ and the covariance matrix C(C -C We have -V (v) (C + -2
B.2 DISTRIBUTION OF Y|
If Y's is sampled from py, let Po be the distribution of [X|]||| Y + (Y) We can bound the 2- Wasserstein distance between p and P.. by sampling X from px, and taking Y + 7 P(Y). where Y N(X,o21), This leads to particular "coupling" (Villani, 2021)
between p and P.. It follows:
We(p. p.)² & Ex,x - X
Given that the conditional expectation E[X|Y] is the optimal estimator for the square loss, the last expression is less than the trival estimator Y. that is,
W2(p.po> 22)2 & Ex,rlly -xpl/--XX---- =11=022
(B.5)
13

---

Published as a conference paper at ICLR 2024
B.3 PROPOSITION 2
The calculations in Sec. B.2 that resulted in (B.5) leads to the following:
Proof for Proposition 2. Using Proposition 1, we have
W2(p, Po.m) W2(p.pm-u1
which is then combined with (B.5):
o² =dd d. W2(p.Pamm) &
C PROOFS FOR SEC. 3
C.1 PROPOSITION 3
Proof for Proposition 3. We start with an expression for log p(y1mm):
m12
log p(g(())))) = log
II
N(x; yu 221) v
N(x; 0,C, C) da
Jzd
d
(yks =xa)²
xa
log
exp
20
2T
dx; + est
(C.1)
B
k=
(I)
a,)²
log
exp
Ys'
dr, cst.
2,3
The expressions for as, 8i, and % are given next by completing the square via matching second, first and zeroth derivative (in that order) of the left and right hand sides below
27
(x; a2,)
20
2T
233
%
evaluated x; 0. The following three equations follow:
1
m
3
σ²
2
m
m
0g
1/i
7 0j
ytir
83
m
+ 2²T
2
2
VEi
Vii
Σ Ycs
2,3
7
20
71
20
20 (m + -²²
1=1
[=1
1=1
It is convenient to define:
Au t O=
(C.2)
Using above expressions, (C.1) simplifies to:
1i 1
d
log p(11:m) =
+ est =
Ams
yes
+ est.
(C.3)
202
20
1=1
1=1
1=1
The energy function can be written more compactly by introducing the matrix Fa.mi
1
logp(yim)) = - 0² Fo = ((1 Ami)dw - Amis (1 Sec o..
2 77 Famwiim
14

---

Published as a conference paper at ICLR 2024
In words, the md x md dimensional matrix Fa.m is block diagonal with d blocks of size m. The blocks themselves capture the interactions between different measurements indexed by t and t': The x blocks of the matrix F Fa,ms indexed by i [d] thus have the form:
2 F(0) a,m
= (1 Am + Ami 1m 1m),
where 1 = (1, (1....... ..... 1) is m-dimensional. straightforward to find the m eigenvalues of the m x matrix F(8) for i E [d]:
degenerate eigenvalues equal to corresponding the eigenvectors
{(1, -0.-1,11 1,0, 0........ (L, -1..0..... (L,0, -.0.,..__0___................... }
eigenvalue equal to 21 mmmi corresponding the eigenvector (I, 1)T.
Since mAmi 0 have:
Amax (Fa =
which is (m 1 d degenerate. The remaining d eigenvalues are fo- mm mm) NE(d). the smallest of which is given by
Amin (Fom 0- -1-1111
711 + o Traax
+ mo T
Thus have:
Ka.m Amax (Fo.m Amin (minn (Fo,m) = 1 mo 22 max
C.2 PROPOSITION 4
Proof for Proposition 4. Using (C.3) and (C.2) have:
-20² log (ogp())1)) = Σ
18u
A1s
ykuk
A=1
k=
Since log p(y1111a-1) = logp(v1:)) log p(8)::11). we have:
1-1
-2 log logp(vllvi=-1) = 11ya
A Yei + Σ yki
+ est
Σ (1 And (uii - Ari
Σ
1-1
yki)) + est,
Aci
k=1
Therefore the conditional density p(yeel11-11) 1) is the Gaussian N(Hee-11 F FTTL1) with a shifted mean
Ass
Pij-1 = ( Ass
Σ yka
sEld
and with an anisotropic, diagonal covariance/precision matrix whose spectrum is given by:
= =1 X² -1- (t ?
Thus:
Amaae (F 1 (1 (1+
Kijc-1 Xuia (F(c-1) 1 1- (t 2
Lastly, to prove monotonicity result (3.3) we do analytic continuation of Kejt-1 to continuous values by defining 2t( Kije-1 and taking its derivative, below R 2 - mmax thus
15

---

Published as a conference paper at ICLR 2024
R >
(t R) (1 - (1-(+R)))((ttr))
(1())
1 r1+ (1 (1-
(1 + (1+r(()--]))))
-
- 1)(++ R)
(t ()++tr11111) (++r- ()+ Ry²
+ (R (1+r-1)((1 R)²
(r- (2 R -1-)
O.
(0+r-1)2(t R)²
D PROOFS FOR SEC. 4
D.I LEMMA
Proof for Lemma 1. We would like to find upper bound for (log By using 2T (logg))(v) - + o (ov
derived in Appendix B, we need to upper bound cov (X ly. It suffices to study [B|| xo 1oll
since
cov(Xlv) cov(X zolv) & E||| xo 1y]
We have by a convexity argument
-((0)= 2tb vll | xo 22
B|||X - xoll211 [b]
(((()--- मे -1 dr
(D.1)
(w)) to ( (0- (r~zo) ||
xi 22
.e --=((()++ T (y-2a) da
Next, we find an upper bound for || zoll2 itself. Using the assumption in the lemma we have:
[v) - (3r b) [
2|3000 vll >/ /a|
=ool
A a
PI1200 -
yll,
leading to /||l zoll < V + (xo - v)1l + o 2|110
yl and thus
|T0 y11? σ² Finally, using (D.1), we only need to find an upper-bound for [f(xx) + o 0 y) under the distribution p(x) e f(x), where f(z) f(x) σ - 2(x Øo) T (v zo). This is achieved with:
n² - zo! 2
V fx)
+
To -
12 31 2 +
f [
-0 / (())) dx = ((x) tr ( & Ld,
where second equality is obtained using integration by parts akin to score matching (Hyvärinen, 2005), and for the last inequality we used our assumption V (²f(z) & LI. Putting all together we arrive at:
3Ld
3A2
3|.00 - v11 I
(log))) & (-1+
1202 +
1202
1206 0
'We consider the exponential family
p(a|v) = eexp-f((>)
1 oe - xo) T (y m0)
1 1| xo
V 11 - zoll - a(v)).
2
2
Since a(v) is convex (Wainwright & Jordan, 2008), we have a'(1) > a'()), which is exactly the desired statement.
16

---

Published as a conference paper at ICLR 2024
D.2 LEMMA 2
Since log ((plliaaa-)) = logg()))) - logp(ynee-),, we have
V log (Ss((ell1e-1) V logp(via)
Start with the score function =):
3 ye ())(())))))))))))))())))) IT plui r)dr =0 (E( yr).
V de log p(y1())
( II p(y/ (r)dr
T + Ar BR, where A, and Br are given by: Next, we derive the Hessian V logp(vi:)) -0-
A. 1x - 31)( - ye) p(x) II p(uilz)dd o - E[( -))) we) [Ma).
fp(x) IT p(ys (x)
B, o 0- [ [(0-v)p(e ()) II p( p() (1(x) Ie p(u dx
[ (x--)p(<)]] v))(() II p(ullz)dr [p(x) II p(u))))
=o- (IE[X(B1ad] ye) (E[X y:)
By simplifying At Bi, the y cross-terms cancel out and the posterior covariance matrix emerges:
Vg logp(y::)) - - I
(D.2)
The lemma is proven since the mean of the posterior covariance Co (Xlv).) can only go down upon accumulation of measurements (conditioning more variables).
D.3 THEOREM
Proof for Theorem 1. We start with the definition of (a, m)-density:
L. p(z) (Ls exp
p(y1m) = | p(=)p(yim lz)dz x
||
o yell?
dx
dz,
2-
20
which we express by integrating out æ. We have
p(y)ml/>) x / exp
1 272 |1 2||2
2x [|s
F1.ml
1
1171.m
11
+
1
llvall
dx
m
exp
2(mr2 + 02) || 7 Iimml 112 +
ll
20
20
We then express pm log p(y)m)) ) in terms of E|Z
[mmm]:
you log logr())mm) @m(m²2 + 02)
mr?
Vim
ym
mn2 + o²
2 [1:m].
Next, we take another derivative:
1
2 Fm
logp(yimm)
@ + 0²)
I
mn² o2
V B V1
Next, we compute V ym E|Z [[
m
ElZZ
[1:ml
E|Z 11.m
F..
E Z 11:mm1m/
Vi
mr² + o²
mr² + o²
mr²
T
1 + o2
[1mml +
m²72 ṭ2
2 ltm
σ2 9m
m 22 cov( Z1312:m)
17

---

Published as a conference paper at ICLR 2024
Putting all together, arrive at:
V2 Um
log p(y)mm)
o² M +
7²
-1)
1+
(mT² + 02)² Cov (1)mm.
(D.4)
Finally, since
|| < R2 almost surely, we have cov cov(Zly):m) x R therefore
V log pyy(lm/0mmm111 = V Vm log logp(()mm) क ((m)
where
1
(
MTT - - -1) +
72
R
((m)
(mT² + 02)2
D.4 DERIVATION FOR (4.6)
In this example, we have X therefore
Z + No. No N(0,721). and Y X + N1, Ni N(0,021).
Y Z + N, N ~ N(O, (o² + 7)),,
(D.5)
where Z ~ p(z),
p(z) = 1 50(z p)) + 8(2 p))
is the Dirac delta function in d-dimensions. Alternatively, we have ply) = J p(=)((y)==))))) where
p(yl=) N(y; (o² + 1))1)..
Using (B.3), adapted for (D.5), we have:
( = - v (logg) o² + 72 (1- (1- 72 cav (2lw) Next, we derive an expression for cov cov(Zly):
(D.6)
02
cov cov(Z y) E|ZZ' [g] [[
We have:
E T lv] = 14
A e
E
+ e
A |1y 11
2(02 + T2)
|1 11l2
B
2(02 +2)
It follows:
cov ud T. (1-
e A - c B 2)
2141
e e-A + c-B c
1 1 cosh (B A)
(D.7)
2p y)
-1
= =2up T.
1 cosh
2 + 72
By combining (D.6) and (D.7) we arrive at (4.6).
"oot required for the proof, but we can also relate cov ( X (lymm and cov (Z|ytm)) directly since as we know (see the proof of Lemma 2):
V fm log 1oyi())).))
a² + ao cov(X m)m)).
(D.3)
By combining (D.4) and (D.3), it follows:
o27²
cov X (l)1mm)
1
mT² + 0z cov(Zlui.m)
m7² + g2
18

---

Published as a conference paper at ICLR 2024
E DETAILED ALGORITHM
Algorithm 2: MCMC, in Algorithm via Underdamped Langevin MCMC by Sachs et al. (2017)
1: Input y(i-1), 8 1:t-1
2: Parameters current MCMC iteration i, current measurement index t, step size 8, friction 7
steps taken ne, estimated smoothed score function (vv: Lipschitz parameter noise level
3: Output r,(
4: Initialize y: (1,0) Unif([O, []4) + N(0,
5: Initialize t
6: for k (0......
do
7: yliiht1) (8.k) + V
8: FI = =Y F + (Y, (r((ib11)_ - -
G t (Y 1=/t-1122 + 02V y((h+l) according to (4.9)
9:
10:
V V
$ 2I G
B N(O.T) I)
12:
V exeppp-90 V +
5 G VL
2I
y((k+t1) y, (1,t+)))
1v v
13:
14: end for
return y,(1 ((i
15:
F ESTIMATORS FOR V log)(y)
The purpose of this section is derive estimators for the smoothed score function
و σ) V(logp) V)
which can be used to run Algorithm 2, in turn running Algorithm 1. We first derive (4.8) by rewriting (B.1) as follows:
1 f(x - y) exp (x-)))x((-((tvv)))e g(v;o)= = 0² [exp(-f(r:y,
σ2 (E - y). X -
where
((5,y,0) = f(x) | 8/1²
20
F.1 THE PLUG-IN ESTIMATOR Iplugin
We can arrive at simple plug-in estimator for g(y; by rewriting (B.1) as follows:
EN [Ne (NNV((+
V G² Ex-N s(N+w)
We then simply estimate the numerator and denominator above using i..d. Gaussian draws:
1 Σ exp exp (-f(ae; + 8))
6, itd N(O, T).
Iplugin (y: 0) o Σ exp (-f(ef + v))
(F.1)
(Above, we used one set of i.idd. draws for estimating the numerator and denominator used in our experiments. Of course, one should take independent draws if computation budget is not an issue.)
Gradient-based plug-in estimator. One can also obtain a different type of plug-in estimator for Vlog p(y) using (B.2) that directly takes the gradient information V f((z) into consideration. This estimator should have better properties, but we did not experiment with it in this paper. Intuitively, such plug-in estimators (gradient-aware not) will suffer from the curse of dimensionality.
19

---

Published as a conference paper at ICLR 2024
Connections to importance sampling. One can also arrive at (F.1) by using importance sampling to estimate [[X] in (4.8), which is insightful. In particular, in importance sampling instead of sampling from the probability measure v associated with exp(-f(:: y, σ)) (which hard to sample from) we sample from easier probability measure μ4, which in Gaussian measure N(y; I)). With this setup, the estimator (F.1) 1) follows through. It is known that the number of samples required to have accurate estimator based on importance sampling is of order exp( (L)))))) where (L(w/|k) is the Kullback-Leible divergence of μ from v (Chatterjee & Diaconis, 2018). This further highlights the limitations of the estimator (F.1).
took it to be the
F.2 THE LANGEVIN ESTIMATOR Ilangevin
better way to estimate [[X] in (4.8) is to use MCMC - in particular, the gradient- saadeent-based Langevin MCMC. For any y, this is done by running Langevin MCMC using the score function
y r²
0( (x)
As example, and to be concrete, for any sample X by discretizing the overdamped Langevin diffusion:
X.- x² ds + 2dB,
Xx, = (-v(()v()))
Given n such independent samples, (X:) drawn by running the Langevin dynamics, we then use (4.8) to arrive at Ilangevin:
Iaangevin (3) = σ² n2
(-
1
-
X; - y).
(F22
The procedure above for drawing a sample X from (x:go)))))) using overdamped Langevin MCMC is given in Algorithm 3.
Algorithm 3: Draw X (((eig)) with Langevin MCMC for the estimator (F.2)
1: Parameter noise level a
2: Input current y for which we need to estimate ( a))
3: Hyperparameters step size 6, number of iterations
4: Output X
5: To
6: for k 1,...., K] do
7:
ÿ = -V((x=11 0- (-x-11 - y)
8: E 1)
9: Fk zk-1 + ōg + V28
10: end for
return 6 FK
We should highlight that in our experiments, instead of using the overdamped Langevin MCMC above, we used the more sophisticated underdamped Langevin MCMC algorithm by Sachs al. (2017).
Connections to "Entropy SGD". The problem estimating smoothed score functions have also been of interest in the neural network optimization literature under the terminology of entropy SGD (Chaudhari et al., 2017). There, the smoothed score function (although the kernel smoothing lexicon not used there; 1/0² is denoted by y and is referred to scope) utilized for optimization, not for sampling, to arrive at flatter minima of the loss landscape with better generalization properties. Due to different motivations, Algorithm 3 differs from the one used in entropy SGD, e.g., we do not have exponential averaging here.
20

---

Published as a conference paper at ICLR 2024
G EXPERIMENTAL DETAIL
The hyperparameters were tuned on a log-spaced grid. We searched the step size over {0.03, 0.1, 1.0), the effective friction ÿ6 over (0.0625, 0.125, 0.25, 0.5, 1.0), per-t MCMC iterations ny over {1, 4, 16), and the Lipsehitz parameter over f1/02, 1.0). We found that the hyperparameter combinations (Ō 0.03, 76 0.0625) and (6 1.0, 76 0.5) worked well for most configurations of test density type, o, d, and MCMC algorithm. For JMS, m ran over {200, 400, 600, 800, 1000}.
We experimented with two initialization schemes in SMS. At each t, the walkers sampling from P(u (( were initialized at warm: [[X E, where N(0, o² ) or (ii) cold: at samples from Unif 111) + N(0, o²). For convergence experiments, we report results from the warm start, as it was more robust the choice of hyperparameters. The tunneling results were obtained with cold start. The hyperparameters were tuned for each algorithm while the total number of iterations was kept fixed at large value. We define each iteration as an MCMC update step. For SMS, the total number of iterations is Σx=1 nz. the number of MCMC iterations for each measurement t, ne, summed up over the m measurements. We had n, E {1, 4, 16) and m = 5, 000. For the remaining three sampling schemes (JMS, m 1, and a 0), the total number of iterations is simply the number of MCMC iterations, fixed to 20, 000, but we found these algorithms to converge much earlier, around 5, 000.
H MCMC ALGORITTM
Our algorithm is agnostic to the choice of MCMC sampling algorithm used in the Markovian phases. In this section, we run SMS sampling with four different Langevin MCMC algorithms. The results presented earlier in Sec. 5 uses an ULD algorithm with an Buler discretization scheme that extends the BAOAB integration using multiple time steps for the O-part ("Sachs et al.") (Sachs et al., 2017). Next, consider two algorithms that operate the integral representations of ULD. Recall that continuous- tie ULD is represented by the following stochastic differential equation (SDE):
dv: = -o;dt a (x;)dt + (ettt++2222222 B..
dr. v.dl,
where Ht, 12 E Rd and B, is the standard Brownian motion in Rd. The solution (Iz, 0r) to the continuous-time iL
11 Uhe - (L/ exeexp(( (-(p (s(r))ds)
((x)))))))(( /2 ( exp exp (->(1 8)) dB,
84 = 2o + S' unds.
(H.I)
Similarly, the discrete ULD is defined by the SDE
dō; -yydtt u dt + (v dy dē, = e,
which yields the solution
ú. -ūne - ((-exp(-((-s)) v((oodds + //77/ 1' exp (xp(-1(t-0)) dB)
2, To + f u,ds.
(H.2)
Shen & Lee (2019) seeks lower discretization error by using 2-step fixed point iteration method, or the randomized midpoint method. The integral in (H.1) is evaluated along uniform random points between and On the other hand, Cheng et al. (2018) computes the moments of the joint Gaussian over (ă,, 5x) in the updates of (H.2). In our comparison we additionally include MALA, an Euler discretization of the overdamped Langevin dynamics represented by the SDE
dx. - + (V 22d)dB..
(H.3)
21

---

Published as a conference paper at ICLR 2024
MILA Stdhoit al.
Chera et al
4.0
Dosett sumples of
2.0-
6.0
200
000
500
1000
Figure 5: W2 vs. m for various MCMC, algorithms used the inner loop of Algorithm
accompanied by Metropolis adjustment to correct for the discretization errors (Roberts & Tweedie, 1996).
The algorithms are compared in Fig. for the Gaussian mixture test density introduced in Sec. 5.2 with d 8. The three (unadjusted) ULD algorithms converge faster than does MALA lower W2. In ULD algorithms, Brownian motion affects the positions through the velocities U,, rather than directly as MALA, resulting in a smoother evolution of z; that lends itself better to discretization. The first two dimensions of the final samples are displayed in Fig. 6. MALA samples fail to separate into the two modes, whereas Cheng et al, Shen and Lee, and Sachs et al have better sample quality, with Sachs et al performing the best and almost approaching the sample variance of W2 when samples are directly drawn from px (the gray band).
(a) MALA (Roberts & Tweedie, 1996)
(b) Cheng et al. (2018)
>
(c) Shen & Lee (2019)
(d) Sachs et al. (2017)
Figure Final X samples for various MCMC, algorithms used in the inner loop of Algorithm
22

---

Published as a conference paper at ICLR 2024
I SCORE ESTIMATION
Langevin MCMC requires the smoothed score function g 0) = V pog) Results presented earlier in Sec. 5 assumed access the analytic score function. Here, we use the plug-in estimator Oplugin presented in (F.1) with varying numbers of MC samples n as well as the Langevin estimator Aaaggvvin presented in (F.2).
To prevent numerical underflow, we implemented (F.1) as follows:
logsumexp (-f(y + os;));
(1.1)
B_ logsumexp 1 (loge; f(y + as;));
(1.2)
B_ logsumexp f- (log(-E fly + oεx)),
(133
1 (.(e (e Bt -A (
Oplugin (w)
(1.4)
where
n1
logsumexpi= (a.) = Qmax + log Σ exp(a: - amax):
4max =
max, a,, and N).,). (I.2) and (1.3,, L,..., n4 and k = I,..., n. denote the indices for which is positive and negative, respectively, with n4 + n_ n. Note that the same Gaussian samples E were used to evaluate the numerator and the denominator.
For Gaussian mixture density introduced in Sec. 5.2 with d 2, both the analytic and the Oplugin score functions converge to a low W2, as shown in Fig. 7 (a, b). The sample quality on par with the analytic score function with small 500 and there is little benefit increasing the past 500. The Dangevin score function struggles in this low-dimensional example, but catches up to the highest-n Oplugin score function at d 8, as Fig. 7 (c, d) shows. All estimated score functions significantly underperform the analytic score. This is preliminary result: for obtaining Glangevin we did not extensively tune the hyperparameters. In addition, for Iplugin variance reduction techniques. such as importance weighting, may help get more mileage from finite n.
J MIXTURE OF CORRELATED GAUSSIANS
In this section, we study a correlated test density, namely a mixture of two 2-dimensional Gaussians with full covariances:
p(x) a N(;;; En) (1 0z) N(x; -4, En).
(J.1)
We choose = 3 1a, Eo = R diag( 4. 4) RT, Σ RD diag(1, 9) R, and the rotation matrix
cos 0
sin cos 0
sin 0
with = x /660 Fig. 8 compares the performance of the same Langevin MCMC algorithm (Sachs et al., 2017) used within the SMS scheme (a) with that used without Gaussian smoothing (b). In both cases here, we initialized the samplers at the origin (5000 particles total). The total number of iterations are the same between two algorithms (see Appendix G for details).
23

---

Published as a conference paper at ICLR 2024
10.0
0.35
Estimater Andlytie
Direct aample of
0.30
250
== 4000
R - 200
Analytic
0.25
R1 = 1000
MCMC
W2
PX
0.20
0.15
0.10
0.1
0.05
0.00
100
200
300
400
2000
4000
Number of MC samples
(a) 22 vs. m (d = 2)
(b) Final W2 vs. for Oplugin (d = 2)
3.0
10.0
2.5
2.0
1.0
Analytic
W2
1000
MCMC
Estimater
M 2000
X ~
Analytic
Direct samples px
0.1
100
200
300
400
250
300
1000
2000
1000
Number MC samples
(c) W2 vs. (d=
(d) Final W2 vs. for Iplugin (d
Comparison of the analytic g(y; o) V log p(y), the plug-in estimator 9plugin (F.1) with varying numbers of MC samples n, and the Langevin estimator Alangevin (F2). The test density was the mixture of Gaussians introduced in Sec. 5.2 with d (a, b) and d = 8 (c, d).
Figure 7:
(a) SMS, = 16, m 1000
(b) Langevin MCMC
Figure 8: Mixture of correlated Gaussians.
24


Reviews:
[{"review_id": "GUk7xM4jFm", "reviewer": "Reviewer_ZJfB", "summary": "This paper presents a framework for sampling from an unnormalized density $p\\propto{\\rm e}^{-f}$ on $\\mathbb{R}^d$ by using a Gaussian smoothing scheme. It constructs $X\\sim p$, $Y_i|X\\sim_{\\rm i.i.d.}\\mathcal{N}(X,\\sigma^2I)$ ($i\\in[m]$) for some suitable $\\sigma$ and $m$. By sequentially (one at a time) sampling the joint distribution of $Y_{1:m}$, the conditional distributions will become \"more log-concave\", and finally it outputs $\\mathbb{E}[X|Y_{1:m}]$ as an approximate sample from $p$, which can be expressed in a function of $\\bar{Y}_{1:m}$. The paper demonstrates the advantages of the sequential sampling scheme through some theoretical and empirical examples.", "questions": "1. The algorithm outputs a *biased* sample of $p_X$, but the bias $W_2^2(p_X,\\hat{p}_{m^{-1/2}\\sigma})$ can be made arbitrarily small by choosing $\\sigma$ and $m$ such that $\\frac{\\sigma^2d}{m}$ is small enough. How can we choose these parameters in a principled way? In the Gaussian example (proposition 4), we can reduce $\\kappa_1$ by increasing $\\sigma$, and in theorem 1, we need a large $\\sigma$ to make $p(y_1)$ strongly log-concave. However, this also increases the bias, so we need a larger $m$ to compensate that. How can we balance this trade-off between complexity and accuracy?\n\n2. In the experiment, how do you choose the projection direction $\\theta$ for estimating the Wasserstein-2 distance? A more fair comparison would be to use $\\theta$ following the uniform distribution on the Euclidean unit ball (i.e., the *sliced* Wasserstein-2 distance). Efficient approximation algorithms are available, for example, in (Nadjahi et al., 2021).\n\nMinor corrections:\n\n1. In A.1, $-2\\sigma^2\\log p(y_{1:m}|x)=\\sum_{t=1}^{m}||y_t-x||^2+{\\rm cst}$. \n \n2. In remark 3, for (4.7) to hold, ${\\rm cov}(Z|y_{1:m})\\preceq{\\rm cov}(Z|y_{1:m-1})$ almost surely; ... is very large, so that ${\\rm cov}(Z|y_{1:m-1})\\approx0$.\n\n3. The exponential family in the third footnote should be $p(x|\\nu)=\\exp(-f(x)+\\frac{1}{\\sigma^2}(x-x_0)^\\top(y-x_0)-\\frac{\\nu}{2\\sigma^2}(||x-x_0||^2+||y-x_0||^2)-a(\\nu))$. \n\n4. In the third paragraph of introduction and the second paragraph of section 2, $Y_t=X+N_t$, $t\\in[m]$. $N_t$ are i.i.d. $\\mathcal{N}(0,\\sigma^2I)$ random variables that are also independent of $X$. In theorem 1, $X=Z+N_0$, $Z$ and $N_0$ are independent.\n\nReferences:\n\n(Nadjahi et al., 2021) Kimia Nadjahi, Alain Durmus, Pierre Jacob, Roland Badeau, Umut Şimşekli. Fast Approximation of the Sliced-Wasserstein Distance Using Concentration of Random Projections. In Advances in Neural Information Processing Systems, 2021.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "Overall, the paper introduces a novel theoretical framework that differs from existing MCMC methods or sequential inference models. The paper is well-written, with clear notations, convincing examples, and mathematically rigorous proofs. The structure is well-organized: the authors first introduce an example of an anisotropic Gaussian target distribution, which leads to the important observation that OAT is easier than AAO in terms of the condition number. Then, they present the theory for general distributions (once log-concave, always log-concave), which helps the reader understand the motivation.", "weaknesses": "Score estimation is essential for the algorithm's efficiency. However, most of the experiments in the paper use the analytic score function, which is not realistic in practice. Moreover, the score has to be estimated at each step of ${\\rm MCMC}_\\sigma$, which will significantly increase the algorithm's complexity.\n\nThe paper only compares the performance of accurate and approximate scores in appendix H, and the result is not very satisfactory. The estimated score performs well when $d=2$, but it deteriorates when $d=8$. These results suggest that the score estimate's precision could be much worse in high dimensions, which I believe is the paper's main weakness.\n\nIn 4.2.1, the proposed method for score estimation is based on importance sampling (which transforms the expectation w.r.t. the probability density $\\nu\\propto\\exp(-f-\\frac{1}{2\\sigma^2}||\\cdot-y||^2)$ to $\\mu=\\mathcal{N}(y,\\sigma^2I)$). However, (Chatterjee and Diaconis, 2018) has shown that the sample size needed for accurate estimation by importance sampling is usually $\\exp({\\rm KL}(\\nu||\\mu))$ for general target functions. For this reason, (Huang et al, 2023, section 3) use Langevin Monte Carlo to sample from $\\nu$ and combine it with importance sampling in a similar task of score estimation. The experiments in your paper indicate that the score estimator also suffers from high variance. The authors should experiment with different score estimators, instead of saying that \"studying the covariance of the plug-in estimator and devising better estimators is beyond the scope of this paper\".\n\nReferences:\n\n(Chatterjee and Diaconis, 2018) Sourav Chatterjee, Persi Diaconis. The sample size required in importance sampling. Ann. Appl. Probab. 28(2): 1099-1135 (April 2018). DOI: 10.1214/17-AAP1326.\n\n(Huang et al., 2023) Xunpeng Huang, Hanze Dong, Yifan Hao, Yian Ma, Tong Zhang. Reverse Diffusion Monte Carlo. ArXiv preprint arXiv:2307.02037.", "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": "qn6i1IhhJb", "reviewer": "Reviewer_9mSy", "summary": "This paper proposes a framework of sampling from unnormalized densities by transforming this problem to more amenable log-concave sampling.", "questions": "1. How to handle the case when the desired probability measure does not have a density?\n2. It would have been insightful to discuss adversarial cases where the proposed sampling technique fails.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "**Originality**: The paper is an extension of the paper by Saremi & Srivastava. The proposed framework is very different from what I could find in the literature.\n\n**Quality**: Once familiar with the literature (see weaknesses), I found the framework really elegant and simple. I learnt a lot from reading this paper.\n\n**Clarity**: The paper is well-written (again see the caveat in the weaknesses section) for someone already familiar with this literature. The calculations are thorough without being burdensome. The explanations are crisp without being terse.\n\n**Significance**: The final algorithm proposed is simple and should be easily adopted by the community. More importantly, I found the paper useful in learning how to think about sampling, and I believe it will be useful for the rest of the community also.", "weaknesses": "The paper could have done a better job at exposition of not so well-known notions like walk-jump sampling. The current exposition is good enough for someone intimately familiar with the literature, but for others, like me, it takes some reading of cited papers to not see ideas as being \"pulled out of a hat\".", "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": "xBy0bWTMNF", "reviewer": "Reviewer_gqez", "summary": "The paper introduces an MCMC-type method that sequentially samples an unnormalized density with Gaussian smoothing. The critical theoretical contribution shows that with a proper choice of noise level and the sequential strategy, the conditional densities with compact support (or similar ones) become more and more log-concave \"healthily.\" Through the process, the algorithm mainly keeps the means, given the properties of the Gaussian kernels. Internally, the algorithm employs an MCMC strategy (like Langevin) and methods for score function estimation.", "questions": "Please take a look at the \"Weaknesses,\" where the comments cover both theory presentation and experiments, as well as some additional questions (under \"minor points\"). \n\nOne additional question: \n- Is the actual code anonymously shared somewhere? It would be good to have it for result reproducibility and smooth adoption of the method.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The submission is a well-written research paper. The authors explained the concepts clearly and presented the critical theoretical contribution elegantly. The authors have done an excellent job of demonstrating the effectiveness of the proposed method and its potential applications. \n\nThe problem itself is also fundamental and worth studying. The authors start by finding that \"all (measurements) at once\" (AAO) is suboptimal and then present a \"one (measurements) at once\" (OAT) algorithm that has solid theoretical guarantees. In particular, Gaussian smoothing helps transform a (nearly) compact density towards log-concavity so that it's feasible to sample from it, and measurements accumulation (in a non-Markovian manner) makes the density more and more log-concave and well-conditioned. Combining both effects makes the proposed OAT strategy superior to AAO ones.\n\nMoreover, the paper demonstrates the value of combining theory and experiments in a research project. The theoretical analysis provides a solid foundation, while the experiments validate the proposed method and help readers gain more insights and understand its accuracy. Combining theory and experiments leads to a more complete understanding of the OAT strategies and opens up new avenues for further research.", "weaknesses": "One weakness is in the result presentation: The paper does not provide a direct performance guarantee on Algorithm 1. I understand the inner-loop method is the object of study here. Still, it should be feasible to make some assumptions around its properties and derive a theorem that captures the algorithm's performance as $m$ and $n_t$ increase. The error metric could ideally be a Wasserstain-type distance for consistency with the experiment section. \n\nThe experiment evaluation partially supports and reflects the theoretical claims, but the settings could be more complex for a thorough evaluation. For example, it would be interesting to understand what magnitude the # sampling iterations must be for some rather complex distributions, e.g., $k$-mixture of Gaussians where $k$ increases. The densities used in Section 5 are well-designed to test specific hypotheses/claims. Yet, experiments should probably cover more than that, such as demonstrating the properties of the main algorithm(s), which can be critical for method adoption.\n\nMinor points: \n- May I know the purpose of \"4.1 Example\" as it seems to be just deriving various quantities for a concrete/standard Gaussian mixture? I don't see how it helps the readers better understand Theorem 1.\n- In Section 5.2, is there an intuitive explanation for why \"the optimal $\\sigma$ here is in fact larger than the noise level needed to make $p(y_1)$ log-concave\"?\n- Would it make sense to add more description for Figure 4, panels (c) and (d)? Currently, it is not entirely clear what they are offering.", "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": "FBd1xPbojs", "reviewer": "Reviewer_3oQn", "summary": "This paper introduces a general framework for sampling  unnormalized probability densities. The main idea is based on smoothing the distribution\nby adding Gaussian noise to the original sample. The process then repeats $m$ times,\neach time independently adding gaussian noise to the original sample.\nThe final step involves computing the Bayes optimal predictor of the original sample based on the i.i.d. samples from the smoothed density.\n\nThe authors first demonstrate via an example that sampling all the smoothed observations jointly at once might not be optimal, since the \nsampling problem might become ill-conditioned as $m$ increases, which prohibits the use of standard MCMC methods. Subsequently, they demonstrate\nthat if instead sequential sampling is used, where each sample is taken conditional on the previous samples, then the condition number\nof the sampling problem can only improve as $m$ increases, enabling the use of MCMC methods. They prove that this process succeeds to sample\nin a general setting where the initial distribution is supported on a bounded subset of $\\mathbb{R}^d$. Additionally, they evaluate their\napproach experimentally, by sampling from elliptical gaussians (which are ill-conditioned) and mixtures of gaussians. \nThe performance is then compared to that of standard Langevin MCMC algorithms without smoothing, showing significant benefits both\nin terms of distributional convergence, as well as finding less dominant modes of the distribution.", "questions": "-Related to a previous point about estimating the score function, how exactly do the authors arrive at the expression (4.8) for estimating the\nscore function? In appendix B, a number of properties is proven about the score function, among which that it is related to $\\mathbb{E}[X|y_{1:m}]$.\nHowever, I don't see how this is captured in the empirical definition of (4.8). What is the meaning of taking a weighted average of\nthe $\\epsilon_i$s?\n\n-As the authors discuss in Remark 2, the convergence of the condition number to $1$ with rate $1/m$ holds as long as the original distribution involves \nsome (even small) additive gaussian noise. Why is it enough that even a tiny bit of Gaussian noise be added to the distribution for this decay to happen?\nDo the authors have some intuition about that statement?\n\n-Related to the previous question, if we care about the expected hessian of the conditional sampling, I'm wondering whether it's possible to have a version of Theorem 1 showing that the expected condition number of the conditional sampling problem\nconverges to $1$ with rate $1/m$ without assuming that the original distribution contains some amount of Gaussian noise. The reason is\nthe following: the proof for all the claims about the \"increase\" in log-concavity basically hinge on equation (D.2), which asserts that\n$\\nabla^2 \\log p(y_{1:m}) = -\\sigma^{-2} I + \\sigma^{-4} Cov(X|y_{1:m})$. The argument then proceeds by the observation that conditioning on more and\nmore observations can only reduce the variance of $X$ in expectation, so the \"positive\" part on the right hand side can only become smaller as\n$m$ increases. The proof of Theorem 1 expresses $Cov(X|y_{1:m})$ as $O(1/ m) + O(1/m^2) Cov(Z|y_{1:m})$ \nand bounds the latter covariance by a constant, since $Z$ is supported on a bounded set (as the footnote in page 18 suggests). \nHowever, it seems conceivable to me that if we do not impose the assumption s $X = Z + N_0$, we could instead\nshow directly that $\\mathbb{E}[Cov(X|y_{1:m})]$ decreases as $O(1/m)$ if $X$ is bounded. The reason I believe this is true is that, by Proposition 2, we know that the posterior\nof $X$ conditioned on $y_{1:m}$ converges in distribution to $X$ at a rate of $1/m$. Thus, the conditional variance of $X$ given $y_{1:m}$ \nshould decrease as $m$ becomes larger. Another intuitive way to see this is that we could for example simply take the average \n$\\overline{y_m}$, which is $O(\\sigma/m)$ close to $X$. Thus, we gain more and more information about $X$ as $m$ increases, which suggests\nthat $\\mathbb{E}[Cov(X|y_{1:m})] = O(1/m)$ always. Have the authors considered this more general statement?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "This is a very original work, introducing an elegant and intuitive general framework for sampling from unnormalized probability densities. The idea of smoothing \nand walk-jump sampling has already appeared in prior work (Saremi & Hyvarinen, 2019), but this is the first time it has been applied to\nthe general problem of sampling from non-log concave densities. I like the intuition given by the authors about how smoothing helps by \"filling in\"\nprobability mass between modes of the original distribution, thus reducing the general problem to the log-concave setting.\nOne strength of this new approach is that the sequential sampling can simply be implemented by keeping a running average of all the previous\nsamples up until this point, so the memory footprint is small. \nMoreover, this paper is written with great clarity, as all theoretical results are explained in detail and they convincingly demonstrate the power of the new approach.\nThe experiments also complement the theory very nicely, by showing the benefit of running Langevin with smoothing. The tunneling phenomena that\nare observed from the simulations are also very interesting. Lastly, most of the proofs follow from standard techniques, but they are at an\nadequate level for ICLR.", "weaknesses": "One thing that is missing from this paper is the final step of computing the bayes estimator $\\mathbb{E}[X|y_{1:m}]$ to obtain samples from the true \ndistribution.\nFor that, it is needed to estimate the score function $\\nabla \\log p(y_{1:m})$ of the smoothed distribution. The authors acknowledge that this is \na separate issue that needs to be addressed. I understand this might be beyond the scope of the current work, which is more focused on \nshowing the benefits of smoothing. However, ultimately estimating the score function is a crucial step for this to be an end-to-end approach\nfor sampling and it is not clear what is the difficulty of this problem compared to the original sampling problem (is it strictly easier in some sense?).\nIn my opinion, this would make clear what the benefit of this approach over standard MCMC algorithms is.", "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)