Sampling with Mollified Interaction Energy Descent
Authors: Lingxiao Li
qiang liu
Anna Korba
Mikhail Yurochkin
Justin Solomon
Venue: ICLR 2023 poster (2023)

Abstract: Sampling from a target measure whose density is only known up to a normalization constant is a fundamental problem in computational statistics and machine learning. In this paper, we present a new optimization-based method for sampling called mollified interaction energy descent (MIED). MIED minimizes a new class of energies on probability measures called mollified interaction energies (MIEs). These energies rely on mollifier functions---smooth approximations of the Dirac delta originated from PDE theory. We show that as the mollifier approaches the Dirac delta, the MIE converges to the chi-square divergence with respect to the target measure and the gradient flow of the MIE agrees with that of the chi-square divergence. Optimizing this energy with proper discretization yields a practical first-order particle-based algorithm for sampling in both unconstrained and constrained domains. We show experimentally that for unconstrained sampling problems our algorithm performs on par with existing particle-based algorithms like SVGD, while for constrained sampling problems our method readily incorporates constrained optimization techniques to handle more flexible constraints with strong performance compared to alternatives.


Published as a conference paper at ICLR 2023
SAMPLING WITH MOLLIFIED INTERACTION ENERGY DESCENT
Lingxiao Li MIT CSAIL lingxi ac@mi
Qiang Liu
Anna Korba
University of Texas at Austin l ange utexas. edu
CREST, ENSAE, IP Paris anna koma.korb@eenae.
Mikhail Yurochkin
Justin Solomon
IBM Research, MIT-IBM Watson Al Lab mikhail, yurochkine ibm. com
MIT CSAIL
jsolomon@mit.. edu
ABSTRACT
Sampling from a target measure whose density is only known up a normal- ization constant is fundamental problem in computational statistics and ma- chine learning. In this paper, we present a new optimization-be based method for sampling called mollified interaction energy descent (MIED). MIED minimizes new class of energies on probability measures called mollified interaction ener- approximations gies (MIEs). These energies rely mollifier functions-smost the Dirac delta originated from PDE theory. We show that as the mollifier ap- proaches the Dirac delta, the MIE converges to the chi-square divergence with respect to the target measure and the minimizers of MIE converge to the target measure. Optimizing this energy with proper discretization yields practical first- order particle-based algorithm for sampling in both unconstrained and constrained domains. We show experimentally that for unconstrained sampling problems, our algorithm performs on par with existing particle- based algorithms like SVGD, while for constrained sampling problems our method readily incorporates con- strained optimization techniques to handle more flexible constraints with strong performance compared to alternatives.
INTRODUCTION
Sampling from unnormalized probability density is a ubiquitous task in statistics, mathemati- cal physics, and machine learning. While Markov chain Monte Carlo (MCMC) methods (Brooks et al., 20 provide way to obtain unbiased samples at the price of potentially long mixing variational inference (VI) methods (Blei et al., 2017) approximate the target measure with simpler (e.g., parametric) distributions lower computational cost. In this work, focus particular class of VI methods that approximate the target measure using a collection of interacting particles. A primary example Stein variational gradient descent (SVGD) proposed by Liu & Wang (2016), which iteratively applies deterministic updates to a set of particles to decrease the KL divergence to the target distribution.
times,
While MCMC and VI methods have found great success in sampling from unconstrained distri- butions, they often break down for distributions supported a constrained domain. Constrained sampling is needed when the target density is undefined outside a given domain (e.g.. the Dirich- let distribution), when the target density is not integrable the entire Fuclidean space (e.g., the uniform distribution), or when we only want samples that satisfy certain inequalities (e.g., fairness constraints in Bayesian inference (Liu et al., 2021)), 1)). A few recent approaches (Brubaker et al., 2012; Byme& & Girolami, 2013; Liu & Zhu, 2018: 8: Shi et al., 2021) extend classical sampling methods like Hamiltonian Monte Carlo (HMC) or SVGD to constrained domains. These extensions, however, typically contain expensive numerical subroutines like solving nonlinear systems of equations and require explicit formulas for quantities such as Riemannian metric tensors or mirror maps to be derived case-by-case basis from the constraints.

---

Published as a conference paper at ICLR 2023
We present optimization-busel method called mollified interaction energy descent (MIED) that minimizes mollified interaction energies (MIEs) for both unconstrained and constrained sampling. An MIE takes the form of double integral of the quotient of mollifier-soool approximation of 60. the Dirac delta at the origin-over the target density properly scaled. Intuitively, minimizing an MIE balances two types of forces: attractive forces that drive the current measure towards the target density, and repulsive forces from the mollifier that prevents collapsing. We show that as the mollifier converges to do. the MIE converges to the x² divergence to the target measure up to an additive constant (Theorem 3.3). Moreover, the MIE -ccnverges to x² divergence (Theorem 3.6), so that minimizers of MIEs converge to the target measure, providing a theoretical basis for sampling by minimizing MIE.
While mollifiers can be interpreted as kernels with diminishing bandwidths, our analysis is funda- mentally different from that of SVGI) where a fixed- bhndwdhad kernel is used to define reproducing kernel Hilbert space (RKHS) on which the Stein discrepancy has a closed-form (Gorham Mackey, 2017). Deriving version of the Stein discrepancy for constrained domains is far from trivial and requires special treatment (Shi et al., 2021; Xu, 2021). In contrast, our energy has unified form for constrained and unconstrained domains and approximates the x² divergence long the band- width is sufficiently small so that short-range interaction dominates the energy: this idea of using diminishing bandwidths in sampling is under-explored for methods like SVGD.
Al Algorithmically, we use first-order optimization to minimize MIEs discretized using particles. We introduce log-sum-exp trick neutralize the effect of arbitrary scaling of the mollifiers and the target density; this form also improves numerical stability significantly. Since we turn sampling into optimization, we can readily apply existing constrained sampling techniques such as reparameteri- zation using a differentiable (not necessarily bijective) map or the dynamic barrier method by Gong & Liu (2021) to handle generic differentiable inequality constraints. Our method is effective as only uses first-order derivatives of both the target density and the inequality constraint (or the repa- rameterization map), enabling large-scale applications in machine learning (see e.g. Figure 4). For unconstrained sampling problems, we show MIED achieves comparable performance to particle- based algorithms like SVGD, while for constrained sampling problems, MIED demonstrates strong performance compared to alternatives while being more flexible with constraint handling.
2 RELATED WORKS
KL gradient flow and its discretization for unconstrained sampling. The Wasserstein gradi- ent flow of the Kullback-Leibler (KI.) divergence has been extensively studied, and many popular sampling algorithms can be viewed as discretizations of the KL-divergence gradient flow. Two pri- mary examples are Langevin Monte Carlo (LMC) and Stein variational gradient descent (SVGD). LMC simulates Langevin diffusion and can be viewed
forward-flow splitting scheme for the KL-divergence gradient flow (Wibisono, 2018). At each iteration of LMC, particles are pulled along -V llog p where is the target density, while random Gaussian noise is injected, although a Metropo- lis adjusting step is typically needed for unbiased sampling. In contrast with LMC, SVGD is a deterministic algorithm that updates collection of particles using combination of an attractive force involving - V log and a repulsive force among the particles; it can be viewed as a kernelized gradient flow of the KL divergence (Liu, 2017) or of the x² divergence (Chewi al., 2020). The connection to the continuous gradient flow in the Wasserstein space is fruitful for deriving sharp convergence guarantees for these sampling algorithms (Durmus et al., 2019; Balasubramanian al., 2022: Korba et al., 2020; Salim et al., 2022).
Sampling in constrained domains. Sampling in constrained domains is more challenging com- pared to the unconstrained setting. Typical solutions are rejection sampling and reparameterization to an unconstrained domain. However, rejection sampling can have a high rejection rate when the constrained domain is small, while reparameterization maps need to be chosen on a case-by-case basis with a determinant-of-JJcooho propose constrained version of HMC for sampling on implicit submanifolds, but their algorithm is expensive as they need to solve a nonlinear system of equations for every integration step in each step of HMC. Byrne& & Girolami (2013) propose geodesic Hamiltonian Monte Carlo for sampling on embedded manifolds, but they require explicit geodesic formulae. Zhang et al. (2020); Ahn & Chewi (2021) propose discretizations of the Hirror-Langevin diffusion for constrained sampling
term that can be costly to evaluate. Brubaker et al. (2012)
2

---

Published as a conference paper at ICLR 2023
provided that mirror map is given to capture the constraint. Similarly, Shi et al. (2021) propose mirror SVGD that evolves particles using SVGD-like updates in the dual space defined by a mirror map. While these two methods have strong theoretical convergence guarantees, they are limited by the availability of mirror maps that capture the constraints. Liu et al. (2021) extend SVGD to incor- porate a population constraint over the samples obtained; their setting is different from ours where the constraint is applied to every sample.
Pairwise interaction energies on particles. Pairwise interaction energies on particles take the form [({=}]]) Σit W (1;; xy) for a kernel W (Ix y). The gradient flow of E on N particles can give rise to phenomena like flocking and swarming, and a long line of mathematical research studies mean-field convergence of such particle gradient flow continuous solutions of certain PDEs as N + 00 (Carrillo et al.. 2014; Serfaty, 2020). Of particular interest to sampling, a separate line of works summarized by Borodachov et al. (2019) demonstrates that for the hypersingular Riesz kernel W(x, y) ||2 yll2 with sufficiently big s, minimizing E over compact set yields uniform sam- ples as N - oc. Moreover, W can depend on p to obtain non-uniform samples distributed according to p. As the hypersingular Riesz kernel is not integrable, their analysis is entirely based on geometric measure theory and avoids any variational techniques. In contrast, we take a different approach by using integrable kernels in MIEs through mollifiers. This enables us to apply variational analysis to establish interesting connections between MIEs and x² divergence. We compare our formula- tion with theirs further in Appendix B. Joseph et al. (2019) propose to minimize the maximum of pairwise interaction akin to the interaction in Borodachov et al. (2019) without gradient information using a greedy algorithm. Recently, Korba et al. (2021) propose sampling via descending on kernel Stein discrepancy which can also be viewed as type of interaction energy, but like SVGD, it is limited to unconstrained sampling and can be slow as higher-order derivatives of log are required. Craig et al. (2022) propose approximating the x² divergence with a functional different from ours that also involves a mollifier. The resulting algorithm needs to evaluate an integral containing the target density over the whole domain each step which can be costly. Consequently, their method is only applied to special target densities for which the integral has closed form. In comparison, our method can handle general target densities with cheap updates in each step.
MOLLIFIED INTERACTION ENERGY
Notation. Let X C R" be the domain we work in. We use Ln to denote the Lebesgue mea- sure on R" and write J fz) dr to denote integration with respect to L. We will assume all measures are Borel. We use Ha to denote the d-dimensional Hausdorff measure. We use C R For x E RR" let 6z be the Dirac delta measure at and dwn be the empirical measure & Ek= Ez.- We denote LP(X) == { : X R Lebesgue measurable with fx [)(z)) dr < oo), and for s E LP(X) let I15lI,
wN to denote a set of points (sI,...,,IN)
(((x))) //8 We use || to denote the essential supremum of [f]. For g RD - R, we define convolution (z) := f))g((z z y) dy provided this integral exists. We use C* (X) to denote continuously k-differentiable functions and C (X) to indicate the smooth functions. The space of continuous k-differentiable functions with compact support on X is C² (X). We denote by P(X) the set of probability measures, and P2(X) the set of probability measures with bounded second moments. We denote by Maign( (X) the set of finite signed measures. Given a Lebesgue measurable map T X + X and μ E (x)).. Typ u is the pushforward measure of μ by T.
Setup. The domain X C R where constrain samples is assumed be closed and full- dimensional, i.e., (X) > 0; the unconstrained case corresponds to X = R The target density, always denoted as p, is assumed to satisfy p E (X),,, p(x) 0 for all E X, and 0 Ix p(x) dx < co. Let jı" be the target probability measure p² = p(x( dz/sx p(n) da for Borel B C X. Our goal is sample from μ²,
3.1 MOLLIFIERS
Our framework relies on the notion of mollifiers originally introduced in Friedrichs (1944).
Definition 3.1 (family of mollifiers). We say (0cleoo C C¹ (RR) is a family of mollifiers if it satisfies
3

---

Published as a conference paper at ICLR 2023
R", 6x(x) 2 0 and 0x(2) Oel (-(- (b) For any c > 0, || = and suPreR (x(2) < OC.
(a) For any
(c) For any > 0, E (1, lim IR" Bs 1e11, O.
As e - O, the distribution with density d. converges to do. and if f E L*(R") and is continuous, then convolution o. * converges to f both pointwise (Proposition A.4) and in L (Proposition A.7).
Our definition is different from the one used in PDE theory (Hörmander, 2015), where some á e C is fixed and the family of mollifiers generated by (<(x) : ="(xz/t). While in PDE theory mollifiers are typically used to improve the regularity of non-smooth functions, in our framework they are used construct interaction energies that approximate the x² divergence.
We will use the following mollifiers, In the following definitions, we include normalizing constants Ze that ensure || 1. We do not write out Z. explicitly as they might not admit clean forms and they are only relevant in theory but not in our practical algorithm. For >> n, the s-Riesz family of mollifiers is defined as (x) == (1|=l12 + < The Gaussiun family of molli- fiers is defined 6 (z) exp (-Il=ll//22)) /Zg. The Laplace family of mollifiers is defined as (x) exp (- l//e) /Z!. Since {dg) and { correspond the densities of centered Gaussian and Laplace random variables, they satisfy Definition 3.1. Proposition A. proves that { also satisfies Definition 3.1.
3.2 MOLLIFIED INTERACTION ENERGIES
Definition 3.2 (mollified interaction energy). Given a family of mollifiers {defeso satisfying Defi- nition 3.1, for each > 0, define symmetric kernel W. R x R + [0, 20 by
W. v) = { (x ((--)(((())p())-- v otherwise.
(1)
00
Define the mollified interaction energy (MIE), E. P(R") - [0, oo]. to be
(c(u) // We (z v) dp(x)) du(y).
(2)
Intuitively, minimizing (4) balances two forces: a repulsive force from d. y) makes sure μ has a good spread, and an attractive force from (p(=)p(y)) -1/2 helps f concentrate on high- density regions. The exponent -1/2 chosen to balance the two forces so that, we will show in Theorem 3.3, (p) approximates the x² divergence with respect to μ for small e.
3.3 CONVERGENCE TO 22 -DDVERRGENCE
Recall that the x2 diverggene between probability measures P, Q is defined by ( - 8
(2 PI
dp iQ P
(3)
otherwise,
where Q & P denotes absolute continuity Q with respect P with density dQ/ar.
Theorem 3.3. Suppose μ E P(R") satisfies x2(p 1²) < oc. Then, (x( < 0o for any > 0.
Furthermore,
lim (1( = (2 em>0
+1. +1.
We provide a succinct proof of Theorem 3.3 in Appendix A.2 using the theory of mollifiers devel- oped in Appendix A.1. In Remark A.10 discuss extension of Theorem 3.3 to cases when X not full-dimensional; in particular Theorem 3.3 still holds when X is "flat", i.., has Hausdorff dimension d and is contained in d-dimensional linear subspace.

---

Published as a conference paper at ICLR 2023
3.4 CONVEXITY AND T-CONVERGENCE
We next study properties of the minimizers of MIEs: Does mi))) & (n) admit a unique min- imum? If so, do minima of {EeJeso converge to μ² as e - 0? In order to answer affirmatively to these questions, we will need the associated kernel (z, == 0x y) to satisfy the following property.
Definition 3.4 (i.s.p.d. kernel), symmetric lower semicontinuous (Ls.c.) kernel K on X X is integrally strictly positive definite (i.s.p.d.) if for every finite signed Borel measure on X, the energy EK(v) == JJ (x. y) d(u x (2, y) well-defined (iee. the integrals over the negative and positive parts of are not both infinite,, and Ex(v) 2 0 where the equality holds only if v 0 on all Borel sets of X.
For the mollifiers we consider, the associated kernel (e, v) d. 4) is i.s.p.d on any compact set (Pronzato & Zhigljavsky (2021, Example 1.2)).
Proposition 3.5. Suppose is compact and ki (x, = Φ2 - y) is i.s.p.d. on X. Then,
(a) The kernel W. u) defined in (1) also i.s.pdd on X.
(b) The fanctional & is strictly convex on Magan (X) and atlains a unique minimum on P(X).
We next show the convergence of minima of {Ee} p² consequence of T-convergence of the sequence {Er).
Theorem 3.6. Suppose ke(s, v) == de v) is i.s.pdd. on compuct sets for every e > 0. Then we have I T-convergence (Definition A.13) E. I x² +1 as e + O. In particular, if X is compact, if we denote μ := arg min min (p). then Me + μ weakly and lim E. Ee())) =1
We prove Theorem 3.6 in Appendix A.3 using Fourier transforms and Bochner's theorem. Theo- rem 3.6 provides basis for minimizing & for small e, since its unique minimum will be a good approximation of y²
3.5 DIFFERENTIAL CALCULUS OF & IN P2(R")
We next study the gradient flow of & in Wasserstein space P2(R") for X ing the gradient flow of functional often provides insights into the convergence of algorithms that simulates gradient flow with time discretization (Ambrosio et al., 2005, Chapter 11) or spatial discretization (Chizat, 2022). Proofs of this section are given in Appendix A.4.
R". Understand-
Let F P2(RR") - R be a functional. The Wasserstein gradient flow of F (Ambrosio et al., 2005, Definition 11.1.1) is defined as the solution {Me (12 of the PDE:
0y.. V (,, where ot
) is Frechét subdifferential of F at μ. Intuitively, gradient flows capture the evolution of the variable being optimized if we were to do gradient descent' on F when the step sizes go to O. We next show that the gradient flow of &, agrees with that of x² || μ²) the sense that their subdifferentials coincide as + O.
WF E (,; R")
Proposition 3.7. Assume μ E P2(RR) has density q E (R"). Then any strong Frechét subdif- ferential (Definition A.20) of & at μ takes the form
(.)(x)) = 2V (p(x)-1/2(0 gt (D)(x)) for μ-a.e. x E R
(4)
Moreover, for sufficiently small e, de has compact support independent of e, then
lim We,p(x)) 2x2 (z), for p-a.e. R. --0
(5)
where W is strong Frechét subdifferential of x²(- (1)).
While simulating X2 divergence gradient flow often intractable (Trillos & Sanz-Alonso, 2020, Section 3.3), our MIE admits straightforward practical algorithm (Section 4).
'The more precise notion is of minimizing movements (Ambrosio et al., 2005, Chapter 2).
We assume that has compact support because W. can be unbounded and cause integrability issues due the presence of (p())p(v)) 1122 term.
5

---

Published as a conference paper at ICLR 2023
We next show that & is displacement convex at μt P2(RR) - O, obtaining result similar to Korba et al. (2021, Corollary 4). This hints that gradient flow initialized near p will have fast convergence.
C. (R"). Suppose that for sufficiently small e, o. has compact support independent of e. Assume ki (z. v) = dz(x y) i.s.p.d. Let & C. RR":R") and
Proposition 3.8. Suppose p
14 (I + #1) Then lime-0
d²
(pr) 2 0.
de2
t=0
4 A PRACTICAL SAMPLING ALGORITHM
We now present a first-order particl--base algorithm for constrained and unconstrained sampling by minimizing a discrete version of (2). Substituting an empirical distribution for μ in (2) gives the discrete mollified interaction energy, for N particles WN = (NN)..
Ea(wN)
d.(x; - (x)o=))))))))
(6)
N2
Denote wir, E arg min C CX Fe((wN) and u4, = arg pinn & (p)), If X compact, by Boro- dachov et al. (2019, Corollary 4.2.9), have weak convergence ówž,. + 1e as - oc. If in addition & (x, y) := (2 (2- y) is i.s.p.d. then by Theorem 3.6, we have weak convergence μr + p²" This shows that minimizing (6) with a large N and a small will result in empirical distribution of particles that approximates 1e", Our sampling method, mollified interaction energy descent (MIED), simply running gradient descent on (6) over particles (Algorithm 1) with update that resembles the one in SVGD (see the discussion in Appendix C.I). Below we address a few practical concerns.
Optimization in the logarithmic domain. In practical applications, we only have access to the unnormalized target density p. The normalizing constant of the mollifier & can also be hard to compute (e.g., for the Riesz family). While these normalization constants do not affect the minima of (6), they can still affect gradient step sizes during optimization. Moreover, ó can be very large when e is small, and in many Bayesian applications p can be tiny and only log is numerically significant. To address these issues, we optimize the logarithm of (6) using the log log-sum-ew, trick (Blanchard et al., 2021) to improve numerical stability and to get rid of the arbitrary scaling of the normalizing constants:
N
N
log (.(wN) = log Σ
Σ
exp log log - xz2) (log p(x;) (logp(x;) + logr(z))) - -2 log N.
(7)
Special treatment of the diagonal terms. Since de goes to the Dirac delta as c - O, the dis- cretization of (2) on a neighborhood of the diagonal {(x,,) extra care. In (6) the diagonal appears as >N o. (0)p(())) which can dominate the summation when is small and then o, (0) becomes too large. For the mollifiers that we consider, different diagonal term Σ Φe (h, ((z./) where h := min || xxyll and Kn 2 1 is constant depending only on the dimension n. Since 0 < (h,/kn) I (0). the energy obtained will be bounded between (6) and the version of (6) without the diagonal terms. Hence by the proof of Theorem 4.2.2 of Borodachov et al. (2019), the discrete minimizers of F. still converge to ne, as N + oo and e + O. Empirically we found the choice fin (1.3n) I/n works well for the Riesz family of mollifiers and we use the Riesz family primarily for our experiments.
= y) needs to be handled with
use
N
Constraint handling. If X * R", we need to minimize (7) subject to constraints WN (Z1...., ,Ny...,,N)].... X. Since the constraint is the same for each particle Fú, we want our algorithm to remain parallelizable across particles.
We consider two types of constraints: (a) there exists a differentiable map f R" + X, with C., (X f(R") = 0; (b) the set X is given by {i E R" g(x) 0} for a differentiable g R" + R. For (a), we reduce the problem to unconstrained optimization log (f(wn) with WN fz1...., 2N} C R" and f(wx) == (f(x))....)) f f(xN). For (b), apply the dynamic barrier method by Gong & Liu (2021) to particles in parallel; see Appendix C.2 for details and an extension to handling multiple constraints.
R using objective

---

Published as a conference paper at ICLR 2023
5 EXPERIMENTS
We compare MIED with recent alternatives on unconstrained and constrained sampling problems. Unless mentioned otherwise, we choose the s-Riesz family of mollifiers {Ør] with n+ 10- and = 10 we found minimizing the MIE with such mollifiers results in well-separated particles so that we can take с be very small as our theory recommends. This is not the case for the Gaussian or the Laplace family as setting c too small can cause numerical issues even when particles are well- separated. In Appendix D.5, we compare different choices of on constrained mixture distribution.
Unless mentioned otherwise: for SVGD (Liu & Wang, 2016), we use the adaptive Gaussian ker- as in the original implementation (adaptive kernels can be prone to collapse samples see Ap- pendix D.2); for KSDD (Korba et al., 2021). )) we use a fixed Gaussian kernel with unit variance. All methods by default use a learning rate 0.01 with Adam optimizer (Kingma Ba, 2014). The source code can be found at https: http://ggttt.tteeee/ com wwww///wwwiwool//////mo
5.1 UNCONSTRAINDE SAMPLING
Gaussians in varying dimensions. We first compare MIED with SVGD and KSDD on Gaussians of varying dimensions and with different numbers of particles. In Figure I, see that the number of dimensions increases, MIED results in best samples terms of W2 distance (Wasscrstein distance computed using linear programming) with respect to 10 i.i.d. samples, while SVGD yields lower energy distance (Szákely & Rizzo, 2013). We think this is because MIED results in more evenly spaced samples (Figure so the W₂ distance is lower, but it is biased since 0 for energy distance SVGD performs better. More details can be found in Appendix D.I.
Dimemniio 2
Dimensian2
Dimersion
Dimsmion
Dimamion B
Dimmmiion 8
SYGO
SVGD
SVCD R5D0 KIED
Erogg Disance
SVCD
SYGD
KSBb MED
KADD MIED
KSRD
5VED
MED
Nomber of parides
Namber of paridles
Numior of pertill
Namber putiils
Marker of patikks
Namber of patidks
Figure 1: Gaussian experiments results in varying dimensions and different numbers of particles. For each method, plot the metric (W2 distance energy distance) versus the number of particles, averaged over 10 trials (shaded region indicates standard deviation).
Product of two Student's :-distributions. Next, we consider a 2-dimensional distribution con- structed as the product of two independent t-distributions with the degree of freedom 2 com- posed with a linear transform. Unlike Gaussians, Student's t-distributions have heavy tails. On the left of Figure 2, visualize the results of each method with 1000 samples after 2000 iterations. MIED captures the heavy tail while SVGD fails. Quantitatively, on the right of Figure 2, while 3, MIED yields lower metries for bigger
SVGD captures the distribution in a.
[-a, a|2 better for
1
SVSD
Figure 2: Left: visualization of samples from each method for the product of two Student's t- distribution (composed with linear transform). Right: metrics of each method when restricting to |-a,a]² As t-distributions have heavy tails, if we draw i..dd samples from the true product distribution, a small number of them will have large norms, making the computation of metrics unstable. Thus we restrict both i.i.d. samples and the resulting samples from each method to [-a, a)² before computing the metrics for [2, 10).
7

---

Published as a conference paper at ICLR 2023
Bayesian logistic regression. We compare MIED with SVGD and KSDD for the Bayesian logistic regression setting in Liu & Wang (2016). We further include a naive baseline, independent particle descent (IPD), which simply runs gradient descent independently posterior probability. In addition to using test accuracy as the metric, we include W/2 distance and en- engy distance between the samples from each method and 104 samples from NUTS (Hoffman al., 2014) after sufficient burn-in steps. We summarize the results in Table 5.1. All methods, including the naive baseline IPD, are comparable in terms of test accuracy. In other words, accuracy is not good metric for comparing the quality of posterior samples. Bayesian inference is typically pre- ferred over maximum a posteriori estimation for its ability to capture uncertainty. When evaluating the quality of the uncertainty of the samples using distances between distributions, MIED provides the best approximation in terms of the W2 distance and all methods (except IPD) are comparable in terms of the energy distance.
each particle to maximize the
NUTS
IPD
SVGD
KSDD
banama (d= 3)
0.55
-6.14/3..880.55
7.811-5.24/0.55
-8.24/-.7600.55
-7.37/-5.0600555
breast cancer
10)
0.64
1.511-L0330.60
-1.62/22.0600.6
-1.71/-2.330.66
11.99/2.188/..
diahetis (d 9)
0.78
-2. 18/ 1,5500777
-3.09/-3 72
-2.91/-3.90/0.7
3/.1 13/0.77
flare solar (d 10)
0.59
3.30/2.6500 48
6.91/4.09/052
1.77/-0008/0.55
7.09/4.25/0.48
germun (d = 21)
0.65
-1 -1.80-1.22 1.25
-1.891-2.6590666 64
-1.27/-2.83/0.6
-1.96/-2.8000.66
heart 14)
0.87
D.40/ 0.56/0.87
-0.41/115500 5570..77
0.100-17660.887
-0,92/ 1.67//1667/0.88
image (d 19)
0.82
5.53//..1/0.83
7.17/4.01/0.83
2.16/-0.500088
1.14-18800 82
ringnorm (d = 22))
0.77
A.5N.77
4.111-5.9800.7
2.107-
4.03//5.7000.7
splice (d 61)
0.85
11 1880885
-1.221-2.6500.85
2.04/-005500.84
1.45/0,00/0
thyroid (d = G)
0.84
1 995/0 53/0.84
1.17/-0.0000.84
2.42/1.57/0.74
0.84/-0.37/0.84
titanic 4)
0.40
-0.46/-0.11/04400
-0.631-0.3900.40
-1.00%-0,45:0.40
twonorm (al = 21)
0.97
-1.21/-113.0.
-1.32/-2.78/0.97
1.55/-0,220.97
-1.44-3.21/00997
waveform (d 22)
0.77
-.677-18770077
-..98/5.2300.78
-2.60V 18800.77
33.9973.17/0.75
Bayesian logistic regression results with 1000 particles. We include the test accuracy for NUTS in the second column. Three numbers A/B/C in the following columns are logarithmic W2 distance, logarithmic energy distance, and test accuracy. Bold indicates the best numbers. We use 80% /20% training/test split. All methods Results are reported after 10a iterations.
Table 1:
with identical initialization and learning rate 0.01.
5.2 CONSTRAINED SAMPLING
Uniform sampling in 2D. We consider uniform sampling in the square [-1, 1|2 with 500 particles. We reparameterize our particles using tanh to eliminate the constraint and show results with vari- ous choices of mollifiers we always choose the smallest while the optimization remains stable. We compare our method with mirror LMC (Zhang et al., 2020) and SVMD/MSVGD by Shi et al. (2021) using entropic mirror map 6(0) Ei=1 ((1 0.) log(1 + 0,) + (1 0,) log(1 0.)). To demonstrate that SVGD and KSDD break down in constrained domains, we implement these two methods adapted to the constrained setting using the same reparameterization as our method. The initial particles are drawn uniformly from (-0.5, 0.5]2. The left plot of Figure 3 shows that quantita-
CAMD
Teall
Figure 3: Convergence of metrics and visualization of samples for uniform sampling from a 2D box.
tively our method achieves much lower energy distance and W22 distance (compared to a 5000 i.i.d. samples uniformly drawn in 1-1, 1)..). On the right of Figure 3, we visualize the samples from each method. We see that samples from SVGD with an adaptive kernel collapse- we investigate this issue in Appendix D.2. Samples from SVGD with a fixed kernel size and KSDD are non-uniform (we choose kernel sizes that produce the best result). While SVMD creates locally clustered arti- facts, MSVGD produces good results. For mirror LMC, the resulting samples are not evenly spaced,
8

---

Published as a conference paper at ICLR 2023
resulting in worse W2 distances. When using tanh the reparameteriiztiion map, MIED produces decent results for all three families of mollifiers (Riesz, Gaussian, Laplace). When using the same entropic mirror map, MIED also produces good results. This highlights the flexibility of our method with constraint handling, whereas, for LMC and SVMD/MSVGIL the mirror map has to be chosen carefully to be the gradient of strongly convex function while capturing the constraints.
In Appendix D.3, we further test SVMD/MSVGD with a different choice of the mirror map where they break down. In Appendix D.4, we conduct similar comparison for sampling from a 20- dimensional Dirichlet distribution using the same setup as Shi et al. (2021). In scenarios where a good choice of a mirror map is available, SVMD/MSVGD can obtain better performance compared to MIED. In Appendix D.6, we conduct additional qualitative experiments for MIED, demonstrating its effectiveness for challenging constraints and multi-modal distributions.
Fairness Bayesian neural networks. We train fair Bayesian neural networks to predict whether the annual income of a person is least $50, 000 with gender as the protected attribute using the Adult Income dataset (Kohavi et al., 1996). We follow the same setup as in Liu et al. (2021) where the dataset D (x((), y 2 2 consists of feature vectors (.)), labels y() (whether the in- come 2 $50, 000), and genders 20)) (protected attribute). The target density taken to be the posterior of logistic regression with two-layer Bayesian neural network i(:; 0) with weights 0, and we put a standard Gaussian prior on each entry of independently. Given ( >> 0, the fair- ness constraint is g0) (x.)(x...)))- [2, [2;; 0)]) S O. On the left of Figure 4, we plot the trade-off curve of the result obtained using our method and the methods from Liu et al. (2021) for é E {10-5 5, 10 0.0001, 0.001, 0.002, 0.005, 0.01}. Details can be found in Appendix D.7. Our method recovers a much larger Pareto front compared to the alternatives. On the right of Figure 4, visualize the curves of the energy and the covariance versus the number of training iterations: as expected we see a smaller / results in bigger MIE (lower log-likelihood) and smaller covariance between the prediction and the protected attribute.
0.16
MIED (Ours)
0.50
Langevi
Demagaphie Parter
0.22
Primal Dual+SVGD
Prinal-Dusl--aaa Hargeinn
0.30
6.08-
0.04
0.02
0.34
0.53
Accuracy
Figure 4: Left: trade-off curves of demographic parity versus accuracy on the test data for MIED and methods from Liu al. (2021). Right: MIEs and (ng,(),) [z, ()); (measured on the training data) versus the number of training iterations, for various t.
6 CONCLUSION
We present new sampling method by minimizing MIES discretized as particles for unconstrained and constrained sampling. This is motivated by the insight that MIES converge to x" divergence with respect to the target measure as the mollifier parameter goes to O. The proposed method achieves promising results on the sampling problems consider.
Below we highlight three limitations. First, as discussed in Remark A.10, our theory only applics when the domain is full-dimensional flat. Extending theory to handle cases where the domain is an arbitrary d-rectifiable set is an important next step as it allows the handling of more complicated constraints such as nonlinear equality constraints. Secondly, when > 0. the minimizer of MIE can be different from the target measure. Finding a way to debias MIE (e.g.. like how Sinkhorn distances debiased (Feydy
al., 2019)) is an interesting direction. Lastly. the connection between the minimizers of the discretized MIE (6) and those of the continuous MIE (2) is only established in the limit N - Do. We hope to investigate how well the empirical distribution of particles minimizing (6) approximates the target measure when N is finite as in Xu et al. (2022).

---

Published as a conference paper at ICLR 2023
Reproducibility statement. We provide self-contained proofs in Appendix A for the theoretical results stated in the main text. The source code for all experiments can be found at attpss// github. cwm http://c//nnggaoo11////1/////
Acknowledgements Qiang Liu would like to acknowledge the support of NSF CAREER 1846421 and ONR. The MIT Geometric Data Processing group acknowledges the generous support of Army Research Office grants W91 INF2010 16 and I91IN 10293, of Air Force Office of Scientific Research award 199550-19-1-05 of National Science Foundation grants IIS-1838071 and CHS- 1955697, from the CSAIL Systems that Learn program, from the MIT-IBM Watson AI Laboratory, from the Toyota-CSAIL Joint Research Center, from gift from Adobe Systems, and from Google Research Scholar award.
REFERENCES
Kwangjun Ahn and Sinho Chewi. Efficient constrained sampling via the mirror-langevin algorithm. Advances in Neural Information Processing Systems, 34:28405-28418, 2021.
Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient flows: in metric spaces and in the spuce of probability measures. Springer Science & Business Media, 2005.
Krishna Balasubramanian, Sinho Chewi, Murat A Erdogdu, Adil Salim, and Shunshi Zhang, To- wards theory of non-log-concar sampling: first-order stationarity guarantees for langevin monte carlo. In Conference on Learning Theory, pp. 1896-2923. PMLR, 2022.
Pierre Blanchard, Desmond Higham, and Nicholas J Higham. Accurately computing the log-sum- exp and softmax functions. IMA Journal of Numerical Analysis, 41 41 2021.
David Blei, Alp Kucukelbir, and Jon D McAuliffe, Variational inference: review for statisti- cians. Journal of the American statistical Association, 112(518):859 877, 2017.
Sergiy V Borodachov, Douglas P Hardin, and Edward B Saff. Discrete energy on rectifiable sets. Springer, 2019.
Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng. Handbbook of markov chain monte carlo. CRC press, 2011.
Marcus Brubaker, Mathicu Salzmann, and Raquel Urtasun. A family of meme methods on implicitly defined manifolds. In Artificial intelligence and statistics, pp. 161-172. PMLR, 2012.
Simon Bymee and Mark Girolami. Geodesic monte carlo on embedded manifolds. Scandinavian Journal of Statisticy, 40(4):825-845 2013.
José Antonio Carrillo, Young-Pil Choi, and Maxime Hauray. The derivation of swarming models: mean-field limit and wasserstein distances. In Collective dynamics from bacteria to crowds, pp. 1-46. Springer, 2014.
Sinho Chewi, Thibaut Le Gouic, Chen Lu, Tyler Maunu, and Philippe Rigollet. Svgd kemel- ized wasserstein gradient flow of the chi-squared divergence. Advances in Neural Information Processing Systems, 33:2098-21000 2020.
Lenaic Chizat. Sparse optimization measures with over-parameterized gradient descent. Mathe- matical Programming. 194(1):487-532, 2022.
Katy Craig, Karthik Elamvazhuthi, Matt Haberland, and Olga Turanova. blob method method for inhomogeneous diffusion with applications to multi-agent control and sampling. arXiv preprint arXiv:2202.129927, 2022.
Alain Durmus, Szymon Majewski, and Blažej Miasojedow. Analysis of langevin monte carlo via convex optimization. The Journal of Machine Learning Research, 20(1):26662271,, 2019.
Jean Feydy. Thibault Séjourné, François- Xavier
-anvir Vialard, Shun-ichi Amari, Alain Trouvé, and Gabriel Peyré. Interpolating between optimal transport and mmd using sinkhorn divergences. In The 22nd International Conference on Artificial Intelligence and Statisties, pp. 2681-2690. PMLR, 2019.
10

---

Published as a conference paper at ICLR 2023
Kurt Otto Friedrichs. The identity of weak and strong extensions of differential operators. Transac- tions of the American Mathematical Sociery, 55(1):133-155) 1944.
Chengyue Gong and Xingchao Liu. Bi-objective trade-off with dynamic barrier gradient descent. NeurlPS 2021, 2021.
Jackson Gorham and Lester Mackey. Measuring sample quality with kernels. In International Conference Machine Learning, pp- 2292-1301. PMLR, 2017.
Matthew Hoffman, Andrew Gelman, al. The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo. J. Mach. Learn. Res., 15(1):1593-162) 2014.
Lars Hörmander. The analysis of linear partial differential operators I: Distribution theory and Fourier analysis. Springer, 2015.
V Roshan Joseph, Dianpeng Wang, Li Gu, Shiji Lyu, and Rui Tuo. Deterministic sampling of expensive posteriors using minimum energy designs. Technometrics, 2019.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint urXiv:1412.6988. 2014.
Achim Klenke. Probability theory: a comprehensive course. Springer Science & Business Media, 2013.
Ron Kohavi al. Scaling up the accuracy of naive-bayes classifiers: decision-tree hybrid. In Kdd, volume 96, pp. 202-207, 1996.
Anna Korba, Adil Salim, Michael Arbel, Giulia Luise, and Arthur Gretton. A non-asymptotic anal- ysis for stein variational gradient descent. Advances in Neural Information Processing Systems, 33:4672 4682, 2020.
Anna Korba, Pierre-Cyril Frankosra,,, Szymon Majewski, and Pierre Ablin. Kernel stein discrepancy descent. In International Conference on Machine Learning, pp. 5719-5730. PMLR, 2021.
Franziska Kühn. Probability and heat kernel estimates for lévy (-type) processes. 2016.
Chang Liu and Jun Zhu. Riemannian stein variational gradient descent for bayesian inference. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018
Qiang Liu. Stein variational gradient descent as gradient flow. Advances in neural information processing systems, 30, 2017.
Qiang Liu and Dilin Wang. Stein variational gradient descent: general purpose bayesian inference algorithm. Advances in neural information processing systems, 29, 2016.
Xingchao Liu, Xin Tong, and Qiang Liu. Sampling with trusthworthy constraints: A variational gra- dient framework. Advances in Neural Information Processing Systems, 34:23557-23568 2021.
Lue Pronzato and Anatoly Zhigljavsky. Minimum-energy measures for singular kernels. Journal of Computational and Applied Mathematics, 382: 13089, 2021.
Adil Salim, Lukang Sun, and Peter Richtarik. A convergence theory for svgd in the population limit under talagrand's inequality tl. In International Conference on Machine Learning, pp. 19139- 19152. PMLR, 2022.
Sylvia Serfaty. Mean field limit for coulomb-type flows. Duke Mathematical Journal, 169(15): 2887-2935, 2020.
Jiaxin Shi, Chang Liu, and Lester Mackey. Sampling with mirrored stein operators. arXiv preprint arXiv:2106. /2506. 2021.
Gábor J Szěkely and Maria L Rizzo. Energy statistics: A class of statistics based on distances. Journal of statistical planning and inference, 143(8): 1249-1222 2013.

---

Published as a conference paper at ICLR 2023
Ryan Tibshirani. Dykstra's algorithm, admm, and coordinate descent: Connections, insights, and extensions. Advances in Neural Information Processing Systems, 30, 2017.
Nicolas Garcia Trillos and Daniel Sanz-Alonso. 'he bayesian update: variational formulations and gradient flows. Bayesian Analysis, 15(1):29-56, 2020.
SR Srinivasa Varadhan. Probability theory. Number 7. American Mathematical Soc., 2001.
Andre Wibisono. Sampling as optimization in the space of measures: The langevin dynamics as composite optimization problem. In Conference on Learning Theory, pp- 2093-3027. PMLR, 2018.
Lantian Xu, Anna Korba, and Dejan Slepcev. Accurate quantization measures via interact- ing particle- sasdell-bbasdd optimization. In International Conference Machine Learning, pp. 24576- 24595. PMLR, 2022.
Wenkai Xu. Generalised kernel stein discrepancy (gksd): A unifying approach for non-parametric goodness-of-fit testing. arXiv preprint arXiv:2106.1 12105, 2021.
Kelvin Shuangjian Zhang, Gabriel Peyré, Jalal Fadili, and Marcelo Pereyra. Wasserstein control of mirror langevin monte carlo. In Conference on Learning Theory, pp. 3814-3841, PMLR, 2020.
A DETAILED ANALYSIS
A.1 PRELIMINARIES ON MOLLIFIERS
Proposition A.1. For 8 the s-Riesz family of mollifiers, defined as (((x) = (1<=112 +
2) s/2 ZZ satisfies Definition 3.1.
In order to prove Proposition A.1, we first prove a lemma. Lemma A.2. Let 2 Then for any e > 0,
gn
z: lvilll dy = H Hn-1(5m-1)
C mtths
(8)
(t2 + 1)0/2
B2(0)
where 15-1(5"-1) is the volume of the (n 1) dimmenssio sphere. Furthermore, assuming > b b n, then for any > 0,8 0,
gn+b-=
z
1w01) dy Hnn (Sn-1) (n + b)
(B.
Proof. For (8), compute
|1
z: B.(u)
dw S0.00 (1v122 2)s/2
= Hn 1
(x² + 22)002 dr
= Hn-1 Hn-1(Sm-1) S'' (el)m+b-1
(et + 2))7 cdl
g-4b-1 (t2 1)0/2 dl uth
HH.-(Sm-11
where we use substitution el. on the second line.
If 8 > b then
/
1
z
R BB(0)
RR
(1lv) 2 + 22)s/2
||vl
/ 0
H., (8
||
Rd (.,BB(
8n+b-8
Hn 1
(n + b)
12

---

Published as a conference paper at ICLR 2023
where the last integral uses b
Proof of Proposition A.I. It is clear that (a) of Proposition A.1 holds for d2(x). Since (2 is bounded we have 0f E L 0 (R") For any ò > 0, we have, for S 8,
(x) Z2 S o2(x) de = Ce Bg(0)
[_
where we use (8) with b 0 and C' is constant depending only on n, On the other hand, using (9) with =0, 0,
%
(x) S C'8m=
(x(BB(())
where C depends only on 8. With 5 = e, we see that z E L (R") so (b) is satisfied.
Since 8 is fixed and > n, we see that
BB (f) dr
Ce
lim
lim
(.))B (2(x) dx
C"8n-s
Since Sn,(o) (1(x) S (1()) (z dr 1, we have shown (c).
In the rest of this section, assume (de is a family of mollifiers satisfying Definition 3.1. Lemma A.3. For any p E [1, oo), bo), for any 0, de (P(RR) and lim IR (B 0. holds.
Proof. Assume & {1, ) since both cases are covered in the assumptions. Then for any > 0,
by Hölder's inequality,
1 Ip = -0--o s 110e11ll@8~1111 = =10ee11llée11?" < oo.
Similarly,
1 10)0010 = =||_R= =100, 1R) 0f
< 11 [R= 10
Il 1 p-1
Letting + and applying (c) gives lim B. =0
Proposition A.4. Let f E LP(R"), pE E oo). Then every e > 0, the integral
I Mi - (y)
exists, so that (f (x)())) finit. Moreover, f continous at x, then
lim (z)(()) fz)).
(10)
Proof. Observe that by Hölder's inequality, for such that 1/p + 1/g 1 (allowing infinity),
/ |Is(e - (v) S II <c.
Hence (f * (x) is integrable and finite.
For any 0, note that
|(5 (x) f(x)) / v)0 dy--(uydy-(a0 1(0)
13

---

Published as a conference paper at ICLR 2023
Since dc())) 1, we have
(()+))(=) f(x)) - (f(2- 1 -(()=))))( dy - f |f(c=
/ y)
f(x))
f())lo
Fix > 0. Continuity of f at x implies there exists > such that If(z - y) - -f(x) < t for all
y E Bs()). Then
J BB((0)
s [B,(0) dy) St.
On the other hand, by Hölder's inequality,
J (R f(x - (()))(=)) dy &
- s 211fllp) f s
Hence
|(f * (x)(x))) f(x) st 2 1
By Lemma A.3, since ||/ co, taking - 0 we get, for any >
lim sup (o)())))))) f(x) st.
c-0
Now let t + we get (10).
Corollary A.S. Let I IP(R"), p E oo), be continuous function. If {Se} a sequence of functions that converge uniformly to J, then for every J,
lim * (0]xx)((a))) -00
(11)
Proof. Note that
l( * (z) - f(x)] s 1(6. * fz)()) - (6. * (z)) 1(6 f)(z) f(x)).
Proposition A.4 shows the second term goes
0. For the first term, have
(( fc)(2) - (o. * M(x) = /ʃ ( 3) - f(x 5)) (x dy
S sup (f(())) - fx)) +
by uniform convergence.
Lemma A.6. For f E LP(R"), p E (1), we have
lim w
Tyf flp 0,
where we use raf denote the translated function Ta Taf(x) = f(x - a).
Proof. Fix O. It is standard fact that CRR)") is dense in (P(R)). Hence there exists E CR") such that | - ellp < e. Since is continuous with compact support, it is uniform continuous. Then there exists 8 0 with 1g(2 y) g(x)) < </ (K) if Ba(x). Hence for such we have [ - =ll"" < with 2, < co. Thus
|| fll, ||r Tuollp + || allp + Ilg f11, < 3e.
Proposition A.7. Let f LPR"), p E 01) Then for every e > 0,
|| dell, < 3/1511
In particular, f o. E (R").. Moreover.
lim || 0. - -111p c00
14

---

Published as a conference paper at ICLR 2023
Proof. For every O, we have
[ + 0r - -1p / (f(x y) (())) (v) dy dx
/ 10(x ( (()))"(((0) dy
/ (/ (115co-)) v) ((x)) dz (c(1) dy
= / |Tb fflf-blliseee f1) dy.
where we use Jensen's inequality with the observation that de (y) dy is a probability measure and Tonelli's theorem to exchange the order of integration. To show the first claim, note that,
Ilf Oc - f11l - / 1 f1))0(()) dy - / (211fl dy (211 fll).
Hence |lf del1, 1l + Ilf o. fll, s 3ll
For the second claim, by Lemma A.6, the function y + Tly fllp is continuous at y 0 with limit Hence by Proposition A.4, we are done by taking e -
A.2 CONVERGENCE TO x2 -DIVERGENCE
We start by giving an alternative formula for x2 divergence with respect to the target measure.
Lemma A.8. Define a functional E :P(R") + [0, oc] by
Jx g())² p(a) D0
dr if dy(x)) - q(x) dLn and p n¹
E(H) =
(12)
otherwise.
Then for every μ E P(R"),
E(H) =
+1.
Proof. If μ. is not absolutely continuous with respect to j²", then both sides are infinity. Otherwise,
((X)= = 1 and z has some density q. We then compute
/ ftt (x) - (x)
2 du (x) (x))+==//
2 [((x) -1) p(x) dx +
x2 (1) +1=
(((x)
= / g(x)²
dx - 2 2/ ((x) do + = / q(x)2
p(x) dx = E(H).
p(x)
3.3. Since 2² || μ²) < co, we know μ < μ², so x(X) We let (()) be the density of which satisfies q(z) = 0 for x E R X. We compute using Tonelli's 's theorem (since our integrand positive):
Proof of Theorem
// ((le-v(((()ppuu)) 112 ((())())) dr
Ex(H)
=
/ (e - (z) (c) / (g. g) dr,
(13)
where we define g(x) = /<x))/p((z) /p(z) X and g(x) 0 for & X. Moreover, by Lemma A.8, x² 1²) +1 1= / g(x)2 dr.
15

---

Published as a conference paper at ICLR 2023
The assumption 2 || μ²*) < oo then implies g E L2(R)). By Proposition A.7, we have O. * g L (R()). Next notice that
p*) - 1)1 f l( * ())()) g(x)l1(2)) dz S ll * 0|lاll91l2 < <o.
(1)
This shows the first claim. Finally, the last expression goes to 0 as + 0 since | 9 gl12 - by Proposition A.7.
Remark A.9. One may ask if similar results as Theorem 3.3 f-divergences other than x² diver- gence can be obtained using similar techniques. We think the proof of Theorem 3.3 is highly specific x2 divergence because (13) there are two copies of q coming from the definition of & as double integral over p. An f-divergence between il and J has the form fx f(g(z)//(())) dp(x), and the only way to have g² showing up is to choose [zz =x² corresponding to x2 divergence.
Remark A.10. It possible to prove versions of Thhorrem 3.3 when X has Hausdorff dimension d < n. in such cases <ivdieergence still makes sense as does (2). When X is "flat", i.e. with Hausdorff dimension d and contained in a d-dimensional linear subspace of R", e.g., when X is defined by a set of linear inequalities, then Theorem 3.3 follows if we adapt the assumptions of Definition 3.1 and calculation in the proof of Theorem 3.3 to be in the subspace. For a general d-dimensional X, (3)) (x) dHL(z)). For similar argument to go through, we will need the normalizing constant fx Óe - y) dHL(()) to the same for all x. This is true when X is a d-dimensional sphere, but for a general X the last integral will depend on the base point x. Proving version of Theorem 3.3 for X with Hausdorff dimension < n is an interesting future direction.
similar calculation yields E((N) = fx (Jx (x v) (3p) (y) dHH(V)
A.3 CONVEXITY AND E-CONVERGENCE
We start by recalling few definitions regarding functionals in P(RR).
Definition A.11. We say a functional F P(R") + (-o0, oc] is proper if there exists μ E P(R") such that F(1) < oc, and is lower semicontinuous (L.sc.) if for any weakly convergence sequence l4k Js have lim inf F(p)) 2 F(P).
Lemma A.12. For any > O, the functional& E, P(R") + (-oc, oo] is proper and l.s.c. Moreover, if X is compact, the minimum min (x) is attained by some measure in P(X).
Proof. Taking any z E X, since &e is bounded and p(x) > O, we see that E, (82) < 0o so & is proper. Moreover, given weakly convergence P4 μ4, by Portmanteau theorem and the fact that W. is nonnegative and I.s.c., we conclude that E, is also l.s.c.
The set of probability distributions P(X) C P(R") is tight by the compactness of X. It is closed since if { C. P(X) weakly converges to μ. then by Portmanteau theorem, μ(X) Z lim sup (X) 1 so that y P(X). Hence by Prokhorov's 's theorem (Theorem 5.1.3 (Ambrosio et al., 2005)), P(X) is (sequentially) compact. It is then an elementary result that any l.s.c. function attains its minimum on compact set.
We next prove Proposition 3.5 regarding the convexity of & and the uniqueness of its minimum.
Proof of Proposition 3.5. Let μ be any finite signed Borel measure. The compactness assump- tion and the fact that p E (X) imply p(x) > 5 for any x E X for some > 0. Hence p(x) 1/2 < 8 8-1/2, so the weighted measure ğ defined by dā(x)) = p p-p--1 (111((2)) du(x) is also fi- have Ew, (n) &c, (ū) > 0 if f is not the zero
nite. By the definition of i.s.p.d. kernels, measure, which is equivalent to f not being zero since p > O. Thus W. is i.s.p.d. on X. Also note that (p(x)p(y))~1) 1/2 < 5-1 for all I, E X, so Ew. is always finite. By Lemma 1.1 of Pronzato & Zhigljavsky (2021), conclude that Ew, is strictly convex in Msgaa the space of finite signed measures, and in particular is convex on P(X). Hence combined with the existence result from Lemma A.12 we conclude & attains unique minimum in P(X).
16

---

Published as a conference paper at ICLR 2023
Our goal in the rest of this section is to prove the r-convergence of & to & P(R") 00 ool defined in (12). By Lemma A.8, we have E()) = x2 (1 ll p²) + 1. Hence together we will have proved Theorem 3.6.
Definition A.13 (T-convergence). A sequence of functionals F. P(R") (-o0, oo] said to I r-converge to F P(RR) (-o0, oo). denoted as F. + F. if:
(a) For any sequence pe P(R") converging weakly to μ E P(R"), lim inf Fe(Mee 2 F(u);
(b) For any μ E )( there exists a sequence He E P(R") converging weakly to p with lim limsup, Fe(ee) S F(u).
will show & & using Fourier transforms and Bochner's theorem.
Definition A.14 (Fourier transform). For E (R.), its Fourier transform j is the complex- valued function defined via
f(5) = / f())
More generally, for a signed finite measure Msignn (R"), its Fourier transform fi is the complex- valued function function defined via
(L{) := / - du(x))
This integral is always well-defined and moreover ja is uniformly continuous; see Borodachov et al. (2019, Section 1.10).
We will prove the following weak version (under the additional assumption that mollifier ø is inte- grable) of Bochner's s theorem suitable for our case. In particular we will need the Fourier inversion formula which is not given in the usual statement of Bochner's theorem. On the other hand, can- not directly use the Fourier inversion formula since it is not obvious how to check the integrability of à when é is mollifier.
Lemma A.15. Suppose & E (RR) is even, bounded, continuous, and k(I,y) := 6(x y) is i.s.p.d. on any compact sets. Then its Fourier transform é is real, nonnegative, and the following inversion formula holds:
((x)=
/
2 for all R
(14)
Proof. The proof adapted from that of Varadhan (2001, Theorem 2.7) and is extended to the multi-dimensional case.
( Since d(E) : 2 (xx) dæ and d(x) = d(-x), with change of variable x' xr, obtain
Ō(C) d(E), so e is real.
Next we show that & is nonnegative. For 0, we compute, for fixed & E R.".
Tm 10,17 11.211
-2 6) pi(t
8) dt
[ |-TTTT"
2-n du
2niu (d() du
T
UII.Iu.1.2T u:l.2T will
2T - 2//ul
II
25ii (()) du
T
|-TTT|"
2
\1=l
|us]]
) T
b( du,
[=]
17

---

Published as a conference paper at ICLR 2023
where have used change of variable t 8, ++ s. Since & E L' (R"), by dominated convergence theorem, we have as T + 0c
I
lim ((E) 2 ((x) dr
2) 8) dl ds.
T-+00 T 10,27 10.37
For 1, 8 R.,, we have
Re (e ( 2 p(t s)))
=cos(2mt E)k(t, 8) cos((2t 8 E) + sin(2nt g)a(t, 8) sin(2 E)-
For fixed T, we define finite measure μ as du [= cos(2 E) dl, then k being i.s.p.d. implies
cos(2nt () 8) cos(2 dt 0,
J0.271 10,27
and similarly for sin. Since is real, conelude that I is nonnegative.
Finally we prove Equation (14). For > O, define (E) = ò(E) -aa--2lle 11el1, Since @ is bounded (because @ E (R"),.)), we see that do E (R.)... We compute, for x E R", using Fubini's theorem,
J @ Sorum520l5 (E) / ò(E)E -0 2xims de
/ (f 2 (()) dy
- 2 de
/ (f 2x1(v 2 E o² 11l13 dE
p(y) dy
/ (m/o 12200///220 + || w1B/ (bv) dy.
where we use the Fourier transform formula (Borodachoy et al., 2019, (4.4.1)) of the Gaussian dis- tribution. Notice that (3)) == 2z/02)nn/2 - is the density of multivariate Gaussian centered at with covariance /²/ (2 ) I. Hence e de = (po ())) Since á is bounded by assumption, with
0 find (E) de
| Taking
0, by monotone
convergence theorem, we have f ((E) de I || so together with the fact that á 2 we have E (RR).. Finally, since do < G. by dominated convergence theorem and Proposition A.4 (note Pa centered at æ), we have
I 2 o(E) dē fim
lim / 2 (5) df = lim (p. o)(0) 6(x).
Proposition A.16. Given d R R. satisfying the assumptions of Lemma A.15, for any Msign (R")), it holds that
[f MT - y) du(x) =
Proof. By Lemma A.15,
II 61 y) du(z) du(u)
(5 dv(x:) du(y)
(
= / -2 2 dv(x))
2 du(e)) 2(E) dg
((E)D(E)O(E) de 1 (=(E) de.
where we use Fubini's theorem (all measures are finite and e 2mi is bounded) to exchange the order of integration.
18

---

Published as a conference paper at ICLR 2023
Proof of Theorem 3.6. First note that Definition A.13(b) is immediate from Theorem 3.3 with fe μ.: {()) = oc, then trivially lim sup,-s (x( & {p)); otherwise we apply Theorem 3.3. So we focus proving criterion Definition A.13(a).
Fix a sequence Me E P(R") that converges weakly to l E P(R"), and our goal to show lim inf,-o (pe) 2 E(u). Without loss of generality, may assume & (e(pe) oo for all (these terms have no effect in lim inf), which implies (X( = 1. By Portmanteau's 's theorem and the assumption that X is closed, we have u(X) limsup Me (X) So all of μe and μ will have support in X.
For m 0, define a nonnegative measure Ve,m by
dve.m (z) = hm (1)p= 1/2 du)((),
where hm R -+ R is continuous monotonically decreasing cutoff function satisfying hm (2) if ||>|| 2 < m and hn (z) = if ||=1| > m + and that hm hm(2)) I hm (a) for m < m'. Then since > 0 is continuous, it is bounded below on any compact set in X, and hence Vem is finite. Also define, for m > nonnegative measure Um by
dum (x) == fim (x)p-1/2)) du(x),
which is again finite following the same reasoning.
Then for any m > 0, denoting dvc(x: == p- (xx) ddcc(a),
(1e) = 1f de y) du (1) dv)
> /ʃ Φr y) due,, (Ix du (u) / m (5) dg,
(15)
where we apply Proposition A.16 for the last equality. On the other hand, note that
De,m - / - due (x) f 2mi 1e -
2
(2(()) (x).
Since Me - weakly and the last integrand is continuous bounded function, we have
lim De,m (E) / 2 a hm (x)p 1/2 (xx du(x) = (E).
c+0
On the other hand, (.) = ʃ --riis-z de (I) dx. Since by Definition 3.1, de dur converges to óg in probability (and in particular weakly). we have lim (e
Applying Fatou's 's lemma to (15), we obtain, for any m > 0,
/ 12m (E) dg.
lif (pce) e
If 1m (5)12 dg for any > 0, then we are done since lim inf,- (p.) = oo. Otherwise, by Kühn (2016. Lemma 2.11), for every m O, V' has density in (R"). This imples μ has density everywhere. Suppose du(z) (xx) (C,. By Plancherel's theorem and the monotone convergence theorem, have
/ pim (E)1 dg
lim inf (He,) 2 lim
m=
( (1mn(=)p11/(())(()))))) dx
lim m-
= (x)g(x)) dz=)) E(v).
This completes the proof of & F E.
Now suppose is compact, and let p == arg mig Ep(.. To establish μz + p² weakly (resp. 1im (x4) E(²") 1), suffices show that every subsequence of (eC) (resp.
19

---

Published as a conference paper at ICLR 2023
(E)() has further convergence subsequence converging μ² (resp. E(2"). With slight abuse of notation, we assume the sequence of is chosen so that { (resp. {E(()))) is already some subsequence of the original sequence. As argued in the proof of Lemma A.12, P(X) is
compact with respect to weak convergence. Hence {u²} has a weakly convergence subsequence
Ve - v for some v E P(X). The r-convergence & 4 & implies
lif Eee (1*) liminf Eee (1T2)) E(v)= 2 lim Eee (v) lim sup &ee (Hez) 1-+00
R=+00
1-000
k-000
where we have used, for each inequality, Mee arg (x) & Definition A. A.)3 the first para- graph of this proof, and again the fact that Me, arg (x) Ee,. Since lim inf Ec, (1*) E(u)) by Theorem 3.3, we have
E(n") lim & (MTt) = E(v) 2 E(1"), k=000
where the last inequality follows because j² the minimizer of Then lima- & (MT,)) E(u)) 1. Moreover, x²(v || v²)) = O. This can only happen if v and μ² agree on all Borel sets, so
A.4 DIFFERENTIAL CALCULUS OF & IN P2()
Lemma A.17. Let Ho R. > and S R. be a compact set. Suppose f : (xo h, zo + h) Ω R is jointly continuous and the derivative f (zo- h, 2o + h) Q + R exists and jointly continuous. Then Sa fx w) dw is differentiable for x E (xo h, + h), and
da J ((r,w)
a
I
- Ou f( w) dw
where the integration is with respect to the Lebesgue measure L.m on Ω
Proof. Fix x E (xo - h, xy + h) and let t > 0 be small enough such that [x 1 f, + t C (ry ho + h). Note that Sn fzz w) dw is well-defined by the dominated convergence theorem since f(I, w) < oo and (2)) < co. Define (ww) = suPre[m-t.o+4] au far w) which is finite since ax (xx, w) is continuous, that ((w) 2 के [f w)r for all x E [r f, x f and ( is integrable since C. (2) < co. Hence by the differentiation lemma (Klenke, 2013, Theorem 6.28) we are done.
0
The following lemma is similar to Korba al. (2021, Proposition 3) but with different assumptions: we do not put any integrability assumptions on W., but we do restrict measures to have compact support. For a differentiable f : R x R + R. we use fif y) to denote the gradient with respect to x. We also use Hif(x, y) to denote the Hessian of f with respect to æ. We similarly denote V y) and fz v)-
Lemma A.18. Assume μ E PP2(R") has density g E (R.).. Let g E C(R":R") Denote s(x, t) = x + tg(x). Then for all t > 0 and all > 0,
de (m)) =2 // F.W. ())))))))))))))))))))))
d
(16)
In particularly,
d {(m)= = /ʃ (E(G) (VW.(,,)) du())
(17)
di
1 1-
Moreover,
d² di2
(14) =2 -2/l // (E(x) vlV + <(x) HH W. (),,))) du(x) du(s).
1=0
(18)
20

---

Published as a conference paper at ICLR 2023
Proof. We compute
(c( = dt (2 lʃ (,(s( t)) du(x)) du(v)
d
Since f has compact support and (x, y, + (.(s(a, s(y, ())(())())))))))) is jointly continuous and its derivative with respect is also jointly continuous, by Lemma A.17, we can push ar the double integral and we obtain (16). Another application of Lemma A.17 shows that if we take derivative with respect again on (16) and evaluate we obtain (18).
d inside
Lemma A.19. Let f E (R"),, p E oo]. Assume Qe has compact support for some 0. Then for all l,..., n, Oe f is differentiable and
0 dc, (o f)(x) (0+
(z).
0x,"
Proof. Since supp(ée) compact and E (R") is bounded on any compact set, know de is well-defined at every x E RR. Note that
① dz
Oxi (ʃ f(x ())
0
(o. f)(x)
(()
!
0f (x (x- ()) dy @.
(x).
Oxg
Since supp(Φ.) is compact and (z f(x y)v(v) is C², by Lemma A.17 we justify the existence the derivative and the exchange of differentiation and integration (?).
L
A.4.1
SUBDIFFERENTIALS OF
Recall the following notion of "Wasserstein gradient" P2(R") from Ambrosio et al. (2005, Definition 10.1.1).
Definition A.20. A vector field w E (2; R*) is a strong Fréchet subdifferential of functional F (2") > - (-00, +000 for all T [); R*), the following holds:
F(T+H) - F(1) / (o(x))(T(())) x) dp(x) + (IIT - - (l
(19)
Note that we cannot apply Ambrosio et al. (2005, Lemma 10.4.1) directly to prove (4) because interaction energies cannot be written in the form of (10.4.1) their setup.
Proof of Proposition 3.7. Let E C (R"; R"). By Lemma A.18, we have E(())= // ([(x) v)) (u(z)d((y)
df t=
= 2 // (<(=) TV, (o.
(0- - vl(p)) g(v) dydp(()
2 {(x) TV (p(z) 11(2))1/((0) g/ (x) du(x),
where the last step follows from applying Lemma A.17 since q has compact support. Now suppose w E (p;; R") is a strong Fréchet subdifferential satisfying (19). For the sequence (Ti}. we have by definition
(x()) - (c( 2 un(x) (x) (T.(()) - x)) du(x)) + (IIT: Taii-Illaase =11
= = wxxx) (1E(x)) du(x:) o(t).
Hence
lim inf (,())) - (1) 2 wxxx) T(()) du(x)) 2 lim inf Es (v(
t
21

---

Published as a conference paper at ICLR 2023
The previous calculation shows d (e exists, and hence it is equal to f (o(x)))(()) du(x). This proves for any & E Ca (R"; R"),
de
/ w(xx) T ((x) du(x)) = / E(x) V (2p(x)-1/2(() * (elo6g///DD (D)(()) du(x))
Hence we have shown (4).
Finally, we show lime-s we (x) = wx²(x) (I) for y.a.e. x under the additional assumption that oe has 2
compact support. By Ambrosio al. (2005, Lemma 10.4.1) with F(x,p(x)) ( - 1) p(x))
we find the strong subdifferential of 2 p²) is given by
wx² (x) = 2 ((x) p(x)
for p-a.e. & E R".
(20)
To show (5), we compute, for y-a.e. T.
(x)xx) v (p(()) 1/2 4/ (z)
(2) q/ (z) + p(z) (122)))-VVV(() * (x) ((p(x)-1/3)b * // vD) (z) + (/3 *
where for the last equality we have applied Lemma A.19. Now taking + 0, by Proposition A.4 using the fact that supp(Ø.) is compact (so that // vp and V (4// /p) bounded x + supp(Ø.) )) we obtain (5).
A.4.2 DISPLACEMENT CONVEXITY OF & J AS
The statement of Proposition 3.8 is similar in form as Korba et al. (2021, Corollary 4) but in our case do not have the second term in Equation (18) vanishing and we need to take the limit -> 0. We also do not resort to RKHS theory in the proof.
Proof of Proposition 3.8. By (18), have
d²
= +.),
1=
where
F. = // (E(x) Vi (x)(x) (p)
G. = // (E(z)) (H, (z)) du (z) dp (y).
We tackle F. first. Observe that successive application of integration by parts using the fact that has compact support gives
8² dx,0y;
F.
Σ +===
W((,,)) dr
a (E(())((z))
Σ
W. (0)) dr dy
dx,
Jy,
5-=1
δ (E(
Σ
v) Oyj
(E, (v)p())) dr
0r.
11=1
Σ
8
Wc (()))()))
Σ
By (((0)())) drdy.
0
dz,
22

---

Published as a conference paper at ICLR 2023
ör, (2)p((() the density of signed measure (it is integrable since it has compact support), and since We is i.s.p.d. on the support of & by Proposition 3.5(a), that each double integral the last expression is non-negative. Hence F. 2 0.
If view Ei=1 2.
Next we show lim Ge O. Since μe has compact support, by Fubini's theorem,
G. = f {(x) ((x) H (WWalss)opl)) dy) {(z)p(z) dr.
To expand the integral inside the Hessian operator, we have
Hi ( Wo
H (p(x)-1/2(0) *VD(
First by the chain rule and Lemma 1.19, we have
d dr (p(x)-1/(6) * vD(a)) dr (p(x)-1/2) (de vp) + p(x) p(2)-1/2
1/2 Øc d dr
vplx))
Diffcrentiating again while applying Lemma A.19, we obtain after rearranging terms,
Hi ( (W.le.w)llo)a))) (v)
d2 dxz
(0. (x)
d de
Φ d
+ 2
p(x)
vD (x)
dr
d²
1/2
vD
(z).
(21)
dr
By Proposition A.4 and the fact that E (R*), we have
lim H.
c+0
dr (c) (x)) vD(x) +2
d2
de D())
dx vD(x) + p(x) 1/2 d?
d
dx vD(x))
(p(x) 1/2 /p(x)) =0. 11/2
2x
Finally, we have
lim G. = lim / {(x) {(x)THL TH,
H. (Fwdlesnhot))) dy {(x)p(z)
c-0
= - {(x)T
(WWreew)rrl)))) dy) ((x))(()) da
lim Hi (
0,
where interchanging the limit and the integral is jusfitied by the dominated convergence theorem and the fact that p and de have compact support (we need compact support assumption of d. make sure convolutions appearing in (21) are uniformly bounded when is sufficiently small).
A.4.3 A DESCENT LEMMA FOR &, WITH TIME DISCRETIZATION
We now consider a time discretization of the gradient flow of & given by, for an initial measure #4 E P2(R") and m Z>o. with a step size 7 O,
Mm+1 = (I - c, ) Mmn
(22)
Similar to Korba et al. (2021, Proposition 14), we can show the following descent lemma for itera- tions (22) with the following additional assumptions:
23

---

Published as a conference paper at ICLR 2023
Assumption A.21. Suppose the target density p (R.) satisfies 1/Cp p(z) < Cp.
Ip(x)]l2 I Cp. [HH(())ll2 S C% for some Cp. Cp: Cp for all I E R", where we use | 2 to
indicate the matrix spectral norm (i.e. ||A]|2 = (max Suppose the mollifier de E C² (R")
satisfies, in addition to Definition 3.1, (c(x) S C. [[Voe( S Co and [H)e( S Co for
some Co. C. C.
Lemma A.22. Under Assumption A.21, there exists a constant such that the function V V.We: R" x Rd + RR is L-Lipschitz in terms of Il-112 in either input.
Proof. Denote r(x) = 111) Then assumptions imply r(z) S C, C. |Vr(z)l2 1 C² = C G, and ||Hr[())]2 < C. = c2 +
f J Cr. We compute, for x, y E R",
V,Wc(z,)) V. (0)(z=)))(=)r()))
V (1 ()
Then
- H)rr())()()) v)(( H1 (,,,( Hob. -))((r))(()))) V r(v) +
and
V2ViW.(Z,v)) = -lo(
Voll - 9)) + 22 - VVV(o
Then we have
HH, 2C,, C.C.Co.. V2 v2. S 2C,,,
Hence conclude that V, W. is L-Lipschitz with
: C² + +C,,C C. C,2).
Remark A.23. Compared with Korba al. (2021, Lemma 1), to ensure W.W, Lipsethiz, we only require uniform boundedness of p and 6. up to second order derivatives instead of up to order 3.
Proposition A.24. Under Assumption A.21, 21, suppose pa E P2(RR) has compact support. Then for any S 2L' with L defined in Lemma A.22 and { Jmzo defined in (22),
((m+)) Ee (pm) S --(1- - 22- We ()(()))) S 0.
(23)
Proof. We first show pm has compact support for m Zz For μu E P2(R") with compact support, the proof of Proposition 3.7 implies (note Proposition 3.7 assumes y has density but it is not necessary to obtain the following formula using the same proof)
We, (x)= 2 V.W.(a, v) du(y)-
Since V W. y) is C³ by the assumptions and y has compact support, by Lemma 171 see that we,a is differentiable with
(x) = H.W. dp(v).
By induction, suppose pm has compact support. Then
supp(pm+1) C Wee,am) y supp(Hm.) supp(um) +y
sup =E
|e l12 (1).
Since W.pp is continuous, the set on the right-hand side is bounded. Hence Mm+1 has compact support.
24

---

Published as a conference paper at ICLR 2023
Now fix Zy0 and we show (23). Define path (pelcejo.11 defined by pr (I-nttwe,arr ) Mm- Let f(t) == (p4). By Lemma A.18, the continuity of We,mm implies that f is differentiable. Moreover, another application of Lemma A.17 implies that f is twice differentiable, and in particular continuously differentiable. Hence f is absolutely continuous on the compact interval [o, 1]. By the fundamental theorem of calculus, we have
E(( 1) - ((1 (()) + S (r(1) (0)) dt.
Observe that by Lemma A.18,
((()) = /f (-yWe, V2)) y) dpm dum
(-yWe,pm (x))
V WiWe(an (2, y) dpm(()) dum (x)
We (I) (x)
where we apply Fubini's theorem to exchange the order of integration together with the fact that pin has compact support and the integrand is continuous. Let s(x, t) ytu (x)). Note that, again by Lemma A.18,
(f((t) - f(()) =2 // [(-rwe (x)) (V. - V (Z B)) dum (x) dum (v)
22y /ʃ |Wee.mm (1) VIW.
W. W., du (x) dy (v)-
Note that, by Lemma A.22, have
( s(y, t)) ViWe(z,))ll
V ViW. s(, V,Wc(s(z, + 11V ViWa(z,v)l.
<T>(([ W,., (x)]l2 (v)ll2))
Thus
(f(()- - f'()) 22 /f ||w0s.mm (z) (1lwe,am (x) m. (1) dlem (x) dum (v)
[1 (emu + (f Ileeemm (2))l2 d(vm(()) 2)
< 2 Tt
S 47 Ltll (((m)))
where we have used the Cauchy-Sudee Schwartz inequality in the last step. Combining everything, we have shown
f( + . (f(tt) ((()) dt
(pm+t)) - Ee =
5 5- - 22:L) We,tem (²((m)) S 0,
1 since 7 < II
B WEIGHTED HYPERSINGULAR RIESZ ENERGY
We recall results most relevant to from Borodachov et al. (2019). Suppose X C R is compact, of Hausdorff dimension d, and d-rectifiable, ie.. the image bounded set in Rd under a Lipschitz mapping. Given non-vanishing continuous probability density p on X, define measure H²
25

---

Published as a conference paper at ICLR 2023
[g p(x) dHa(() for any Borel set B X. The target measure Ry is defined to be h2(B) (2(B)/ The weighted s-Riesz energy is defined as the interaction energy
Fa(wN) = (p(x,)p(zz,) s/2d
(24)
|| Tz l1
The following result states that minimizers of the weighted s-Riesz energy approximate the target measure when the s is sufficienlly large.
Theorem B.1 (Theorem 11.1.2, Borodachov al. (2019)). Suppose d. For w/ arg min E. with w { xN ). we have weak convergence dwž. -hy N
similar result with a slightly different assumption holds for 8 = d (Theorem 11.1.3 of Borodachov al. (2019)).
Comparison of (24) with discretized MIE (6). Our theory is only valid for the full-dimensional case, i.e., n d (see a few exceptions discussed in Remark A. 10).). When this is the case and if the mollifier is taken to be from the s- s-Riesz family, (6) becomes, for s >>
N E (.(wN) = 6-=1 - =2,12 + 8/2
(p(x;) (((x))((2;)))-
Compared with (24), in our case there is an in the denominator, so that the continuous version of the energy & does not blow up all the time (this is in contrast with the continuous version of al. (2019, Theorem 4.3.1)). Moreover, the exponential scaling on p is different: in our case we use 1/2 whereas in (24) it is -s/2n. The diagonal i = j is included our energy, but this inconsequential as another valid discretization to discard the diagonal term (Borodachov et al., 2019, Theorem 4.2.2). On the other hand, our theory allows bigger class of mollifiers that are not necessarily of Riesz families. We also allow X to be non-compact. An interesting future research direction is to extend our theory to cases where d < n, e.g., when X is an embedded d-dimensional submanifold in R"
(24)-see see Borodachov
C ALGORITHMIC DETAILS
In this section we provide algorithmic details of MIED and compare with the updates of SVGD (Liu & Wang, 2016). The negative gradient of (7) with respect 2, is
-v, log Ee(wN) =2
=2
> 371
clus
(IVlogo.)o)- x1) + L log log
ehu
+
Eiyet cts Vlogp(r;)
elay
(2V log ( - x:) + V- Vlogp(a:)).
En, Eiis
where
I5 log Φ. (x; (z,-2;)) (log p(,,) + log
(25)
and get the last equality used the fact that Vo(0) 0 thanks to the assumption o(z) o(-z). Then gradient descent on (7) gives our algorithm (Algorithm 1). The special treatment of the diagonal terms described in Section 4 amounts modifying only the diagonal T
26

---

Published as a conference paper at ICLR 2023
Algorithm 1: Mollified interaction energy descent (MIED) in the logarithmic domain.
Input: target density p. mollifier &e, initial partieles {.]) i==1 learning rate 7), total steps T. for f <- to do
for i < 1 to N do
2 x + n E
(2Vlodde - xa) V Vlogg(ss)) where lig
defined in (25);
end
end
return final particles { Ji=
C.I COMPARISON WITH SVGD
The update formula in Algorithm is similar to the one in SVGD: if use 6. y) in place of the kernel k(z, y) in the SVGD update and rewrite:
N
+1
2
(Vocl a.)) + delag 23 2 Vogg(()
N x Σ (e x;)) (V log - xx) + V Vlgg(zz)))
(SVGD)
For both algorithms, the update formula for each particle consists of attraction and repulsion terms and the total time complexity of each update iteration is ( We note the following differences. First, in our formulation we have scaling factors Σn, e which help stabilizing the optimization (as a by by-product of working in the logarithmic domain) and put more weight on nearby particles as well as particles in low-density regions, whereas (SVGD) the scaling factors are de (2, - xx)) which are not adapted to prioritize low-density regions. Second, in MIED, the attraction force for particle i only comes from V log p(x:), whereas in (SVGD) the attraction force comes from V log p(x,) for all j's. Third, for each j. in formulation the repulsive force has additional factor of 2 in front of V log Φ. (xx z2).
Empirically, the additional scaling factors in MIED help produce samples with good separations compared to SVGD, since closer pairs of points will have large weights. Additionally, since MIED optimizes dimen-dimeniionll objective (7),
employ accelerated gradient-based based optimizers like Adam (Kingma & Ba, 2014), which we used in our experiments. In contrast, SVGD does optimize any finite-dimensional objective. While practical SVGD implementations also optimizers like Adam. it is unclear how the resulting particle dynamics is related to the gradient flow of KL divergence.
C.2 HANDLING CONSTRAINTS WITH DYNAMIC BARRIER METHOD
The dynamic barrier method is introduced in Gong & Liu (2021) which solves min, f(x) subject to g(x)) < where y valued, Intuitively, their method computes update directions by either decreasing g(x) when g(z) > 0, following fz) if the constraints satisfied, balancing both gradient directions.
In order to handle multiple constraints such as in Figure/ 11, consider generalized version of their dynamic barrier method. In this generalized setting, : R" + R is vector-valued and constraints are (x) < where the S sign interpreted coordinate wise.. Suppose are at iteration with current solution x' Then the next update direction vz is taken to be the arg min of
min | - V[(")] s.t. Vi -i= l..., m, (9i("' aggi(x"),
(26)
where ag 0 are fixed hyperparameters; in our implementation simply choose a, 1. Then 1-- x² net with learning rate 17. In our implementation we use Adam (Kingma & Ba, 2014) that modulates the update directions. Observe that the optimization problem (26) is the
27

---

Published as a conference paper at ICLR 2023
same projecting a point fx") onto the polyhedron formed by the intersection of the halfspaces n {z E R Vg.(=') 2 azgi(z")). To solve (26), we use Dykstra's algorithm (Tibshi- rani, 2017) which can be interpreted rumning coordinate descent fixed number of 20 iterations for the Dykstra's algorithm which we found to be sufficient for our experiments; in the of single constraint, only need to use iteration.
the dual of (26). We use
D EXPERIMENT DETAILS AND ADDITIONAL RESULTS
D.1 GAUSSIANS IN VARYING DIMENSIONS
We generate the n-dimensional Gaussians used to produce Figure as follows. We generate matrix A VAVA / det (VA) This way (A) 1. We then use A as the covariance matrix for the Gaussian (centered at 0). We use Adam with learning rate 0.01 for all methods for a total of 2000 iterations. This is enough for SVGD and MIED to converge, while for KSDD the convergence can be much slower.
V E R with i.i.d. entries uniformly in [-1, 1]. Then
We visualize the samples from each method for n = in Figure 5. We notice that MIED capable of generating well-separated samples, while for SVGD there is gap between the inner cluster of samples and a sparse outer ring. For KSDD we see the artifact where too many samples concentrate on the diagonal.
K500
Figure 5: Visualization of samples from each method for a 2D skewed Gaussian.
D.2 COLLAPSED SAMPLES WHEN THE KERNEL WIDTH TOO BIG
In Figure 3, we see that samples from SVGD collapse with adaptive kernel where the variance taken to be half of the median of the squared distance among all pairs of points (Liu & Wang, 2016); at termination the median of the squared distance is greater than 1 in that experiment. Here we investigate this issue further. In Figure 6, for the same uniform sampling setup, we visualize the results of SVGD with fixed kernel width and MIED with Gaussian mollifiers with the same kernel width: when the kernel width (i.e. 2e in (x)) exp(
/29) is too big, both SVGD and MIED result in collapsed samples. This is because since the target density is constant, the only force in the updates is the repulsive force. When the kernel width is too large, the repulsive force coming from points in the same collapsed cluster is dominated by the repulsive force coming from points from other clusters-this evident in the update directions shown in the leftmost column of Figure 6 (SVGD with kernel width 0. When using the dynamic barrier method Gong & Liu (2021) to enforee the square constraints instead of reparamcterization with tanh, we obtain similar results as in Figure 6.
| 1 2x
This pathological phenomenon not only limited to constrained sampling: when sampling the 2D Gaussian from Figure using too big kernel width also result in collapsing (Figure 7).
We emphasize that our theory of MIED suggests that in practice we need to choose the kernel width very small in order sample from the correct target measure according to Theorem 3.3 and Theorem 3.6. In comparison, the theory of SVGD has no such implication.
D.3 UNIFORM SAMPLING WITH AN ALTERNATIVE MIRROR MAP
In this section, show that for sampling from uniform distribution in the square [-1, 1|2, the results of SVMD/MSVGD (Shi et al., 2021) depend heavily on the choice of the mirror map. In- stead of the entropic mirror map used produce results in Figure 3, here the mirror map
28

---

Published as a conference paper at ICLR 2023
MED0H
Figure 6: SVGD and MIED with fixed-size Gaussian kernels for uniform sampling in square. In each cell of the grid, we plot the samples along with the update direction (black arrows) at that iteration. Rows correspond to iterations 100, 200, 1000. Columns correspond to each method with varying kernel widths (twice the variance of the Gaussian kemnel/mollifier) indicated in the paren- theses.
MEDOO
Figure 7: SVGD and MIED with fixed-size Gaussian kernels for sampling a 2D Gaussian in
Figure 5. We see that using too big kernel size can lead to collapsed samples for both SVGD and MIED.
6(0) -t SVMD/MSVGD fail draw samples near the boundary; we suspect this because the gradient of the conjugate VΦ*(n) = (/1+72-1//, (coordinate-wel wise arithmetic) requires coordinates of n to go to to land near the boundary. We verify this phenomenon by using Vø" (7)) the reparametriza- tion map in MIED (rightmost figure Figure 8): indeed with such reparameterization MIED also struggles the boundary.
i-1 log pー + log T+ as in Ahn & Chewi (2021). The results are shown in Figure 8.
29

---

Published as a conference paper at ICLR 2023
SVMD
MSVGD
MIED Riesz) w/ mirror
1.0-
1.0-
1.0
0.5
0.5
0.5
0.0
Figure 8: Visualization of samples for uniform sampling from 2D box when using suboptimal
mirror map. All three methods fail to draw samples near the boundary of the box [_1,
D.4 20-DIMENSIONAL DIRICHLET DISTRIBUTION
We sample from the 20-dimensional Dirichlet distribution in the same setup in Shi al. (2021) with 50 particles. Results and visualization are shown in Figure 9. We see that unlike sampling from box (Figure 3), both MSVGD and SVMD by Shi et al. (2021) perform well in this experiment. This is due to the fact that the entropic mirror map used here is well-tested choice for simplex constraint, yet obtaining good mirror map for a generic constraint, even if it is linear, can be challenging. Our method does not have such a limitation, as it can easily incorporate existing constrained optimization tools.
SYMD
If
SVDD
MOVDD
MITT/Beo)
MEEX man)) MVED(aaplon
MIED mmen)
12"1
307
9//0
MSW77)
18 We0e
Mi MFD0/oomsoe
M FFD02
6105
610
605
Figure 9: Top: Metrics vs. the number of iteration for sampling the 20-dimensional Dirichlet distri- bution. Bottom: visualization of samples from each method.
D.5 EFFECT OF FOR RIESZ MOLLIFIERS
When we use the s-Riesz families of mollifiers in MIED, we have the freedom of choosing the hyperparameter so long as s n. In this section, we study the effect of on the resulting samples. We consider the problem of sampling from mixture of four 2D Gaussians centered at (#1, =1), each with diagonal variance 0.3 and constrained to the [-1, box. We vary number of particles N in {100, 200, 500, 1000, 20003. All runs use total of 1000 iterations with learning rate 0.01. In the top of Figure 10, we plot the W2 distance and energy distance as functions of for each N. Interestingly, we see the best performing is in [3, 5.0] and depends on N. This suggests that our choice of = n + 10 in Section 5 may not be optimal and there is room for hyperparameter tuning to further improve the performance of MIED with Riesz kernel. At the bottom of Figure 10 we visualize the samples of MIED with s = and of SVGD with kernel width 0.01 (adaptive kernels would result in collapsed samples and other widths tested on would result
[2, 10] and the
30

---

Published as a conference paper at ICLR 2023
in worse samples). While SVGD samples form visible artifacts, the samples of MIED are evenly distributed and the four modes of the mixtures emerge as N increases.
10-1
Energy Distance
W; Distance
10
N = 100
N = 100
10-2
N 200
= 500
10
500
N 1000
N 1000
2000
2000
3
10-1
5.0
7.5
10.0
12.5
15.0
2.5
5.0
7.5
10.0
12.5
15.0
100
500
N 1000
N 2000
100
100
1.0
1.0
1.0
0.5
0.5
0.5
0.5
0.5
0.0
0.0
0.0
0.0
0.0
-0.5
-0.5
-0.5
-0.5
-0.5
-100
1.00
1.00
-.0
-1.0
N 100
. = 200
2V 500
2V 1000
NV 2000
Figure 10: Ablation study of the hyperparameter for MIED with s-Riesz families of mollifiers. Top: metrics as functions of 8. Bottom: visualization of samples MIED with Riesz mollifier with and of SVGD with kernel width 0.01.
D.6 MORE CONSTRAINED SAMPLING EXPERIMENTS
In this section we test MIED on more low dimensional constrained sampling problems and qualita- tively assess the results. Note that mirror LMC (Zhang et al., 2020; Ahn & Chewi, 2021) or mirror SVGD (Shi et al., 2021) cannot be applied due non-convexity of the constraints. In Figure 11, we consider uniform sampling of a challenging 2D region with initial samples drawn from the top- right comner the number of fterations increases, MIED gradually propagate samples to fill up the entire region. In Figure 12, we consider sampling from a von Mises-Fisher distribution on unit sphere. Although our theory focuses on sampling from full-dimensional distribution, as discussed in Remark A. 10, we can extend Theorem 3.3 to the case of a sphere due to its symmetry. We see the samples visualized in Figure 12 capture the two modes that emerge by restricting the Gaussian to unit sphere.
D.7 DETAILS ON FAIRNESS BAYESIAN NEURAL NETWORK EXPERIMENT
We use 80% 80% training/test split as Liu et al. (2021). We use the source code provided by Liu et al. (2021) with default hyperparameters. The source code provided by the authors of Liu et al. (2021) does not implement the calculation of g() faithfully as written in the formula, so we corrected it. All methods use 2000 iterations for training- For our method we leaming rate 0.001. One of their four methods (Control+SVGil got stuck at initialization (with accuracy around 0.75), so we omit its result from the plot in Figure
31

---

Published as a conference paper at ICLR 2023
62
Figure 11: Uniform sampling of the region {x, y) E [-1,
(cos(3TTz) + cos(3Ty))² < 0.3} using MIED with Riesz mollifier 3) where the constraint is enforced using the dynamic barrier method. The plot in row i column shows the samples at iteration 100 + 200(62 + j). The initial samples are drawn uniformly from the top-right square [0.5, 1.0
Figure 12:
Sampling from the von Mises-Fisher distribution obtained by constraining the 3- dimensional Gaussian from Appendix D.1 to the unit sphere. The unit-sphere constraint is enforced using the dynamic barrier method and the shown results obtained using MIED with Riesz kemell and 3. The six plots are views from six evenly spaced azimuthal angles.
32


Reviews:
[{"review_id": "uBQrbld-0kp", "reviewer": "Reviewer_cXXg", "summary": "", "questions": "", "limitations": "", "rating": 8, "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": "", "summary_of_the_paper": "The paper introduces algorithms for sampling constrained distributions. The authors introduce an algorithm MIED that transform this sampling problem into an optimization problem, and practically, to optimize the mollified interaction energy in a first-order particle-based solution. They also proved that by minimizing the energy function, the MIE coverges to the target measure in the sense of chi-square divergence. The authors experimentally validate their algorithms on simulated data to do unconstrained sampling and constrained samplign.", "summary_of_the_review": "This paper addresses an interesting problem and its development looks very exciting. The work is very complete and addresses both theoretical and experimental aspects of the proposed methods. The theory proposed looks sound with a strong background and the experiments look convincing. ", "strength_and_weaknesses": "Strengths:\n\n- The motivation of the paper is clear and interesting: sampling from a constrained domain either contain expensive numerical subroutines or requires explicit formulas for quantities, which are invalid in most of the cases. \n-The paper is overall well-written and the authors tries to make the main idea easy to follow.\n\nWeakness:\n- Since the performance in uncontrained sampling is close to the previous methods, the authors could put more emphasize on the contrained sampling part to show the effectiveness of the method. Some of the constrained sampling results could be moved to the appendix.\n- How to choose the hyperparameter s in practice? Is the method sensitive to s?\n\nTypo:  $\\int_x f(x)^p >0$.\n", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "nFoQctvoZNE", "reviewer": "Reviewer_oXmc", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "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": "", "summary_of_the_paper": "This paper presents a new particle based probabilistic approximate inference method that (similar to Stein Variational Gradient Descent algorithm) is based on an optimisation task involving attractive and repulsive forces. In the present work, these forces are defined differently and via a function called mollified interaction energy (MIE). \nMinimising the proposed MIE function involves balancing a repulsive force to increase the distance between the particles and  an attractive force that pushes the particles to high-density regions. It is argued that the presented work might be suitable for constrained density functions. \n", "summary_of_the_review": "This paper presents a new particle based inference method that minimised a function called MIE. \nThe algorithm is interesting but the experimental results are not thorough (see the previous section).", "strength_and_weaknesses": "Strength: Up to my knowledge the presented algorithm is novel and correct. The paper is overall well-written (though, possibly the measure theoretic jargon could be explained simpler and more intuitively to increase the impact).\n\nWeaknesses: \n1. Most their experimental models are low-dimensional or toy models. It would add to the value of the paper if they could show some results on more complicated models (e.g. multimodal models) and high-dimensional models (to represent the role of the dimension on the attractive/repulsive forces.  \n2. The importance (and novelty) of their theoretical findings is not clearly explained.  \n3. The convergence of the proposed algorithm to the target distribution is not guaranteed meaning the algorithm is not asymptotically unbiased. \n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and 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": {}}, {"review_id": "aVte0WGeiE", "reviewer": "Reviewer_FquR", "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": "", "summary_of_the_paper": "This paper studies the problem of sampling from a target density without knowing its normalization constant. It focuses on one class of variational inference methods called Stein variational gradient descent (SVGD), which formulates the sampling problem into an optimization problem that minimizes the KL distance between the target distribution and a distribution that defines the mollified interaction energy. Based on it, the authors propose the mollified interaction energy descent (MIED) method to run gradient descent on the minimization problem. They proved that the exact version of MIED asymptotically converges to the target distribution. Empirical evaluations are provided to compare MIED with previous methods.\n", "summary_of_the_review": "As described in the previous sections, I think this paper is studying a very interesting problem and the algorithm proposed here seems to be competitive over existing methods, especially SVGD. However, it is not discussed in full detail what the difference between them at an algorithmic level is. It would be better to have some discussion on the difference between the algorithm design in SVGD and MIED in the algorithm section (Section 4).\n \nMoreover, the theoretical results in this paper are all in the asymptotic sense. Moreover, the problem in (2) is not directly solvable and the paper proposes to minimize a discrete version of (2) in practice instead, which brings additional approximation error. It is not clear how large these errors could affect the sampling error in the final results. It is a natural question that does the proposed algorithm enjoy any non-asymptotic convergence? This differs from the LMC based approaches, where finite-time and fast convergence could be established for sampling from log-concave and non-log-concave densities (see [1, 2, 3] for example). Then what are the difference and advantage of using MIED instead of other LMC methods?\n\n[1] Dalalyan, Arnak S. \"Theoretical guarantees for approximate sampling from smooth and log concave densities.\" Journal of the Royal Statistical Society: Series B. 2017.\n[2] Zou et al. \"Faster convergence of stochastic gradient Langevin dynamics for non-log-concave sampling.\" Uncertainty in Artificial Intelligence. 2021.\n[3] Muzellec et al. \"Dimension-free convergence rates for gradient Langevin dynamics in RKHS.\" Conference on Learning Theory. 2022.\n\nIn Figure 3, why do the samples from SVGD collapse?\n\nAlthough the experiments show that MIED can deal with different constraints unlike SVGD, could you explain what causes the performance drop of MIED in Figure 6?\n\nWhat is the dataset used in the training of fair Bayesian neural networks?\n", "strength_and_weaknesses": "This paper is well-written and clearly structured. The proposed method is more general than the previous ones compared in the paper. The experimental results also show the advantage of MIED in constrained and unconstrained sampling problems. The main weakness of the paper is its lack of discussion about the comparison of MIED with LMC-based algorithms both theoretically and empirically. The theoretical results in this paper are mainly asymptotic and it is not clear whether MIED will have good convergence properties in various empirical applications. \n", "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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "X6xT9mcM8nJ", "reviewer": "Reviewer_gruH", "summary": "", "questions": "", "limitations": "", "rating": 8, "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": "", "summary_of_the_paper": "The paper proposes a new method to sample from a target distribution that is known up to a normalizing constant. \nIn essence, the proposed algorithm updates a set of particles by gradient descent and the loss is a smart pick that asses the energy of the pairwise interactions between the particles. When it's small, the particles follow the target distribution by design.\nAs I see it, the way to do it is to \"smooth\" the target distribution through some kernel with a bandwidth \\eps, for which the problem is tractable, and then see that the smaller the bandwidth get, the close we get to the target. This \"smoothing\" operation is called \"mollifying\".\n\nAll in all, the paper is nice to read, provides a lot of useful pointers and the proposed method looks surprisingly easy to implement. ", "summary_of_the_review": "As a short summary, I think the authors do propose a principled and nice method to sample from a target distribution whose density is known up to a constant, that looks like an interesting alternative to MCMC. The paper is very inspiring to me, and leaves me with some comments. \n\nSince I don't exactly know where to write them, here I list them below in no particular order\n* 2/ pairwise interaction: you should maybe mention that $p$ is still the target here, or maybe it's obvious\n* looks to me you forgot |f(x)|^p instead of f(x)^p in the notation for ||f||_p\n* definition 3.2 looks super central. However, I was disapointed there was so few discussion about it, at least to provide some \"feel\" or rationale for it. It actually looked counter intuitive until I understood you were about to _minimize_ it, which was not obvious to me at first. Then I had to check appendix A.2 to see it in action.\n* As some food for thought and just to know your opinion on this (in this thread, not necessarily in the paper), I must ask whether you think you may come up with a way to define that kind of energy for any particular divergence between distributions ? Here we have chi^2, but it looks like maybe you may find some trick to generalize your idea to any other one, or at least ones with some features to be identified.\n* How restrictive is this assumption about ispd ? I am not used to it, and I wonder what it means in practice that may limit my choice of \"mollifier\"\n* you really don't want to write just kernel ? The historical reference is nice but do we need another name ? Is it used in some community ?\n* In equation (7), do we need N^2 terms ? Isn't there some symmetry that allows us using only N^2/2 ?\n* if N is too large, could you consider mini batches of pairs ?\n* the algorithm is nice, but do we know when to stop ? Do we have some clue whether when the samples actually match the target distribution ?\n* The \"special treatment of the diagonal terms\" is definitely not the most elegant part of the paper. Do you think you could handle the diagonal and off-diagonal terms in a different way, maybe by balancing them in some way ?\n* In this paragraph: \"which can dominant\", \"for mollifiers we consider\" -> \"for the mollifiers we consider\"\n* I have no clue about the \"dynamic barrier method\". For the paper to be self contained, it could help at least describing it a bit.\n* page 9 \"the mirror map has to chosen\" -> to be chosen\n* conclusion: \"the target mesaure\" -> measure\n", "strength_and_weaknesses": "Strengths of the paper include:\n* The method is well grounded theoretically and although I could not follow the proofs in detail, I trust the authors on this\n* The paper is a very good balance between motivations, background, a clear explanation of the method with all the hard maths in appendix, and a very nice part about how to implement the idea in practice\n* I personally found the paper very inspiring and it leaves me with many questions / ideas that I will just write down later here.\n* The method seems to work well and to offer a lot of ways of improvement.\n\nOn the downsides:\n* Since the method seems to work just as well as all others (see the experiments) on many aspects, I am wondering what would be the interesting features ? I must say that I like it that a paper is not bragging about how better it is everywhere, but I would appreciate it if the authors could at least provide some hints regarding when it is / is not appropriate. What about computing time ?\n* I am not used to this kind of experiments and I don't exactly know if such low dimensional problems may be considered state of the art. If so, it's all right. But my core concern is whether a particle-based system would really be able to scale to high dimensions, since it would maybe mean an extreme large number of particles ? How would that work ?", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and 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": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster