summaries: 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.   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. The paper proposes a novel way to sample from unnormalized distributions. This is helpful when calculating or estimating the normalizing constant is untractable.  The 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.  The benefits of the new algorithm are demonstrated through several numerical simulations. 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.
abstractText: 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 mixture distributions and Bayesian logistic regression on real datasets demonstrate that the REGS outperforms the state-of-the-art sampling methods included in the comparison.
references.authors: Luigi Ambrosio Nicola Gigli Giuseppe Savaré
references.title: Gradient Flows: in Metric Spaces and in the Space of Probability Measures
references.venue: Springer Science & Business Media,
references.authors: Christophe Andrieu Nando de Freitas Arnaud Doucet Michael I. Jordan
references.title: An introduction to MCMC for machine learning
references.venue: Machine Learning,
references.authors: Matthew J Beal
references.title: Variational Algorithms for Approximate Bayesian Inference
references.venue: PhD thesis,
references.authors: David M. Blei Alp Kucukelbir Jon D. McAuliffe
references.title: Variational inference: A review for statisticians
references.venue: Journal of the American Statistical Association,
references.authors: Steve Brooks Andrew Gelman Galin L. Jones Xiao-Li Meng
references.title: Handbook of Markov Chain Monte Carlo
references.authors: Changyou Chen Ruiyi Zhang Wenlin Wang Bai Li Liqun Chen
references.title: A unified particle-optimization framework for scalable bayesian sampling
references.authors: Tianqi Chen Emily Fox Carlos Guestrin
references.title: Stochastic gradient Hamiltonian Monte Carlo
references.venue: Proceedings of the 31st International Conference on Machine Learning,
references.authors: A Philip Dawid
references.title: The geometry of proper scoring rules
references.venue: Annals of the Institute of Statistical Mathematics,
references.authors: Simon Duane A.D. Kennedy Brian J. Pendleton Duncan Roweth
references.title: Hybrid Monte Carlo
references.venue: Physics Letters B,
references.authors: Andrew Duncan Nikolas Nüsken Lukasz Szpruch
references.title: On the geometry of Stein variational gradient descent
references.venue: arXiv preprint arXiv:1912.00894,
references.authors: Yuan Gao Yuling Jiao Yang Wang Yao Wang Can Yang Shunkang Zhang
references.title: Deep generative learning via variational gradient flow
references.venue: Proceedings of the 36th International Conference on Machine Learning,
references.authors: Yuan Gao Jian Huang Yuling Jiao Jin Liu Xiliang Lu Zhijian Yang
references.title: Deep generative learning with Euler particle transport
references.venue: In Proceedings of Machine Learning Research
references.authors: S. Gershman M. Hoffman D. Blei
references.title: Nonparametric variational inference
references.authors: Tilmann Gneiting Adrian E Raftery
references.title: Strictly proper scoring rules, prediction, and estimation
references.venue: Journal of the American statistical Association,
references.authors: W. Keith Hastings
references.title: Monte Carlo sampling methods using Markov chains and their applications
references.authors: Matthew D. Hoffman Andrew Gelman
references.title: The No-U-Turn Sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo
references.venue: Journal of Machine Learning Research,
references.authors: Richard Jordan David Kinderlehrer Felix Otto
references.title: The variational formulation of the Fokker–Planck equation
references.venue: SIAM Journal on Mathematical Analysis,
references.authors: Takafumi Kanamori Masashi Sugiyama
references.title: Statistical analysis of distance estimators with density differences and density
references.venue: ratios. Entropy,
references.authors: Anna Korba Adil Salim Michael Arbel Giulia Luise Arthur Gretton
references.title: A non-asymptotic analysis for stein variational gradient descent
references.venue: Advances in Neural Information Processing Systems,
references.authors: Randall J LeVeque
references.title: Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-state and Time-dependent Problems, volume 98
references.authors: Chang Liu Jingwei Zhuo Pengyu Cheng Ruiyi Zhang Jun Zhu
references.title: Understanding and accelerating particle-based variational inference
references.venue: Proceedings of the 36th International Conference on Machine Learning,
references.authors: Chang Liu Jingwei Zhuo Jun Zhu
references.title: Understanding MCMC dynamics as flows on the Wasserstein space
references.venue: Proceedings of the 36th International Conference on Machine Learning,
references.authors: Qiang Liu
references.title: Stein variational gradient descent as gradient flow
references.venue: Advances in Neural Information Processing Systems,
references.authors: Qiang Liu Dilin Wang
references.title: Stein variational gradient descent: A general purpose bayesian inference algorithm
references.venue: Advances in Neural Information Processing Systems,
references.authors: Jianfeng Lu Yulong Lu James Nolen
references.title: Scaling limit of the Stein variational gradient descent: The mean field regime
references.venue: SIAM Journal on Mathematical Analysis,
references.authors: Nicholas Metropolis Arianna W Rosenbluth Marshall N Rosenbluth Augusta H Teller Edward Teller
references.title: Equation of state calculations by fast computing machines
references.venue: The journal of Chemical Physics,
references.authors: Radford M. Neal
references.title: MCMC Using Hamiltonian Dynamics, chapter 5
references.venue: CRC Press,
references.authors: Gareth O Roberts Osnat Stramer
references.title: Langevin diffusions and Metropolis-Hastings algorithms
references.venue: Methodology and computing in applied probability,
references.authors: Gareth O. Roberts Richard L. Tweedie
references.title: Exponential convergence of Langevin distributions and their discrete approximations
references.authors: Adil Salim Anna Korba Giulia Luise
references.title: The wasserstein proximal gradient algorithm
references.venue: arXiv preprint arXiv:2002.03035,
references.authors: Adil Salim Lukang Sun Peter Richtárik
references.title: Complexity analysis of stein variational gradient descent under talagrand’s inequality t1
references.venue: arXiv preprint arXiv:2106.03076,
references.authors: Luke Tierney
references.title: Markov Chains for exploring posterior distributions
references.venue: The Annals of Statistics,
references.authors: Cédric Villani
references.title: Optimal Transport: Old and New, volume 338
references.venue: Springer Science & Business Media,
references.authors: Martin J. Wainwright Michael I. Jordan
references.title: Graphical models, exponential families, and variational inference
references.venue: Foundations and Trends in Machine Learning,
references.authors: Max Welling Yee Whye Teh
references.title: Bayesian learning via stochastic gradient Langevin dynamics
references.venue: In Proceedings of the 28th international conference on machine learning,
references.authors: Michael Zhu Chang Liu Jun Zhu
references.title: Variance reduction and quasi-Newton for particle-based variational inference
references.venue: Proceedings of the 37th International Conference on Machine Learning,
references.authors: Jingwei Zhuo Chang Liu Jiaxin Shi Jun Zhu Ning Chen Bo Zhang
references.title: Message passing Stein variational gradient descent
references.venue: Proceedings of the 35th International Conference on Machine Learning,
sections.heading: 1 INTRODUCTION
sections.text: Sampling from unnormalized distributions plays a 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 a 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 on 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 et al., 1953; Hastings, 1970; Tierney, 1994; Dunson & Johndrow, 2019) provides a flexible framework for general MCMC sampling. To implement a Metropolis–Hastings algorithm, one needs to specify a proposal density and an acceptance policy. However, without a careful design of these two aspects, the Metropolis–Hastings algorithm can be inefficient due to strong correlations, slow mixing, or 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 (Welling & Teh, 2011; Chen et al., 2014).
  Variational Bayesian inference (Beal, 2003), often simply referred to 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 to trade off the parameterization flexibility of variational posteriors against the optimization complexity of ELBO in practice.
  In the spirit of VI, particle-based variational inference (ParVI) (Liu & Wang, 2016; Chen et al., 2018; Zhu et al., 2020) iteratively optimizes a set of particles to mimic a functional gradient descent for minimizing the KL divergence. ParVI seeks to move a variational distribution towards the unnormalized target distribution, along a steepest descent direction of the KL divergence. In a 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 of 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 a reproducing kernel Hilbert space (RKHS), which is further recognized as a gradient flow under the Stein geometry (Liu, 2017; Lu et al., 2019; Duncan et al., 2019). A drawback of SVGD is that it tends to collapse at part of the modes of the target, due to a negative correlation between the data dimensionality and the repulsive force in the RKHS (Zhuo et al., 2018).
  In this work, we propose a relative entropy gradient sampler (REGS) for sampling from unnormalized target distributions. To approximate a 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 of the relative entropy gradient flow. Therefore, to sample with REGS, we 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 on this observation, we propose a novel logarithmic density ratio estimation method for unnormalized distributions. By alternating between particle evolution and velocity field estimation, we can collect a set of stable particles which are approximately distributed as the target distribution. Our contributions can be summarized as follows:
  (1) Building upon the relative entropy gradient flow, we propose the relative entropy gradient sampler (REGS) for unnormalized target distributions. REGS preserves high efficiency and strong stability with respect to increasing singularity in mixtures of Gaussians, when the number of components increases (as shown in Figure 1), 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 a nonparametric approach to estimating the density ratio between an unnormalized density and an underlying density represented by samples, which is of independent interest.
  (4) We present experimental comparisons on varieties of multi-mode synthetic data and benchmark data and demonstrate that REGS is a 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 a 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
  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 Euler-Maruyama 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 a Metropolis-Hastings 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 an 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 P2(X ) be the space of Borel probability measures on a support space X ⊂ Rd with a finite second moment, and let Pa2 (X ) be a subspace of P2(X ) whose measures are absolutely continuous w.r.t. the Lebesgue measure. All probability measures we considered thereinafter are assumed to belong to Pa2 (X ). To ease the notation, we use probability density functions such as q(x), p(x),x ∈ X to express probability distributions in Pa2 (X ). Let (Pa2 (X ),W2) denote the metric space Pa2 (X ) 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.
sections.heading: 2 PROBLEM FORMULATION
sections.text: Consider an unnormalized probability density function u : X → [0,∞), where X ⊆ Rd is the support of u. Suppose u has an intractable normalizing constant Z = ∫ X u(x)dx <∞. Our goal is to generate random samples from the underlying distribution p ∈ Pa2 (X ), whose probability density function is only known up to proportionality, i.e., p(x) = u(x)/Z,x ∈ X . The basic idea is to gradually optimize samples from a given distribution q ∈ Pa2 (X ) to approximate samples from p, where it is easy to sample from q. Optimizing samples leads to functional optimization of distributions. We then introduce the classical relative entropy as the functional optimization objective. The relative entropy, a.k.a., the Kullback–Leibler divergence, for q, p ∈ Pa2 (X ) is the average logarithmic density ratio, which is defined as
  Dre(q‖p) = ∫ X q(x) log ( q(x) p(x) ) dx. (1)
  It holds that Dre(q‖p) ≥ 0 and Dre(q‖p) = 0 iff q(x) = p(x) a.e. x ∈ X . Moreover, we denote the relative entropy functional as
  F [·] := Dre(·‖p) : Pa2 (X )→ [0,∞]. (2) To sample from the unnormalized density u = pZ, we consider the functional minimization problem
  min q∈Pa2 (X )
  F [q], (3)
  where F [q] is always minimized at the underlying target distribution p, i.e., q(x) = p(x) a.e. x ∈ X . In a nutshell, problem (3) is an energy functional minimization problem in a metric space. To minimize the energy functional F , it suffices to move along the corresponding gradient flow in a metric space until the flow converges. For example, a 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 a curve that points in the steepest descent direction of a given energy functional. When equipped with the 2-Wasserstein distance, minimization of the energy functional F naturally corresponds to a 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.
sections.heading: 3 RELATIVE ENTROPY GRADIENT FLOW
sections.text: 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 (Pa2 (X ),W2). Recall that F in (2) is the relative entropy functional defined on (Pa2 (X ),W2). One can show that a curve {qt}t≥0 in (Pa2 (X ),W2) is a relative entropy gradient flow of F if it satisfies the continuity equation (Ambrosio et al. (2008), page 295 and Villani (2008), page 631),
  ∂tqt = Div ( qt∇
  δF [qt] δqt
  ) , (4)
  where qt(x) = q(t,x) evolves over time, δF [qt] δqt = log qtp is the first variation of the energy functional F at qt, and∇ δF [qt]δqt is the Euclidean gradient of δF [qt] δqt
  . Here, we identify the gradient as the relative entropy gradient, which is defined by
  ∇W2F [qt] := ∇ δF [qt] δqt = ∇ log qt p . (5)
  Moreover, the relative entropy F dissipates along the relative entropy gradient flow {qt}t≥0 at the rate (Ambrosio et al. (2008), page 295)
  ∂tF [qt] = −Eqt [‖∇W2F [qt]‖2]. (6) Therefore, the relative entropy gradient flow {qt}t≥0 eventually converges to the target distribution p as t→∞. As pointed out in Ambrosio et al. (2008) (Page 175), under mild conditions the continuity equation (4) concerning {qt}t≥0 determines a time-inhomogeneous Markov process {Xt}t≥0 that starts at a random particle X0 ∼ q0 and follows the particle evolution dynamics
  dXt dt = vt(Xt), Xt ∼ qt, t ≥ 0. (7)
  Note that the velocity fields
  vt = −∇W2F [qt] = ∇ log p
  qt , t ≥ 0 (8)
  drive the evolution of the particle Xt in the Euclidian space, which results in the transport of qt in (Pa2 (X ),W2). An important observation is that
  vt = ∇ log p qt = ∇ log u qt , t ≥ 0. (9)
  Therefore, the velocity fields do not involve the unknown normalizing constant Z. This is the key motivation for us to use the relative entropy gradient flow in the proposed method.
sections.heading: 4 SAMPLING AS PARTICLE EVOLUTION
sections.text: As indicated by the energy dissipation of relative entropy F in (6), running the relative entropy gradient flow {qt}t≥0 dynamics can provide a nice approximate solution to the functional minimization problem (3) when time t is large enough. Therefore, to sample from the target distribution p, it is appropriate to simulate the relative entropy gradient flow {qt}t∈[0,T ] 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 (LeVeque, 2007) as follows,
  Xk+1 = Xk + svk(Xk), X0 ∼ q0, k = 0, 1, . . . ,K − 1, (10) with the velocity field at step k
  vk = −∇W2F [qk] = ∇ log p
  qk , (11)
  where s > 0 is a tunable small step size, K = bT/sc is the number of iterations and qk is the corresponding discretized gradient flow at step k, i.e., Xk ∼ qk. Combining the expressions in (10) and (11), we have that the iterations progress according to
  Xk+1 = Xk + s∇ log p
  qk , X0 ∼ q0, k = 0, 1, . . . ,K − 1. (12)
  In principle, it is necessary to evaluate the velocity field vk = ∇ log(p/qk) each iteration in (12). By (9), the velocity field of the relative entropy gradient flow can be simplified to
  vk = ∇ log p qk = ∇ log u qk , k = 0, 1, . . . ,K − 1, (13)
  where u = Zp is the given unnormalized density of the target distribution p. Then only the density qk remains unknown for evaluating the velocity field. Ideally, qk can be estimated by evolving a large number of particles {Xik}Ni=1. However, direct estimation of qk is difficult due to the curse of dimensionality and the potential expensive computation cost for different ks. Our solution is to approximate the velocity field (13) as a whole.
  Assuming a nice approximation v̂k of the velocity field (13) is provided, then one can implement the following iterations for approximately sampling from qK with no effort,
  X̃k+1 = X̃k + sv̂k(X̃k), X̃0 ∼ q0, k = 0, 1, . . . ,K − 1. (14)
  Through the iterations above, we can collect X̃k ∼ q̃k ≈ qk, k = 1, 2, . . . ,K. We will discuss approximation of the velocity field vk = ∇ log(u/qk) from the perspective of estimating the logarithmic density ratio log(u/qk) in the next section.
sections.heading: 5 LOGARITHMIC DENSITY RATIO ESTIMATION AND THE RELATIVE ENTROPY GRADIENT SAMPLER
sections.text: In this section, we first propose a novel estimation procedure of the logarithmic density ratio log(u/q) based on an unnormalized density u and random samples from q.
  We use a model ratio R : X → [0,∞) to fit the true ratio R?uq = u/q between a density q and an unnormalized density u. Let g : R → R be a differentiable and strictly convex function. A Bregman score (Dawid, 2007; Gneiting & Raftery, 2007; Kanamori & Sugiyama, 2014) with the base probability measure q ∈ Pa2 (X ) to measure the discrepancy between R and R?uq is defined by
  B(R) = EX∼q[g′(R(X))R(X)− g(R(X))]− EX∼w [ u(X)
  w(X) g′(R(X))
  ] ,
  where w ∈ Pa2 (X ) 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 of w. Additionally, B(R) ≥ B(R?uq), where the equality holds iff R(x) = R?uq(x) (q, u)-a.e. x ∈ X .
  In this work, we take g(x) = x log(x)− x. We use this function for two reasons: (a) convexity, this 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 as
  B(R) = EX∼q[R(X)]− EX∼w [ u(X)
  w(X) log(R(X))
  ] . (15)
  Recall that the true density ratio R?uq can be factorized as R ? uq = u/q = Z(p/q). Thus the numerical scale of the true density ratio R?uq hinges on two factors, i.e., the normalizing constant Z of u and the standard density ratio p/q. 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 R as a byproduct.
  From now on, we denote D?uq = log(R ? uq), D = log(R) : X → R. Then B(D) can be rewritten as B(D) = EX∼q[exp(D(X))]− EX∼w [ u(X)
  w(X) D(X)
  ] . (16)
  It can be shown that the logarithmic density ratio D?uq is identifiable at the population level by minimizing (16) with respect to D. Theorem 1. For B(D) defined in (16), we haveD?uq ∈ arg minDB(D). In addition, for anyD with EX∼w [ u(X) w(X)D(X) ] <∞, B(D) ≥ B(D?uq), with equality iffD(x) = D?uq(x) (q, u)-a.e. x ∈ X .
  Algorithm 1: REGS: Relative entropy gradient sampler Input: u = Zp // unnormalized target density step size s > 0, an integer K > 0, // step size, maximum loop count X̃i0 ∼ q0, i = 1, 2, . . . , n // initial particles w ∈ Pa2 (X ) // reference distribution k ← 0 while k < K do
  Y ik ∼ w, i = 1, 2, . . . , n // reference samples D̂φk ∈ argminDφ 1 n ∑n i=1 [ exp(Dφ(X̃ i k))− u(Y ik )
  w(Y i k ) Dφ(Y
  i k ) ]
  // log density ratio
  v̂k(x) = ∇D̂φk (x) // velocity field X̃ik+1 = X̃ i k + sv̂k(X̃ i k), i = 1, 2, . . . , n // update particles
  k ← k + 1 end Output: X̃iK ∼ q̃K ≈ p, i = 1, 2, . . . , n // output particles
  Based on Theorem 1, we can estimate the unknown logarithmic density ratio D?uqk = log(u/qk) with a deep neural network Dφ with parameter φ through the sample version of (16). Let {X̃ik}ni=1 be i.i.d. samples from q̃k ≈ qk and {Y ik}ni=1 be i.i.d. samples from a reference distribution w. We solve the following deep nonparametric estimation problem via stochastic gradient descent (SGD) for D̂φk
  D̂φk ∈ arg min Dφ
  B̂(Dφ) = 1
  n n∑ i=1 [ exp(Dφ(X̃ i k))− u(Y ik ) w(Y ik ) Dφ(Y i k ) ] . (17)
  With the logarithmic density ratio estimator D̂φk , the velocity field vk in (13) can be approximately computed by v̂k = ∇D̂φk . By considering sampling as a particle evolution process discussed in Section 4, REGS updates the initial particles {X̃i0}ni=1 with iterations in (14) as follows:
  X̃ik+1 = X̃ i k + sv̂k(X̃ i k), X̃ i 0 ∼ q0, i = 1, 2, . . . , n, k = 0, 1, . . . ,K − 1. (18)
  We summarize the proposed REGS for sampling from an unnormalized density in Algorithm 1.
sections.heading: 6 NUMERICAL EXPERIMENTS
sections.text: We evaluate REGS on a large number of 1D and 2D mixture distributions and test its stability in the high-dimensional setting with multivariate Gaussian distributions. We also use REGS to perform Bayesian logistic regression on benchmark datasets. For comparison, we consider three existing methods including SVGD (Liu & Wang, 2016), ULA (Roberts & Tweedie, 1996) and MALA (Roberts & Tweedie, 1996). All experiments are done using a 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 https://github.com/anonymous/REGS.
sections.heading: 6.1 MIXTURE DISTRIBUTIONS
sections.text: We run REGS and SVGD, ULA and MALA to generate 2000 particles for mixtures of 2, 8 and 9 Gaussians (see Scenarios 4, 5, 6 in Appendix B), and 5000 particles for a mixture of 25 Gaussians (see Scenario 9 in Appendix B). The sampling qualities of these algorithms are compared by scatter plots with density contours of target mixture distributions. We classify all scatter points with labels according to the nearest mode, and plot the histograms of the label counts.
  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 is able to 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 2(c) and 2(d) show that MALA and ULA with 50 chains find all modes but with unequal weights, especially as the variance of each component decreases.
  Gaussian mixtures with unequal weights Figure 3 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 unequal weights (1, 1, 1, 1, 3, 3, 3, 3)/16. Figure 3(a) shows that the samples generated by REGS have the correct weights. Figures 3(c) and 3(d) indicate that ULA
  and MALA assign particles to 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 a test function h in Table 1, where h(x) = αTx, (αTx)2, and 10 cos(αTx+ 1/2) with α ∈ R2, ‖α‖2 = 1, and X is distributed as various Gaussian mixtures with unequal weights. By comparing the Monte Carlo estimates of E[h(X)] using the samplers with the values based on target samples, we see that REGS performs better and is more stable than SVGD, ULA and MALA, especially when h(x) = 10 cos(αTx+ 1/2). We include additional numerical results including more scatter plots and histograms (Figure 6-10) and Monte Carlo estimates with equal wieights (Table 6) in Appendix D.
sections.heading: 6.2 MULTIVARIATE GAUSSIAN DISTRIBUTION
sections.text: Let the target distribution be a d-dimensional Gaussian distribution with mean µ = (1, 1, · · · , 1) ∈ Rd and covariance matrix Σ ∈ Rd×d,Σi,j = ρ|i−j| with ρ = 0.7. We consider four test functions h(x), i.e., h(x) = αTx (the first moment), h(x) = (αTx)2 (the second moment), h(x) = exp(αTx) (the moment generating function), and h(x) = 10 cos(αTx + 1/2) with α ∈ Rd satisfying ‖α‖2 = 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 is 5000 for each sampler, where ULA_k and MALA_k denote the ULA and MALA with k 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 4 presents these Monte Carlo estimates as d increases from 10 to 300 with step size 10. The logarithm of the estimated E[exp(αTX)] is shown. As shown in Figure 4, the estimates using REGS and SVGD have smaller fluctuations than those using ULA and MALA, although all four methods can estimate E[αTX] and E[(αTX)2] well. Moreover, the third and the fourth panels in Figure 4 show that REGS outperforms SVGD, ULA and MALA when h(x) = exp(αTx) or 10 cos(αTx+ 1/2).
sections.heading: 6.3 BAYESIAN LOGISTIC REGRESSION
sections.text: 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 Liu & Wang (2016) and the first four datasets had been analyzed in Gershman et al. (2012). We consider a similar setting to that in (Liu & Wang, 2016; Gershman et al., 2012), which assigns a Gaussian prior π(β|α) = N (0, α−1I) to the regression coefficient β (including the intercept). We specify the prior of α as π(α) = 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 on test data with 5000 particles from the posterior. Table 2 lists the averages and standard errors (in parenhteses) of test accuracy. From Table 2 we can see that REGS is comparable with SVGD, ULA and MALA. For the Covertype dataset, MALA failed to converge, so no results from it are included in Table 2.
sections.heading: 6.4 DISCUSSION OF THE EXPERIMENTAL RESULTS
sections.text: 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 a multimodal distribution even if they are initialized in a different regions. An intuitive explanation is as follows. The movement of the REGS particles is determined by the velocity field. If the velocity field is 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 it 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 a principled way to move particles from an initial reference distribution to a multimodal distribution, albeit with a higher computational cost.
sections.heading: 7 CONCLUSION
sections.text: We have introduced REGS, a novel gradient flow based method for sampling from unnormalized distributions. Extensive numerical experiments demonstrate that REGS performs better than several existing popular sampling methods in 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.
sections.heading: A APPENDIX
sections.text: In this appendix, we prove Theorem 1, give a detailed description of the models in the numerical experiments, the neural network architecture used in implementing REGS, and additional numerical results.
sections.heading: B PROOF OF THEOREM 1
sections.text: Proof. By the definition of Bregman score B(R) between u and q, it is easy to check
  R?uq ∈ arg min R B(R),
  B(R) ≥ B(R?uq) with equality iff R(x) = R?uq(x) (q, u)-a.e. x ∈ X .
  Since D?uq = log(R ? uq), D = log(R),B(R) = B(D), when EX∼w [ u(X) w(X)D(X) ] <∞, we have
  D?uq ∈ arg min D B(D),
  B(D) ≥ B(D?uq) with equality iff D(x) = D?uq(x) (q, u)-a.e. x ∈ X .
sections.heading: C GAUSSIAN MIXTURE DISTRIBUTIONS
sections.text: Let the density of a d-dimensional multivariate Gaussian distribution with mean µ ∈ Rd and covariance matrix Σ ∈ Rd×d be
  f(x;µ,Σ) = 1
  (2π) d 2 √ det(Σ) exp
  { −1
  2 (x− µ)TΣ−1(x− µ)
  } .
  We consider several scenarios below for the unnormalized density function u(x).
  Scenario 1. 2Gaussians_1d1. One-dimensional mixture of 2 Gaussians,
  u(x) = 1
  3 f(x;µ1, σ
  2 1) +
  2 3 f(x;µ2, σ 2 2),
  where µ1 = 1, µ2 = −2, and σ21 = 0.25, σ22 = 2; Scenario 2. 2Gaussians_1d2. One-dimensional mixture of 2 Gaussians,
  u(x) = 1
  3 f(x;µ1, σ
  2 1) +
  2 3 f(x;µ2, σ 2 2),
  where µ1 = 3, µ2 = −3, and σ21 = 0.25, σ22 = 2; Scenario 3. 2Gaussians_1d3. One-dimensional mixture of 2 Gaussians,
  u(x) = 1
  3 f(x;µ1, σ
  2 1) +
  2 3 f(x;µ2, σ 2 2),
  where µ1 = 3, µ2 = −3, and σ21 = 0.03, σ22 = 0.03; Scenario 4. 2Gaussians. Two-dimensional mixture of 2 Gaussians,
  u(x) = f(x;µ1,Σ1) + f(x;µ2,Σ2),
  where µ1 = (r, 0)T, µ2 = (−r, 0)T, and Σ1 = Σ2 = σ2I. Scenario 5. 8Gaussians. Two-dimensional mixture of 8 Gaussians,
  u(x) = 8∑ j=1 f(x;µj ,Σj),
  where µj = r(sin(2(j − 1)π/8), cos(2(j − 1)π/8))T, and Σj = σ2I for j = 1, · · · , 8.
  Scenario 6. 9Gaussians. Two-dimensional mixture of 9 Gaussians,
  u(x) = 3∑ j=1 3∑ k=1 f(x;µjk,Σjk),
  where µjk = 4(j − 2, k − 2)T, Σjk = σ2I for j = 1, 2, 3 and k = 1, 2, 3. Scenario 7. 16Gaussians_1c. Two-dimensional mixture of 16 Gaussians,
  u(x) = 16∑ j=1 f(x;µj ,Σj),
  where µj = 4(sin(2(j − 1)π/16), cos(2(j − 1)π/16))T, and Σj = 0.03I for j = 1, · · · , 16.
  Scenario 8. 16Gaussians_2c. Two-dimensional mixture of 16 Gaussians,
  u(x) = 8∑ j=1 f(x;µj ,Σj) + 8∑ k=1 f(x;µk,Σk),
  where µj = 4(sin(2(j − 1)π/8), cos(2(j − 1)π/8))T, µk = 2(sin(2(j − 1)π/8), cos(2(j − 1)π/8))T, and Σj = Σk = 0.03I for j = 1, · · · , 8 and k = 1, · · · , 8.
  Scenario 9. 25Gaussians. Two-dimensional mixture of 25 Gaussians,
  u(x) = 5∑ j=1 5∑ k=1 f(x;µjk,Σjk),
  where µjk = 2(j − 3, k − 3)T, Σjk = σ2I for j = 1, · · · , 5 and k = 1, · · · , 5. Scenario 10. 49Gaussians. Two-dimensional mixture of 49 Gaussians,
  u(x) = 7∑ j=1 7∑ k=1 f(x;µjk,Σjk),
  where µjk = 32 (j − 4, k − 4) T, Σjk = 0.03I for j = 1, · · · , 7 and k = 1, · · · , 7.
  Scenario 11. 81Gaussians. Two-dimensional mixture of 81 Gaussians,
  u(x) = 9∑ j=1 9∑ k=1 f(x;µjk,Σjk),
  where µjk = 32 (j − 5, k − 5) T, Σjk = 0.03I for j = 1, · · · , 9 and k = 1, · · · , 9.
  Scenario 12. 1circle. Let µi = 4(cos(2iπ/N), sin(2iπ/N))T, i = 0, 1, · · · , N − 1 with N = 400. Consider three noise to be added to each point µi, including the uniform distribution U(µi, 1/30) on a disc with center at µi and radius 1/30, Gaussian distribution N(µi, 0.03I), and mixed these two distributions. Scenario 13. 2circles. Let µ1i = 2(cos(2iπ/N), sin(2iπ/N))T and µ2i = 4(cos(2iπ/N), sin(2iπ/N))T, i = 0, 1, · · · , N − 1 with N = 200. Consider three noise to be added to each point µki including uniform distribution U(µki, 1/30) on a disc with center at µki and radius 1/30, Gaussian distribution N(µki, 0.03I), and mixed these two distributions, k = 1, 2.
  Scenario 14. 1spiral. Let µi = 4iπN cos(4iπ/N), sin(4iπ/N)) T, i = 0, 1, · · · , N − 1 with N = 400.
  Consider three noise to be added to each point µi, including the uniform distribution U(µi, 1/30) on a disc with center at µi and radius 1/30, Gaussian distribution N(µi, 0.03I), and mixed these two distributions.
  Scenario 15. 2spirals. Let µ1i = 3iπN (cos(3iπ/N), sin(3iπ/N)) T and µ2i =
  − 3iπN (cos(3iπ/N), sin(3iπ/N)) T, i = 0, 1, · · · , N − 1 with N = 200. Consider three noise to be added to each point µki including uniform distribution U(µki, 1/30) on a disc with center at µki and radius 1/30, Gaussian distribution N(µki, 0.03I), and mixed these two distributions, k = 1, 2.
  Scenario 16. moons. Let µ1i = (8i/N − 6, 4 sin(iπ/N))T and µ2i = (8i/N − 2, 4 sin(iπ/N))T, i = 0, 1, · · · , N − 1 with N = 200. Consider three noise to be added to each point µki including uniform distribution U(µki, 1/30) on a disc with center at µki and radius 1/30, Gaussian distribution N(µki, 0.03I), and mixed these two distributions, k = 1, 2.
sections.heading: D EXPERIMENTAL SETTING
sections.text: 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(− 1h‖x − x
  ′‖22) and set the bandwidth as h = med2/ log n, where med is the median of pairwise distances between the particles {xi}ni=1. We set the learning rate of neural networks as 5e-4 for density ratio estimation in REGS. The network structures are presented in Table 4. The initial particles in REGS are sampled from Gaussian distributions.
  D.2 NEURAL NETWORK ARCHITECTURE
  Covertype: https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/ binary.html Banana, German, Image, Ringnorm: http://theoval.cmp.uea.ac.uk/ matlab/default.html
  E ADDITIONAL NUMERICAL RESULTS