Relative Entropy Gradient Sampler for Unnormalized Distributions
Authors: Xingdong Feng
Yuan Gao
Jian Huang
Yuling Jiao
Xu Liu
Venue: ICLR 2022 Submitted (2022)
Keywords: Density ratio estimation
gradient flow
neural networks
particles
velocity fields

Abstract: We propose a relative entropy gradient sampler (REGS) for sampling from unnormalized distributions. REGS is a particle method that seeks a sequence of simple nonlinear transforms iteratively pushing the initial samples from a reference distribution into the samples from an unnormalized target distribution. To determine the nonlinear transforms at each iteration, we consider the Wasserstein gradient flow of relative entropy. This gradient flow determines a path of probability distributions that interpolates the reference distribution and the target distribution. It is characterized by an ODE system with velocity fields depending on the  density ratios of the density of evolving particles and the unnormalized target density.  To sample with REGS, we need to estimate the density ratios and simulate the ODE system with particle evolution. We propose a novel nonparametric approach to estimating the logarithmic density ratio using neural networks. Extensive simulation studies on challenging multimodal 1D and 2D distributions and Bayesian logistic regression on real datasets demonstrate that the REGS outperforms the state-of-the-art sampling methods included in the comparison.

Under review as a conference paper at ICLR 2022
RELATIVE ENTROPY GRADIENT SAMPLER FOR UN- NORMALIZED DISTRIBUTIONS
Anonymous authors Paper under double-bling review
ABSTRACT
We propose a relative entropy gradient sampler (REGS) for sampling from unnor- malized distributions. REGS is a particle method that seeks sequence of simple nonlinear transforms iteratively pushing the initial samples from a reference dis- determine
tribution into the samples from an unnormalized target distribution. the nonlinear transforms at each iteration, we consider the Wasserstein gradient flow of relative entropy. This gradient flow determines path of probability dis- tributions that interpolates the reference distribution and the target distribution. It is characterized by an ODE system with velocity fields depending on the density ratios of the density of evolving particles and the unnormalized target density. To sample with REGS, we need to estimate the density ratios and simulate the ODE system with particle evolution. We propose novel nonparametric approach to esti- mating the logarithmic density ratio using neural networks. Extensive simulation studies on challenging multimodal ID and 2D mixture distributions and Bayesian logistic regression on real datasets demonstrate that the REGS outperforms the state-of- sate-offthe-aad sampling methods included in the comparison.
INTRODUCTION
Sampling from unnormalized distributions plays fundamental role in statistical inference and machine learning. This problem is frequently encountered in Bayesian statistics. Conducting Bayesian analysis requires evaluation of multi-dimensional integrals where analytical expressions for unnormalized posterior distributions are usually not available. Consequently, sampling is necessary for Monte Carlo approximation of these integrals. In this work, we propose general purpose sampling algorithm for unnormalized distributions.
Markov chain Monte Carlo (MCMC) methods (Andrieu et al.) 2003 Brooks et al. 2011) are widely used to sample from unnormalized distributions. Sampling with MCMC relies defining an appropriate transition kernel to construct a Markov chain whose equilibrium distribution is precisely the target distribution. Based on rejection sampling, the Metropolis Hastings algorithm (Metropolis) ctt al. 1953; Hastings 1970 Tierney| 1994] Dunson & Johndrow 2019) provides a flcxible framework needs to specify proposal density and an acceptance policy. However, without a careful design of these two aspects, the Metropolis-Hasting algorithm can be ineflicient due strong correlations, slow mixing, low acceptance rates, especially in the large-scale and high-dimensional settings, Moreover, proposals through discretizing some continuous processes like Langevin diffusion and Hamiltonian dynamics are introduced (Roberts & Tweedie 1996) Roberts & Stramer 2002; Duane et al., 1987, Neal 2011 Hoffman & Gelman 2014) and further enhanced by stochastic gradient estimation (Wellig & Teh
for general MCMC sampling. To implement a Metropolis-Hasting algorithm,
201 Chen et al. 2014).
Variational Bayesian inference Beal 2003), often simply referred as variational inference (VI) (Wainwright& & Jordan 2008 Blei et al. 2017), is another prominent approach to sampling from unnormalized distributions. VI approximates the unnormalized posterior distribution with a restricted parametric variational posterior distribution by minimizing the Kullback-Leibler (KL) divergence between them. Since the true posterior distribution is intractable, VI turns to maximize a surrogate variational objective called the evidence lower bound (ELBO). However, one is required trade off the parameterization flexibility of variational posteriors against the optimization complexity of ELBO in practice.

---

Under review as a conference paper at ICLR 2022
(a) 9 Gaussians
(b) 25 Gaussians
(c) 49 Gaussians
(d) 81 Gaussians
Figure 1: Scatter plots of generated samples and histograms of generated sample counts according to the nearest neighbor mode by REGS for mixtures of 9, 25, 49, and 81 Gaussians with equal weights. As the plots indicate, generated samples by REGS cover every component of the mixture distributions and are nearly equally allocated all components.
In the spirit of VI, particle-based variational inference (ParVI) (Liu & Wang 2016 Chen et al. 2018 Zhu al. 2020) iteratively optimizes for minimizing the KL divergence. ParVI seeks to move variational distribution towards the unnormalized target distribution, along steepest descent direction of the KI. divergence. In continuous view, these movements of variational distributions can be understood as a gradient flow in probability measure spaces (Liu et al. 2019a b). A key part ParVI is how to estimate the desired steepest descent direction (i.e., functional gradient) from the evolving random particles. An elegant approach is the Stein variational gradient descent (SVGD) Liu & Wang (2016). In SVGD, the functional gradient descent is embedded in reproducing kernel Hilbert space (RKHS), which is further recognized as a gradient flow under the Stein geometry (Liu 2017 Lu et al.) 2019) Duncan al. 2019). A drawback of SVGD is that tends to collapse part of the modes of the target, due to negative correlation between the data dimensionality and the repulsive force in the RKHS (Zhuo) al. 2018
particles mimic functional gradient descent
In this work, we propose a relative entropy gradient sampler (REGS) for sampling from unnormalized target distributions. To approximate target distribution, we consider the Wasserstein gradient flow of relative entropy (or KL divergence). named relative entropy gradient flow. The relative entropy gradient flow represents a path of probability distributions that follows the functional gradient descent direction of relative entropy. There exists an ODE system of random particles that uniquely determines the spatial and temporal dynamics the relative entropy gradient flow. Therefore, to sample with REGS, only need to simulate the ODE system with particle evolution. Evaluating the velocity fields of this ODE system can be transformed into estimating the logarithmic density ratio between the density of evolving particles and the unnormalized target density. Based this observation, we propose a novel logarithmic density ratio estimation method for unnormalized distributions. By collect of stable
alternating between particle evolution and velocity field estimation, particles which are approximately distributed as the target distribution. Our contributions can be summarized follows:
(1) Building upon the relative entropy gradient flow, propose the relative entropy gradient sampler (REGS) for unnormalized target distributions. REGS preserves high efficiency and strong stability with respect increasing singularity in mixtures of Gaussians, when the number of components increases (as shown in Figure []. the variance of each component decreases, and the distance between any two components increases.
(2) We propose to directly estimate velocity fields of the relative entropy gradient flow as gradients of logarithmic density ratios, that is computationally stable and efficient.
(3) We develop nonparametric approach to estimating the density ratio between an unnormalized density and an underlying density represented by samples, which of independent interest.
(4) We present experimental comparisons on varieties of multi-mode synthetic data and benchmark data and demonstrate that REGS is more accurate sampler than the popular samplers including ULA, MALA and SVGD.
Related work The proposed REGS is most related to sampling methods based on the relative entropy gradient flow, in particular, the recently proposed SVGD (Liu & Wang 2016 Liu 2017). which estimates the velocity fields of the relative entropy gradient flow in reproducing kernel Hilbert space. See also Korba et al. (2020); Salim et al. (2021 2020) for theoretical analysis of SVGD. In contrast, REGS approximates the velocity fields based on a novel logarithmic density ratio estimation approach with deep neural networks. The undesirable mode collapse feature of SVGD is not inevitable for
2

---

Under review as a conference paper at ICLR 2022
REGS since the approximation and expressive powers of deep neural networks is known to surpass those of kernel methods.
MCMC algorithms constructed from overdamped Langevin diffusion can be studied as discretization of the relative entropy gradient flow (Jordan et al. 1998). Based on the Suler--aruyama discretization of overdamped Langevin diffusion, unadjusted Langevin algorithm (ULA) (Roberts & Tweedie 1996) aims at generating samples from an approximation of the unnormalized target, but is biased for fixed step size. When Metropolis-Hastinng step is included, Metropolis-adjusted Langevin algorithm (MALA) (Roberts & Tweedie, 1996) is capable of correcting the bias, but leaves a large number of intermediate samples rejected. In REGS, one only needs to estimate the deterministic velocity fields of the relative entropy gradient flow, which differs from running ULA and MALA with randomness from diffusion processes. All particles produced by REGS are generated from an approximation of the target distribution.
Another line of work (Gao et al. 2019 2021) uses Wasserstein gradient flows of f-divergences for generative learning with samples from the underlying target distribution. In their work, evaluating velocity fields of gradient flows also boils down to estimating density ratios. However, our current problem is to sample from unnormalized target density. Furthermore, we propose a novel density ratio estimation procedure when the target distribution is only known up to a normalizing constant.
Notation
Let (2) be the space of Borel probability measures support space X' C Ed with a finite second moment, and let Pn (2t) be a subspace of P2(x)) whose measures are absolutely continuous w.r.t. the Lebesgue measure. All probability measures considered thereinafter assumed to belong to P2 (2)). To ease the notation, we use probability density functions such as q(x), px)), ) X X to express probability distributions in P2 Let (P2(),, W2) denote the metric space P2 (2') endowed with the 2-Wasserstein distance W2. which is referred to as the quadratic Wasserstein space. We use and Div to denote the gradient operator and the divergence operator, respectively.
2 PROBLEM FORMULATION
Consider an unnormalized probability density function u : X + 00 oo). where X' C Rd is the support of u. Suppose has an intractable normalizing constant Z fx u(x)xx < 20. Our goal is to generate random samples from the underlying distribution p E Pp (A'), whose probability density function is only known up to proportionality. i.e., p(x) = u(x) x E X. The basic idea is to gradually optimize samples from given distribution q E Pg (X) to approximate samples from p. where is easy to sample from q. Optimizing samples leads to functional optimization of distributions. We then introduce the classical relative entropy the functional optimization objective. The relative entropy, a.k.a., the Kullback- divergence, for q.P E PP (X) ratio, which is defined as
is the average logarithmic density
(al(l)))) = / ((x)) log (x()
dx.
(1)
(p(x)
holds that We (gllp) and We (gllp) 0 iff q(x) = =((x) a.e. X E X. Moreover, we denote the relative entropy functional as
FH = (+Ipp): (x) + 10, oc].
(2)
To sample from the unnormalized density u = pZ, we consider the functional minimization problem
min Flgl. (EP)(X)
(3)
where Flq] is always minimized the underlying target distribution p. i.c., qxx) = p(x) a.e. X. In a nutshell, problem (3) is an energy functional minimization problem in metric space. To minimize the energy functional F, it suffices to move along the corresponding gradient flow in metric space until the flow converges. For example, gradient flow in the Euclidean space refers to a curve whose tangent space contains the steepest descent direction of a given function. Analogously, a gradient flow in the space of probability measures means curve that points in the steepest descent direction of given energy functional. When equipped with the 2- Wasserstein distance, minimization of the energy functional J naturally corresponds to continuous path on the quadratic Wasserstein space of distributions. which is commonly known as a Wasserstein gradient flow of the relative entropy. We call this flow a relative entropy gradient flow for briefness.
3

---

Under review as a conference paper at ICLR 2022
3 RELATIVE ENTROPY GRADIENT FLOW
In this section, we briefly review the formulation of relative entropy gradient flow and its connections to differential equations. We consider the properties of gradient flows in the quadratic Wasserstein space (P2 (,)), W2). Recall that F in (2) is the relative entropy functional defined on (P2 (2)), W2).
One can show that a curve (2: (P4 (x)), W2) is relative entropy gradient flow of J if satisfies the continuity equation (Ambrosio et al. (2008), ) page 295 and Villani (008). page 631),
8F(ye)
Oeg: Div q:
ōqc
where g: (x) = g(1, x) evolves over time,
6Fleel
log 2: the first variation of the energy functional
8gg
F at qc. and V & is the Euclidean gradient of ōqc
6Fludl Here, we identify the gradient as the relative
dgg entropy gradient, which is defined by
V Fwg = =V 8F(4
V log
(5)
dge
p
Moreover, the relative entropy F dissipates along the relative entropy gradient flow (4 heo the rate (Ambrosio et al. (2008), page 295)
0.Fla: = E.g II Flae]l
(6)
Therefore, the relative entropy gradient flow {ge 122 eventually converges to the target distribution p as ( + oo. As pointed out in Ambrosio et al. (2008) (Page 175), under mild conditions the continuity equation (4) concerning f2: determines time-inhomogeneuus Markov process (Xeh20 that starts at a random particle Xo ~ qo and follows the particle evolution dynamics
dx, dt
vx())), Xe ~ qt t 0.
Note that the velocity fields
P t 0
V. = Vw, = V log
(8)
drive the evolution of the particle X, in the Fuclidian space, which results in the transport of q in (P)) W2). An important observation is that
P log
V1
V log
2 0.
(9)
q:
Therefore, the velocity fields do not involve the unknown normalizing constant Z. This is the key motivation for us use the relative entropy gradient flow in the proposed method.
4 SAMPLING AS PARTICLE EVOLUTION
As indicated by the energy dissipation of relative entropy F in (6). running the relative entropy gradient flow {ge 1220 dynamics can provide a nice approximate solution to the functional minimiza- tion problem 33) when time / is large enough. Therefore, to sample from the target distribution is appropriate to simulate the relative entropy gradient flow {4theeu with the time horizon T sufficiently large. A natural strategy is to discretize the particle evolution form of relative entropy gradient flow in 7 with forward Euler iterations Veque 2007) as follows,
X
Xx + sVx(X+), Xo qo. k = 1... K K-1 1.
(10)
with the velocity field at step le
log P q6
V& -V,
(11)
where & > O is a tunable small step size, K = [T/s] is the number of iterations and qu is the
corresponding discretized gradient flow step k: i.e., X₁ Ok- Combining the expressions in (10) and (T)). we have that the iterations progress according to
X = X + V P Xo ~ 9o, k = 0, 0,... K - 1.
(12)
Ok
4

---

Under review as a conference paper at ICLR 2022
In principle, it is necessary evaluate the velocity field va V log(p/ (k)pp//))))))) each iteration in 12 By L the velocity field the relative entropy gradient flow can be simplified to
V log Q
2= V
log - -
1,
(13)
q
where Zp is the given unnormalized density of the target distribution p. Then only the density Qs: remains unknown for evaluating the velocity field. Ideally, (k can be estimated by evolving a large number of particles {X However, direct estimation of O is difficult due the curse of dimensionality and the potential expensive computation cost for different ks. Our solution is to approximate the velocity field as a whole.
Assuming a nice approximation Va of the velocity field 13 is provided, then one can implement the following iterations for approximately sampling from 9к with no effort,
Xp+1= Xx + (Xx), Xo qo, 0,
(14)
Through the iterations above, we can collect Xs q qk: 1.2.... K. We will discuss approximation of the velocity field V2 = V log(u// gx) from the perspective of estimating the logarithmic density ratio log(u/qu) in the next section.
LOGARITHMIC DENSITY RATIO ESTIMATION AND THE RELATIVE ENTROPY GRADIENT SAMPLER
In this section, we first propose a novel estimation procedure of the logarithmic density ratio log(u/1) based on an unnormalized density u and random samples from q.
We use a model ratio R : X + [0, oo) to fit the true ratio R u/qq between a density g and an unnormalized density Let g
E be a differentiable and strictly convex function. A Bregman score (Daiid/ 2007 Gneiting & Raftery 2007 Kanamori & Sugiyama (014)) with the
base probability measure q E (2)) to measure the discrepancy between R. and Kig is defined by
[ o' (R(()))) w (X)
((X)
B(R) (8 ((R(X))
Exoo
where Pg (2') is an introduced and reference distribution for calculating the integral involving u. It should be easy to sample from w and the support of u should be included in the support w. Additionally, 2 (R) 2 B (RRig). where the equality holds iff R(x) Rag u4,u))a.e.
X.
In this work, we take g(x) log(x) is to satisfy the basic requirement of the Bregman score; (b) cancellation of the unknown normalizing constant Z of u. Simple calculation shows that B (R) can be written
We use this function for two reasons: (a) convexity, this
B(R)= = (Exx_[((X) - Ex-
Tu xu(X) wx))
log (((X))
(15)
Recall that the true density ratio Ri, can be factorized Rig u/q Z(p/q) Thus the numerical scale of the true density ratio Rt hinges on two factors, i.c., the normalizing constant Z of u and the standard density ratio p/g. Since numerical scales of these factors are difficult to determine in applications, the induced numerical instability can deteriorate the density ratio estimate. In order to prevent the density ratio estimation from such instability, we consider the model ratio R on the logarithmic scale. This will also release the nonnegative constraint on
byproduct.
From now on, we denote D uu = log( ( D = log (R) : X + R. Then B(D) can be rewritten as
u X w
B B(D) = Ex Ex~g (exe(((((X)))))) Ex
D(X)
(16)
can be shown that the logarithmic density ratio Dig is identifiable at the population level by minimizing (16) with respect to D. Theorem 1. For B(D) defined in (16,, we have Di,, E minp B(D). In addition, for any D with
Ex-e
wxX D(X) < x0, 2 (D) 2 B( Dr. with equality iff D(x) = Diy(x) (g, u)-u.e x E X.
vxX)
5

---

Under review as a conference paper at ICLR 2022
Algorithm 1: REGS: Relative entropy gradient sampler
Input: u 7p
// innormalized target: density // step size, maximum 1oop count // initial particles // reference distribution
step size 8 > 0, an integer K
X qu, i= 1,,, .,.,..,
w E (x))
k t- 0
while k K do
Vi ~ w, i = ,.,,...,
// reference samples
Dos E arg ming D Li=1 (Xi)) +0 Da(Y
// log density ratio
Vx(x) = Vi (x)
// velocity field
Xh+1 Xt + Vial (Ki), 1, .2.,
// update particles
k k+1
end
Output: X ūx 23 p, 1,2....
output particles
Based on Theorem we can estimate the unknown logarithmic density ratio D = llollu//// / (s)) with a deep neural network De with parameter & through the sample version of (16). Let ( (Xg] be i.i.d. samples from dı ~ qk and {YE)== be i.i.d. samples from reference distribution w. We solve the following deep nonparametric estimation problem via stochastic gradient descent (SGD) for Do
Do, E arg min B(Dx) = 1 Σ exp( (x+))
(Y w (YE)
(17)
Dg
With the logarithmic density ratio estimator Doc the velocity field vı in 13 can be approximately computed by Vk VD... By considering sampling particle evolution process discussed in
Section REGS updates the initial particles { with iterations 14 as follows:
= = X S
X
go, i 1,1,2.....
k = 0.1 K K-1 1.
(18)
We summarize the proposed REGS for sampling from unnormalized density in Algorithm
NUMERICAL EXPERIMENTS
We evaluate REGS large number of ID and 2D mixture distributions and test its stability in the high-dimensional setting with multivariate Gaussian distributions. We also use REGS to perform consider three existing
benchmark datasets. For comparison, methods including SVGD (Liu & Wang 2016), ULA (Roberts & Tweedie 1996) and MALA (Roberts) Tweedie 1996). All experiments are done using NVIDIA Tesla K80 GPU and common CPU computing resources. The neural network architecture, hyperparameter values, dataset descriptions, and additional experimental results are given in the appendix. The python code of REGS is available at httttps://qquu //tp:bbggttbbbb/// com htt///nnoon.ooooo /EGG
Bayesian logistic regression
6.1 MIXTURE DISTRIBUTIONS
We run REGS and SVGD, ULA and MALA to generate 2000 particles for mixtures of 2, 8 and Gaussians (see Scenarios 4, 5, in Appendix B), and 5000 particles for a mixture of 25 Gaussians (see Scenario 9 in Appendix B). The sampling qualities of these algorithms plots with density contours target mixture distributions. We classify all scatter points with labels according to the nearest mode, and plot the histograms of the label counts.
compared by scatter

---

Under review as a conference paper at ICLR 2022
(a) REGS
(b) SVGD
antii (d) MALA with 50 chains
(c) ULA with 50 chains
Figure 2: Mixtures of 8 Gaussians with equal weights: scatter plots and histograms of generated samples by (a) REGS, (b) SVGD, (c) ULA with 50 chains, and (d) MALA with 50 chains. From left right in each subfigure, the variance of Gaussians varies from o² 0.2 (first column), o² 0.1 (second column), o² = 0.05 (third column), to x² = 0.03 (fourth column).
Gaussian mixtures with equal weights Figure 2 shows the scatter plots and histograms of samples generated by (a) REGS, (b) SVGD, (c) ULA with 50 chains, and (d) MALA with 50 chains from mixtures of 8 Gaussians with equal weights. It shows that REGS able explore all the components in the mixture distribution nearly equally. However, SVGD is only able to find part of the modes, as indicated in Figures 2(b) Figures ((c))es[((c]ree) and ((d) show that MALA and ULA with 50 chains find all modes but with unequal weights, especially as the variance of each component decreases.
(a) REGS
(b) SVGD
(c) ULA with 50 chains
(d) MALA chains
Figure 3: Mixtures of 8 Gaussians with unequal weights: scatter plots and histograms of generated samples by (a) REGS, (b) SVGD, (c) ULA with 50 chains, and (d) MALA with 50 chains. From left to right in each subfigure, the variance of Gaussians varies from o² 0.2 (first column), 0.1 (second column), ṭ² 0.05 (third column), to o 0.03 (fourth column).
Gaussian mixtures with unequal weights Figure shows the seatter plots and histograms of samples generated by (a) REGS, (b) SVGD, (c) ULA with 50 chains, and (d) MALA with 50 chains
from mixtures of 8 Gaussians with unequal weights (1, 1, 1, 3, 3, 3)/6. Figure 3(a) shows that the samples generated by REGS have the correct weights. Figures 3(c) and 3(d indicate that ULA
7

---

Under review as a conference paper at ICLR 2022
Figure 4: Monte Carlo estimates of [[x(X) versus d for d-dimensional multivariate Gaussian distributions of X, For increasing from 10 300 with lag 10. From left to right, h(x) a²a, (a"x)² exp(a"z) z). and 10 cos(a"x + 1/2) with a Rd ||a|l2 = 1. The curves represent the estimates using the target samples ("true", " blue solid line) and the generated samples by REGS (red solid line), SVGD (green dash line). ULA_ 1: gray dotted line, MALA 1: pink dotted line, ULA_50 oringe dotted line, and MALA 50: orchid dotted line.
and MALA assign particles modes with incorrect weights. Moreover, the quality of the samples generated by SVGD, ULA and MALA deteriorates as the number of modes increases, while the performance of REGS remains stable. We also included the results from ULA and MALA with a single chain in Figures 7 and 10 in Appendix D, which show that these samplers have difficulty with multimodal distributions if only a single chain is used.
To further analyze the performance, we report the Monte Carlo estimates of E[h(X)] using test function h in Table 1 where h(x)) = axx, x, (a (aTx)² and 10 cos((a 1/2) with a E R ||a = 1, and X is distributed as various Gaussian mixtures with unequal weights. By comparing the Monte Carlo estimates of E[](X)] using the samplers with the values based on target samples, we see that REGS performs better and is
stable than SVGD, ULA and MALA, especially when h(z)= = 10 cos x + 1/2). We include additional numerical results including more seatter plots and histograms (Figure 6-10) and Monte Carlo estimates with equal wieights (Table 6) in Appendix D.
Table 1: Monte Carlo estimates of [X)] with four samplers for 2D mixtures of Gaussians random vectors X with unequal weights. "Target" denotes the Monte Carlo estimate with target samples. ULA_k A and MALA_K k denote the ULA and MALA with & chains, respectively.
Tayst B00K
ULA.T
LEX21
Tage
3V0D
TEA.1
MALLT
ULC.0
MALKCE
Taga
ULAT
0.1
-3.71
4147
07
-2.03
2.83
0.45
0.30
212
1.39
32.30
$:11
1.99
$.33
8.13
3.48
2.50
-1.1
-6.44
-6.35
-6.12
-2.81
-8.50
4.53
+2.33
3383
0.45
2.33
31
3.55
4.25
-6.33
Spariis
002
-1.30
-1.88
3.55
0U48
1133
DU66
-130
8.55
8.33
3,05
8.60
8.63
3.55
3.54
3.55
-316
1.65
-534
-578
-331
3.6
-3.21
-1.3
(340
2.87
1350
5.56
0.55
K34
1.45
4.55
0.31
1.33
0000
3.55
1550
8.5
3.10
8.65
3.55
3.45
1411
314
Mpra
1.00
1.0
211
056
3.55
357
204
5300
7.2
46.3
0.5
106
110
1991
8.0
1.05
116
000
790
133
110
1564
4.54
730
011
411
3.90
1118
6.2 MULTIVARIATE GAUSSIAN DISTRIBUTION
Let the target distribution be a d-dimensional Gaussian distribution with mean μrr = (1, 1,, 1) Rd and covariance matrix Σ E ddxd Σi, = |-33 with p = 0.7. We consider four test functions h(x). i.e., h(x)) a² (the first moment), h(x)) = (2 (the second moment), h(xx) exp(a">> (the moment generating function), and h(xx) = 10 cos(a"", Tx +- 1/2) with a E Rdd satisfying ||l] = 1. For reference, we provide the Monte Carlo estimates of E[h(X) using target samples. We compare REGS with SVGD. ULA 1, MALA 1. ULA 50, MALA 50 in Figure 4 the number of particles 5000 for each sampler, where ULA_ & and MALA 3: denote the ULA and MALA with le chians. For ULA and MALA, because of large variations of the estimates, we repeat the process 10 times and compute the average as the final estimate. Figure presents these Monte Carlo estimates as increases from 10 to 300 with step size 10. The logarithm of the estimated Clexp(a" X)] is shown. As shown in Figure the estimates using REGS and SVGD have smaller fluctuations than those using ULA and MALA, although all four methods can estimate Eja" X] and E x)² well. Moreover, the third and the fourth panels in Figure 4 show that REGS outperforms SVGD, ULA and MALA when h(x)) exp(a"2) 10 cos(( 1/2).
8

---

Under review as a conference paper at ICLR 2022
6.3 BAYESIAN LOGISTIC REGRESSION
We apply REGS to Bayesian logistic regression for binary classification on five datasets, including Banana, German, Image, Ringnorm, and Covertype. These datasets were analyzed in in/Liu & Wang (2016) and the first four datasets had been analyzed in [eershaan et al.] (2012). We consider a similar setting
that in (Liu & Wang 2016; Gershman al.) 2012), which assigns Gaussian prior (B(a) = =N x "1] to the regression coefficient / (including the intercept). We specify the prior of a as x (a) Gamma(1, 0.01). ). For comparison, we consider SVGD, ULA and MALA. The
inference is based on the posterior data).
These datasets are partitioned randomly into two parts, the training sets (80%) and the test sets (20%). We repeats the random partition 10 times. We evaluate the classification accuracy test data with 5000 particles from the posterior. T lists the averages and standard errors (in parenhteses) of test accuracy. From Table we can see that REGS is comparable with SVGD, ULA and MALA. For the Covertype dataset, MALA failed to converge, so no results from are included in Table 2
Table Averages and standard errors (in parenhteses) of classification accuracy on test data from five datasets, d: number of features, N sample size.
Averages of Accuracy (%)
datasets
d
N
REGS
SVGD
U..A1
MALAA
ULA.A 50
MALA 50
Banana
2
5300
(3.1)
(2.9)
(1.1
(1.9)
(1.1
55.2 (2.1)
German
20
1000
(2.2)
(1.6))
(1.8)
76.6 (2.2)
(2.6
Image
2086
(1.5
(1.7
(2.3
82.9 (2.3)
(2.3
82.8 (1.4())4)
Ringnorm
20
7400
(0.9)
(1.0)
(1.4)
(1.4)
71.4
Covertype
54
581012
(1.2)
75.6 (0.8)
74.1 (0.3)
7..2.0044
6.4 DISCUSSION OF THE EXPERIMENTAL RESULTS
The experimental results reported above and in the appendix indicate that REGS is capable of generating better quality samples than SVGD, ULA and MALA from Gaussian mixture distributions. Also, the results suggest that particles generated by REGS can cross valleys in the landscape of multimodal distribution even they are initialized in a different regions. An intuitive explanation is as follows. The movement of the REGS partieles is determined by the velocity field. If the velocity field not zero at a particle, the particle will continue to evolve towards the target distribution. Moreover, all particles interact with each other through the velocity field, which is beneficial in sampling from multimodal distributions. For ULA and MALA, there are no interactions among particles or incentives for particles to cross valleys between two modes, thus is more difficult for these methods to sample from multimodal distributions. A possible remedy is to use multiple chains as we did in the above experiments. To some extend, this alleviates the problem encountered in sampling from multimodal distributions. However, the success of this strategy depends on the initial samples being near the modes as well as having the correct proportions of the initial samples being close to each mode. In comparison, REGS uses principled way move particles from initial reference distribution to a multimodal distribution, albeit with higher computational cost.
CONCLUSION
We have introduced REGS, novel gradient flow based method for sampling from unnormalized distributions. Extensive numerical experiments demonstrate that REGS performs better than several existing popular sampling methods the setting of challenging multimodal mixture distributions. In future work, we hope to establish the convergence properties of REGS generated sampling distributions as the numbers of iterations and particles increase.
As with any sampling algorithms, there is a trade-off between sampling quality and computational efficiency. On one hand, as our numerical experiments demonstrate, REGS can generate samples with better quality than the three existing methods we considered in the challenging mixture model settings. On the other hand, REGS is computationally more expensive, as it involves neural network training in the iterations, compared with existing methods such as ULA and MALA that can be implemented more quickly. As computational power continues to increase rapidly, REGS can be a useful addition to the toolkit of sampling methods for multimodal distributions.

---

Under review as a conference paper at ICLR 2022
REFERENCES
Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient Flows: in Metric Spaces and in the Space of Probability Measures. Springer Science & Business Media, 2008.
Christophe Andrieu, Nando de Freitas, Arnaud Doucet, and Michael Jordan, An introduction to MCMC for machine learning. Machine Learning, 50(1):5-444 2003.
Matthew J Beal. Variational Algorithms for Approximate Bayesian Inference. PhD thesis, University College London, 2003.
David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859- 8778859977777-- 2017.
Steve Brooks, Andrew Gelman, Galin L., Jones, and Xiao-Li Meng. Handbook of Markov Chain Monte Carlo. CRC Press, 2011.
Changyou Chen, Ruiyi Zhang, Wenlin Wang, Bai Li, and Liqun Chen. A unified particle-opimmizatin framework for scalable bayesian sampling. In UAI, 2018.
Tianqi Chen, Emily Fox, and Carlos Guestrin. Stochastic gradient Hamiltonian Monte Carlo. In Eric P. Xing and Tony Jebara (eds.), Proceedings of the 31st International Conference on Machine Learning, volume 32 of Proceedings of Machine Learning Research, pp. 1693-1691. PMLR, 22-24 Jun 2014.
A Philip Dawid. The geometry of proper scoring rules. Annals of the Institute of Statistical Mathematics, 59(1):77-93, 2007.
Simon Duane, A.D. Kennedy, Brian J. Pendleton, and Duncan Roweth. Hybrid Monte Carlo. Physics Letters 195-2):216-222, 1987.
Andrew Duncan, Nikolas Nüsken, and Lukasz Szpruch. On the geometry of Stein variational gradient descent. arXiv preprint arXiv: 19/2:0092000 2019.
D B Dunson and E Johndrow. The Hastings algorithm at fifty. Biometrika, 107(1): 2019.
Yuan Gao, Yuling Jiao, Yang Wang. Yao Wang, Can Yang, and Shunkang Zhang. Deep generative learning via variational gradient flow. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 2093-2101. PMLR, 09-15 Jun 2019.
Yuan Gao, Jian Huang, Yuling Jiao, Jin Liu, Xiliang Lu, and Zhijian Yang. Deep generative learning with Euler particle transport. Proceedings of Machine Learning Research vol 145:1-33, 2021 2nd Annual Conference on Mathematical and Scientific Machine Learning, 2021.
S. Gershman, M. Hoffman, and D. Blei. Nonparametric variational inference. ICML, 2012.
Tilmann Gneiting and Adrian E Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association. 102(477):359-378 2007.
W. Keith Hastings. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1):97-109, 1970.
Matthew D. Hoffman and Andrew Gelman. The No-U-Turn Sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(47):1593-1623 2014.
Richard Jordan, David Kinderlehrer, and Felix Otto. The variational formulation of the Fokker-Planck equation. SIAM Journal on Mathematical Analysis, 29(1):1-17, 1998.
Takafumi Kanamori and Masashi Sugiyama. Statistical analysis of distance estimators with density differences and density ratios. Entropy, 16(2):921-942, 2014.
Anna Korba, Adil Salim, Michael Arbel, Giulia Luise, and Arthur Gretton. A non-asymptotic analysis for stein variational gradient descent. Advances in Neural Information Processing Systems, 33, 2020.
10

---

Under review as a conference paper at ICLR 2022
Randall J Le Veque. Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-state and Time- dependeen Problems, volume 98. SIAM, 2007.
Chang Liu, Jingwei Zhuo, Pengyu Cheng, Ruiyi Zhang, and Jun Zhu. Understanding and accelerating particle-based variational inference. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 4082-4092. PMLR, 09-15 Jun 2019a.
Chang Liu, Jingwei Zhuo, and Jun Zhu. Understanding MCMC dynamics flows on the Wasser- stein space. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 4003-4103. PMLR, 09-15 Jun 2019b.
Qiang Liu. Stein variational gradient descent as gradient flow. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, Vishwanathan, and Garnett (eds.), Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
Qiang Liu and Dilin Wang. Stein variational gradient descent: A general purpose bayesian inference algorithm. In D. Lee, M. Sugiyama, U. Luxburg, Guyon, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 29, Curran Associates, Inc., 2016.
Jianfeng Lu, Yulong Lu, and James Nolen. Scaling limit of the Stein variational gradient descent: The mean field regime. SIAM Journal on Mathematical Analysis, 51(2) 6711, 2019.
Nicholas Metropolis, Arianna Rosenbluth, Marshall N Rosenbluth, Augusta Teller, and Edward Teller. Equation of state calculations by fast computing machines. The journal of Chemical Physics, 21(6):1087-1092, 1953.
Radford M. Neal. MCMC Using Hamiltonian Dynamics, chapter 5. CRC Press, 2011.
Gareth O Roberts and Osnat Stramer. Langevin diffusions and Metropolis-Hastings algorithms. Methodology and computing applied probability, 4(4):337-357, 2002.
Gareth O. Roberts and Richard L. Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 2(4):341 - 363, 1996.
Adil Salim, Anna Korba, and Giulia Luise. The wasserstein proximal gradient algorithm. arXiv preprint XXXx:20022200033 2020.
Adil Salim, Lukang Sun, and Peter Richtárik. Complexity analysis of stein variational gradient descent under talagrand's inequality tl. arXiv preprint 2xXi:2106.000770 2021.
Luke Tierney. Markov Chains for exploring posterior distributions. The Annals of Statistics, 22(4): 1701-1728, 1994.
Cédric Villani. Optimal Transport: Old and New, volume 338. Springer Science Business Media, 2008.
Martin J. Wainwright and Michael I. Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning. 111):1-305, 2008
Max Welling and Yee Whye Teh. Bayesian learning via stochastic gradient Langevin dynamics. In Proceedings of the 28th international conference machine learning, ICMI.'1, 11, pp. 681-688. ACM, 2011.
Michael Zhu, Chang Liu, and Jun Zhu. Variance reduction and quasi-Newton for particle-based variational inference. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 11576-11587. PMLR, 13-18 Jul 2020.
Jingwei Zhuo, Chang Liu, Jiaxin Shi, Jun Zhu, Ning Chen, and Bo Zhang. Message passing Stein variational gradient descent. In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 6018-6027. PMLR, 10-15 Jul 2018.

---

Under review as a conference paper at ICLR 2022
A APPENDIX
In this appendix, we prove Theorem 1, give detailed description of the models in the numerical experiments, the neural network architecture used in implementing REGS, and additional numerical results.
B PROOF OF THEOREM
Proof. By the definition of Bregman score ((R) between u and is easy to check
Rig E arg min B(R), ft
B(R) 2 ((Rgg) with equality iff R(x) Rx)((x) u.. X.
Since D=g = Hog(H)))) ; D = log( R). B( 2 (D), when Ex
1 D(X)
c0, we have
D E arg min B(D), D
B(D) B(Dig) with equality iff D(x) (x) (1, X X.
C GAUSSIAN MIXTURE DISTRIBUTIONS
Let the density of d-dimensional multivariate Gaussian distribution with mean μ Rd and covariance matrix 2 E BRd be
1 (x; (2x) det(E)) exp - 1- (2: ~) m)
We consider several scenarios below for the unnormalized density function u(x).
Scenario 1. 2Gaussians_la IdI. One-dimensional mixture of 2 Gaussians,
1 fx; p1, of) + 5 f(= (2;
2
u(x) -
where = 1, = 2, and of 0.25, of = 2;
Scenario 2. (Gaussians_Id2. One-dimensional mixture of 2 Gaussians,
1 u(x) ( f(; (1; f(z; 12, (2, oz),
2
where #1 3, 12 3, and of 0.25, o3 = 2: Scenario 2Gaussians_Idt One-dimensional mixture of 2 Gaussians,
1 u(x) = fx; fx; R1, of) + f(x: p2, (1; of),
where #1 3, μl 3, and of 0.03, of 0.03; Scenario 4. 2Gaussians. Two-dimensional mixture of Gaussians,
u(z) f(:; M1 En) ff (2;
U2, Ez),
where μi = (r. === (-r,0) and Σ, Σソ =
Scenario 5. 8Gaussians. Two-dimensional mixture of Gaussians,
8 u(x) Σ
f(r
j=1
where ny r(sin(2() - 1)m/8), cos (20 1)T T8) and E, o²I for 1,... ,8
12

---

Under review as a conference paper at ICLR 2022
Scenario 6, 9Gaussians. Two-dimensional mixture of 9 Gaussians, 3
u(x)
Σ faz Pikk 15kk) j=
where Myg = 4(j -, , 22 for = 2, and 1,2,3.
Scenario 16Gaussians_le Ic. Two-dimensional mixture of 16 Gaussians.
16
u(x)= Σ
f(m
==1
where ,4 = 4(sin(2() - 1,... .6.
- /)///))) and Ej 0.031 for j =
Scenario 8. 16Gaussians 2c. Two-dimensional mixture of 16 Gaussians,
fz; M.. 21, + Σ fa: 2s), u(x) Σ 5=1
where Hy 1)m/8), cos(2(j - ((sin(2( nT 8). cos cos (2(j - 1)T/8))² Pk = 2(sin(2( and Σn E4 0.03I for j = 1,,,,, ,8,,,,ad and k .....
Scenario 9. 25Gaussians. Two-dimensional mixture of 25 Gaussians,
5 5
u(x) ΣΣ
(a; Pik, Estt:
j=1
where Myk = 2(j 33 - ,,
I for = 1, and 1,,..,, 1,.
Scenario 10. 49Gaussians. Two-dimensional mixture 49 Gaussians,
7
u(x) ΣΣ
fæ; Pyk Ezk).
11=1
where Mys = 1 - 4,4 T)T, Enx = 0.03I for I,... 77 and k 1,... ,7
Scenario 11. Gaussians. Two-dimensional mixture 81 Gaussians,
u(x) ΣΣ fl; Myg. Eju).
j=1==1
where MyR = (j k, E,, 0.031 for 1.... and k 1,... ,9.
Scenario 12. Icircle. Let μ4 4(cos sin (i)n/NN)TT. 0, 1,..... N with N 400. Consider three noise to be added to each point μú, including the uniform distribu- tion U (H4 1/30) on a disc with center at μ4 and radius 1/30, Gaussian distribution N14,0003II, and mixed these two distributions.
Scenario 13. 2circles.
= 2(cos(2im/N) /N), sin (i) and 4( ((cos(2iz/N). (iniN))ssin((iii/)e sin(2it/N)) i 0, 1. N with N 200. Consider three noise to be added to each point μąs including uniform distribution U(44s 1/30) on dise with center at μ4s and radius 1/30, Gaussian distribution N(us,, 0.03I), and mixed these two distributions, k- ==1 1.2.
Let
Scenario 14. Ispiral. Let μ4 = 4im N cos(4it/NV sin (in(4i//N))" ( i= 0, o.......
- N with = 400. Consider three noise to be added each point f4i including the uniform distribu- tion U (p, 1 1/30) on a disc with center at and radius 1/30, Gaussian distribution N(u, 0.03I), and mixed these two distributions.
Scenario 15. 2spirals.
(cos(3im siin/NNsssiiiiiiis /N))T and P22 0,1....,45 N 1 with N = 200. Con- sider three noise to be added each point Mąs including uniform distribution
Let Phi
N' (cos(3in//a) /N), sin (3iz/N)),, 3im
U Up 1/30) on a dise with center at Pu and radius 1/30, Gaussian distribution N N(us; 0.03T), and mixed these two distributions, k == 1,2.
13

---

Under review as a conference paper at ICLR 2022
Scenario 16. moons. Let M1 (8i/N 6, 6.44imm si and μuu (8i/N - 2, sin(in/ (
== 0,1, 1,,,,..... .N with = 200. Consider three noise to be added to each point pki including uniform distribution U ('(ua 1/30) on disc with center Mkx and radius 1/30, Gaussian distribution N (/ 0.03T), and mixed these two distributions, 1,
D EXPERIMENTAL SETTING
D.1 HYPERPARAMETER
We burn in the first 1000 particles for ULA and MALA in all experiments. We initialize the particles in SVGD, ULA and MALA as zeros or random samples from Gaussians. We provide the step size settings in Table 3 For SVGD, we use RBF kernel k(x, x') = exp(-
x²12) and set the bandwidth as h = med² log n,, where med the median pairwise distances between the particles {ri))=1- We set the learning rate of neural networks 5e-4 for density ratio estimation in REGS. The network structures are presented Table The initial particles REGS are sampled from Gaussian distributions.
Table 3: Step size settings for REGS, SVGD, ULA and MALA. "BGIR" denotes four datasets including Banana, German, Image and Ringnorm.
Methods
2Gaussians
8Gaussians
9Gaussians
25Gaussians
BGIR
Covertype
REGS
Se-4
5e-4
Se-4
5e-4
2e-3
2e-3
SVGD
2e-2
2e-2
2e-2
2e-2
Se-2
5e-2
ULA
2e-2
5e-2
le-1
5e-2
le-3
lc.4
MALA
2e-1
Se-1
5e-1
le-3
D.2 NEURAL NETWORK ARCHITECTURE
Table 4: Neural network architecture for log-densigg ratio estimation: feedforward neural networks with eiuhh-width hidden layers and Leaky ReLU activation. Depth for 2Gaussians Idl, 2Gaussians _d2, 2Gaussian__ia ld3, and 2Gaussians. = 4 for 8Gaussians, 9Gaussians, Icircle, 2circles, Ispiral, 2spirals, and moons. 6 for 16Gaussians Ic, 16Gaussians 2c, 25Gaussians, 49Gaussians, and 87 Gaussians.
Layer
Details
Output size
Linear, LeakyReLU (0.2)
128
Linear
D.3 DATASETS
Covertype: binary html Banana, German, Image, Ringnorm: http:://theov/ cmp. uea aac. uk/
www. www.:ssie..nt.ed..ee... http://wwhaaattuddd.///////////a
mat vatlab/ae aal.d//aaullll fa..lat..t.. html
E ADDITIONAL NUMERICAL RESULTS
14

---

Under review as a conference paper at ICLR 2022
Table 5: Monte Carlo estimates of [X)]. Here h(x) a²x, (a"x)² and cos(( T= 1/2) with E R², ||a||2 = 1. "true" denotes the Monte Carlo estimate with target samples.
1((ar
a'x
h(())
(=x==)²²
h((r) cos( (a' + 1/2)
Distributions
REGS
REGS
true
REGS
2Gaussians Idl
-0.9886
-0.9887
3.7105
3.7093
0,46838
0.5418
2Gaussions_Ide ld2
0.99972
-1.1630
9.0032
9.0297
-8,3053
-8.2956
2Gaussians Id3
0.96022
-0.728
9.8727
9.8016
-6.3855
-6.1266
2Gaussians
+0.0019
0.0448
8.0243
8.0225
-8.2351
-8,2513
8Gaussians
-0.0021
-0.021
8.0118
8.0276
-3,3925
-3,3727
9Gaussians
0.0003
0.0100
10.6771
10,8063
0,7773
0.7532
16Gaussians Jc
+0.0008
-0.0077
8.0269
8.02889
-3,4393
-3,4225
16Gaussians _c
+0,0006
-0.0060
5.3544
5.2780
-1,4260
-1.877
25 Gaussians
0.0013
-0.0141
8,0276
8.0227
0.1248
0.1189
49Gussianss
0.0006
0.0001
9.0285
9.0215
0.2043
0.1990
8 Gaussians
-0.0011
-0.0009
15.0280
15.0261
0.4126
0.4196
8Gaussians 5=5
-0.0011
-0.0037
12.5291
12.5356
-1,2160
-1.2214
8Gaussians 1010
0.0021
-0.0030
50.0278
49.4783
3,3872
3,4238
&Gaussians 1=15
0.0003
-0.0306
112.5023
111.4456
-1.1116
-1.0930
Gaaussians o² 0.01
-0.0022
-0.0009
0.5119
0.5082
6,6390
6.6504
2Gusssians o² 0.005
-0,0012
0.0027
0.5019
0,4976
6.6694
6.6718
2aasssians 2² 0.0001
0.0005
0.0034
0.5043
0.5049
6.6558
6.6406
Table 6: Monte Carlo estimates of [[) by four samplers for 2D mixtures of Gaussians of X
with equal weights. Here h(x)) = a=xx, (a (x=x)² x²² or co 1/2) with a E IR,2, ||ll|/ "true"
denotes the Monte Carlo estimate with target samples. "ULA A² and "MALA k" denote the ULA and MALA with chains, respectively.
20CS
EVCD
1..A.
MALA.
(LA.5)
MALA.95
REGS
SVGD
ULA.
MALAL
LEA.55
MALA
SVED
DEA. MALA
IILA5
MALA
Spartts
1.3
232
1.30
-117
122
$111
13.0
2722
2.1
$.:1
234
137
116
125
411
-.11
12.0
0.85
5.00
0.82
6/66
1.42
-3.5
2.89
212
2.10
3.19
5.05
1.30
135
-332
5.57
4.30
15:5
-1.01
0.43
0.10
-1.41
1.42
8:31
43.25
232
230
2.11
2384
2.00
530
8.11
1,42
1.35
1.12
3.35
-3.35
-5.59
+2.5
Spassiba
11.2
11:11
141
-2.30
-228
13:1
040
3.5
835
0985
(301
6.68
010
3.53
3.53
0.11
0.35
5.6
8.66
831
205
830
3.30
610
5.3
-338
135
0.5
635
3.50
4.33
40.0
-8.0
4.59
4338
137
6.5
-225
111
1.91
0.90
016
8:15
-335
1U/5
5.11
8.99
211
742
8.99
5.00
150
016
341
8.99
-315
011
4.80
4880
035
4.30
LXT
4.7
337
0.0
15

---

Under review as a conference paper at ICLR 2022
REGS
(b) SVGD
(c) ULA with one chain
dd) MALA with one chain
(e) ÚLA with chains
MALA with 50 chains
Figure 5: Scatter plots with contours 1000 generated samples from unnormalized mixtures Gaussians with equal weights by (a) REGS, (b) SVGD, (c) ULA and (d) MALA with chain, (e) ULA and (f) MALA with 50 chains. From left to right in each subfigure, the variance of Gaussians varies from o² 0.2 (first column), o² 0.1 (second column), d² 0.05 (third column), to x2 = 0.03 (fourth column).
16

---

Under review as a conference paper at ICLR 2022
REGS
(b) SVGD
(c) ULA with one chain
(d) MALA with one chain
annll iI) MALA with 50 chains
(e) ULA with 50 chains
Figure 6: Scatter plots with contours 2000 generated samples from unnormalized mixtures Gaussians with equal weights by (a) REGS, (b) SVGD, (c) ULA and (d) MALA with chain, (e) ULA and (f) MALA with 50 chains. From left to right in each subfigure, the variance of Gaussians varies from o² 0.2 (first column), o² 0.1 (second column), d² 0.05 (third column), to x2 = 0.03 (fourth column).
17

---

Under review as a conference paper at ICLR 2022
REGS
(b) SVGD
(c) ULA with one chain
(d) MALA with one chain
(e) ULA with 50 chains
(I) MALA with 50 chains
Figure 7: Scatter plots with contours of 5000 generated samples from unnormalized mixtures of 25 Gaussians with equal weights by (a) REGS, (b) SVGD, (c) ULA and (d) MALA with chain, (e) ULA and (f) MALA with 50 chains. From left to right in each subfigure, the variance of Gaussians varies from o² 0.2 (first column), o² 0.1 (second column), d² 0.05 (third column), to x2 = 0.03 (fourth column).
18

---

Under review as a conference paper at ICLR 2022
REGS
(b) SVGD
(c) ULA with one chain
dd MALA with one chain
(e) ULA with chains
MALA with 50 chains
Figure 8: Scatter plots with contours of 1000 generated samples from unnormalized mixtures of 2 Gaussians with unequal weights (0.75, 0.25) by (a) REGS, (b) SVGD, (c) ULA and (d) MALA with one chain, (e) ULA and (f) MALA with 50 chains. From left to right in each subfigure, the variance of Gaussians varies from o² 0.2 (first column), a² 0.1 (second column), o² 0.05 (third column), to 0.03 (fourth column).
19

---

Under review as a conference paper at ICLR 2022
REGS
(b) SVGD
(c) ULA with one chain
d MALA with one chain
(e) ULA with 50 chains
(1) MALA with 50 chains
Figure 9: Scatter plots with contours 2000 generated samples from unnormalized mixtures
Gaussians with unequal weights (1, 1, 1, 3, 3, 3) 16 by (a) REGS, (b) SVGD, (c) ULA and (d) MALA with one chain, (e) ULA and (f) MALA with 50 chains. From left to right in each subfigure, the variance of Gaussians varies from o² = 0.2 (first column), の² 0.1 (second column), σ² 0.05 (third column), to 0.03 (fourth column).
20

---

Under review as a conference paper at ICLR 2022
(a) REGS
(b) SVGD
(c) ULA with one chain
(d) MALA with one chain
(c) ULA with 50 chains
MALA with 50 chains
Figure 10: Scatter plots with contours of 5000 generated samples from unnormalized mixtures of 25 Gaussians with unequal weight by (a) REGS, (b) SVGD, (c) ULA and (d) MALA with chain, (c) ULA and (f) MALA with 50 chains, where each the first 12 components has weight 1/51, and each of the rest has weight 3/51. From left right in each subfigure, the variance of Gaussians varies from 6² = 0.2 (first column), d² = 0.1 (second column), o² = 0.05 (third column), to o2 = 0.03 (fourth column).
Figure 11: KDE plots for 1D mixtures of Gaussians. Green lines stand for target samples and pink areas represent generated samples by REGS. From left to right, the means and variances of the components changed and the unnormalized densities are given in Scenarios 1, in Appendix B.
21

---

Under review as a conference paper at ICLR 2022
Figure 12: KDE plots of target samples (first row) and generated samples (second row) for two- dimensional mixtures of Gaussians with variance 0.03. The target samples are from unnormalized density functions u(x) of mixtures of 2 Gaussians in Scenario 4, 8 Gaussians in Scenario 5 and Gaussians in Scenario 6.
Figure 13: KDE plots target samples (first row) and generated samples (second row) for 2D
mixtures of Gaussians with component variance 0.03. The corresponding unnormalized densities presented in Scenarios 7, 8, in Appendix B.
22

---

Under review as a conference paper at ICLR 2022
-3
-3
-3
Figure 14: Scatter plots of initial samples (first row), generated samples (second row) for two- dimensional mixtures of 2 Gaussians, and scatter plots of target samples (last row). The target samples are from unnormalized density functions u(x)) of mixtures of Gaussians with variance x² = 0.01 (left column), o² = 0.005 (middle column) and o² 0.001 (right column).
23

---

Under review as a conference paper at ICLR 2022
33
-3
3
3
-3
Figure 15: Scatter plots of initial samples (first row), generated samples (second row) for two- dimensional mixtures of multiple Gaussians with variance o² 0.03, and scatter plots of target samples (last row). The target samples are from unnormalized density functions u(x) of mixtures of 25 Gaussians (left column), 49 Gaussians (middle column) and 81 Gaussians (right column).
24

---

Under review as a conference paper at ICLR 2022
Figure 16: Scatter plots of initial samples (first row), generated samples (second row) for two- dimensional mixtures of Gaussians with variance o² = 0.03 and varying radius, and scatter plots of target samples (last row). The target samples are from unnormalized density functions u(x) of mixtures of Gaussians with radius being (left column), 10 (middle column) and 15 (right column).
25

---

Under review as a conference paper at ICLR 2022
(a) Target samples (first row) and generated samples (second row) from the uniform distribution
(b) Target samples (first row) and generated samples (sccond row) from Gaussian distribution
(c) Target samples (first row) and generated samples (second row) from mixed Gaussian and uniform distributions
Figure 17: Scatter plots from left to right are one circle (/circle, Scenario 12), two circles (2circles, Scenario 13), one spiral (Ispiral, Scenario 14), two spirals (2spirals, Scenario 15), and moons (moons, Scenario 16).
26


Reviews:
[{"review_id": "zfqaImku6o", "reviewer": "Reviewer_dwfT", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "Wasserstein gradient flow has proved to be a useful tool for sampling from an unnormalized distribution. Section 2 to 4 of this work follow the standard derivation of the work along this research line and well explain the particle evolution strategy. Since the underlying velocity field of the Wasserstein gradient flow requires the access to the variable distribution $q_t$ which is in general not available, a key step in methods along this research line is to estimate such a quantity. To estimate such a quantity, this work proposes to estimate the log density ratio between the potential function and the variable distribution $q_t$ by minimizing the Bregman score which is described in Section 5. \nHowever, I find Section 5 difficult to understand. It would be vary helpful if the authors could explain the intuition of the Bregman score. In fact, I think there should be an individual section in the preliminary that describes the Bregman score and all the statements below equation (15) so that the reader can follow this very important step. \n\nI think section 5 is the part that differs this work from previous work like [1] and is where the novelty of this paper lies. It need to be very clearly explained.\n\n\n[1] Degond, Pierre, and Francisco-José Mustieles. \"A deterministic approximation of diffusion equations using particles.\" SIAM Journal on Scientific and Statistical Computing 11, no. 2 (1990): 293-310.", "summary_of_the_paper": "This paper considers the problem of sampling from an unnormalized distribution. The unnormalized target distribution can be regarded as a stationary point of the Wasserstein gradient flow of the corresponding relative entropy functional, which can be equivalently identified from a microscopic perspective by defining a time-varying velocity field of the particles. While the exact time-varying velocity field is not exactly available, the authors propose to estimate such a quantity by approximating the corresponding logarithmic density ratio through minimizing the Bregman score. Such an approximation requires only samples from the variable distribution which can obtain by simulating particles following the estimated velocity field.", "summary_of_the_review": "This paper leverage the microscopic equivalence of the Wasserstein gradient flow of the relative entropy to sample from an unnormalized distribution, but the derivation of the key step in the proposed approach is not well explained. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "dxzNIx6l0K5", "reviewer": "Reviewer_g4HT", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "The paper appears to miss the fact that the 2-Wasserstein gradient flow for the relative entropy defines a Markov process, which is exactly the Langevin dynamics; this can be seen by comparing Eq. (4) to the Fokker–Planck equation for an Ito diffusion (e.g., Eq. (4.1) in [1]; see also section 3.5 in [2] for a relevant discussion from the ML literature). Indeed, the Fokker–Planck equation determines the diffusion uniquely because the differential operator in the Fokker–Planck equation is the adjoint of the infinitesimal generator of the diffusion. Thus, the proposed algorithm is a rather unorthodox way of approximating a Langevin distribution. \n\nThe paper makes the said approximation more difficult than it should be by using a particle approximation that requires estimating density ratios, a notorious tricky problem. The algorithm ends up having some ability to handle multimodality because the use of weighted particles and density ratio estimation allows the algorithm to effectively compute the relative volumes of different modes of the distribution. The proposed method for estimating the density ratios appears to be the same as Geyer’s reverse logistic regression method [3], with a neural net replacing the inner product. Thus, I expect similar results could be obtained more directly, and with only a single round of volume estimation, by (1) running MALA many times and (2) then estimating the relative volume of each chain (note there are numerous other methods for doing this other than reverse logistic regression). Such an approach should work well on the kinds of low-dimensional examples considered in the numerical experiments. \n\nFurther issues arise in the experimental evaluation. First, the experiments seem to show very similar performance to MALA, all within the standard errors when provided (e.g., in Table 2). Second, I’m concerned about the quality of the MALA implementation, for which code was not included. The lack of convergence in one example suggests MALA was not run with appropriate step size adaptation targeting the optimal acceptance rate [4,5], as is standard in the literature. If so, then the comparison is not appropriate. Third, for a fair comparison, the MALA chains should also be reweighted based on volume estimates for each chain, as described above. \n\nIs it possible there are some gains from using the proposed method on multimodal distribution? Yes. But I remain skeptical. Moreover, if the goal is prediction, I expect combining MCMC with stacking will be more effective [6,7]. \n\n\n[1] Pavliotis, G. A. Stochastic Processes and Applications. (Springer, 2014).\n\n[2] Liu, Q. Stein Variational Gradient Descent as Gradient Flow. In NeurIPS (2017).\n\n[3] Geyer, C. Estimating normalizing constants and reweighting mixtures. Technical Report (1994). \n\n[4] Roberts, G. O. & Rosenthal, J. S. Optimal scaling of discrete approximations to Langevin diffusions. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 60, 255–268 (1998).\n\n[5] Roberts, G. O. & Rosenthal, J. S. Optimal scaling for various Metropolis-Hastings algorithms. Statistical Science 16, 351–367 (2001).\n\n[6] Yao, Y., Vehtari, A. & Gelman, A. Stacking for Non-mixing Bayesian Computations: The Curse and Blessing of Multimodal Posteriors. arXiv.org, arXiv:2006.12335 (2020).\n\n[7] Yao, Y., Vehtari, A., Simpson, D. & Gelman, A. Using Stacking to Average Bayesian Predictive Distributions. Bayesian Analysis 13, 917–1007 (2017).", "summary_of_the_paper": "The paper suggests a method for approximating the 2-Wasserstein gradient flow for the relative entropy. The proposed particle-based method uses a neural network function approximation-based approach to estimating the necessary density ratios. Experiments verify reasonable performance compared to MALA and ULA.  ", "summary_of_the_review": "The paper seems to have a fundamental misunderstanding of the Wasserstein gradient flow for the relative entropy, and the experimental evaluations may not be appropriate. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "aYbIHqkaSk1", "reviewer": "Reviewer_opqD", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "I enjoyed reading the paper. It is well written, the motivation is clear and it is easy to follow the main ideas.\n\nHowever, I find it hard to assess the actual contribution of the paper.\nOn one hand, while the proposed algorithm makes sense, there is no guarantee in the paper, either for the sampling accuracy or even for the fact that the algorithm will converge to the target measure. For example, is there any guarantee that the discretized flow does not add bias to the obtained measure?\nThere are many tunable parameters in the algorithms, $s$ the discretization step, $K$, the time horizon, $n$, the number of particles. What is the interplay between those parameters, given some distribution, how should I choose those?\nI would have expected to see a bit more of the underlying theory behind this algorithm.\n\nOn the other hand, from the perspective of actual results, I find that the numerical experiments are somewhat restricted and artificial. Coupled with the computational overhead, it's not clear to me when one will actually prefer to use the new algorith,\n", "summary_of_the_paper": "The paper proposes a novel way to sample from unnormalized distributions. This is helpful when calculating or estimating the normalizing constant is untractable.\n\nThe main idea is to track the gradient flow of the relative entropy in the Wasserstein space of probability distributions. It is known that the flow converges to the target distribution and the paper introduces a variational characterization of the discretized steps. The main benefit of this characterization is that it bypasses the need to know the normalizing constant as well as being amenable to estimation by using a combined particle evolution.\n\nThe benefits of the new algorithm are demonstrated through several numerical simulations.", "summary_of_the_review": "The idea is elegant interesting but the paper lacks evidence for its usefulness, both from theoretical and applied perspectives.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "1: The main claims of the paper are incorrect or not at all supported by theory or empirical results.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "_l7F_PfUPhk", "reviewer": "Reviewer_r2Jt", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "Overall, the paper is well written: contributions are clearly stated, relation to previous is presented, the proposed method is well explained and a somewhat extended comparative numerical evaluation of the proposed method is given.\n\nI find the overall approach interesting - the difficulty of sampling being cast as a density ratio estimation problem. This new problem, however, is not easy to solve and your approach of using a neural network to estimate the density ratio seems to work, at least in the considered examples.\n\nThat being said, I do have some issues with the paper. You mention in the conclusion that you hope to establish the convergence properties of the proposed method. I understand that it is not trivial to establish convergence. However, not presenting at least some intuition about the convergence of the algorithm is a strong drawback. The way I see it is that there are two sources of error that could hamper convergence: the discretization error of the numerical implementation of the continuous gradient flow and the approximation error of the density ratio. I didn't find anything about either in the paper. If not a proof, at least some intuition about how they affect the results, about how they interact, etc.\n\nIn the numerical experiments section, you present a fair amount of examples which show that the proposed method is capable of outperforming the competing algorithms. The results are interesting, however, there are no results to show how the performances of the proposed method vary with the different parameters, notably the number of considered particles, and the choice of distribution w. Such results would cast some light into the inner workings of the proposed method and would be useful for anyone interested in using it.\n\nIn section 6.4 you do mention that the improved performances of the proposed method come with a higher computational cost. However, you do not perform any analysis of the trade-off between computation time and performances. It would have also been interesting to compare the performances of the different algorithms for the same computational budget.\n\nAnother aspect that struck out to me is the choice of competing algorithms. The algorithms that you choose as competitors are valid, however their choice is questionable. My first remark was why didn't you choose the SMC algorithm as a competitor? It also uses particles in order to estimate the unnormalized target distribution. Also, HMC could have also been considered.\n\nI spotted some typos here and there, for example chians instead of chains on page 8 in the bottom paragraph \"... denote the ULA and MALA with k chians\", repeats instead of repeat on page 9 \"We repeats the random partition 10 times.\". \n\nAnother small issue is with figure 4, it's hardly readable. I understand that there is a limit on the page count, however, that's not a justification for having figures that are hard to read. More so, as there are some redundancies in the text that could have been eliminated, for ex. equations (11) and (13) are the same, is the presence of both necessary for understanding the idea that is presented in section 4?", "summary_of_the_paper": "The paper addresses the issue of sampling from an unnormalized distribution. The sampling problem is cast as the numerical simulation of the gradient flow associated with the KL divergence between the target unnormalized distribution and the approximating distribution. The challenging part is to estimate the density ratio that appears in the gradient term. The authors propose to use a deep neural network to estimate the density ratio. Numerical results show the usefulness of the proposed method.", "summary_of_the_review": "I find the approach interesting. However, there are some issue with the paper as it is, both theoretical and empirical. From a theoretical point of view, there is no discussion about conditions for convergence of the algorithm. From an empirical point of view, the numerical experiments are not complete enough, with respect to the comparison analysis that is carried out, but also with respect to compensating missing theoretical analysis. Overall, the paper is interesting, but in its current form is not ripe enough for publication.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes a sampling technique for unnormalized distributions. The main idea is to gradually transform particles by following the gradient flow of the relative entropy in the Wasserstein space of probability distributions. The paper tackles an important problem and provides an interesting new perspective. However, even putting aside the concerns on the theoretical analysis raised by the reviewers, the experimental evaluations does not seem sufficient to demonstrate the benefits of the proposed approach.
Decision: Reject