# Sampling with Mollified Interaction Energy Descent

Lingxiao Li
MIT CSAIL
<EMAIL_ADDRESS>&Qiang Liu
Computer Science Department
The University of Texas at Austin
<EMAIL_ADDRESS>&Anna Korba
CREST, ENSAE
IP Paris
<EMAIL_ADDRESS>Mikhail Yurochkin
IBM Research, MIT-IBM Watson AI Lab
<EMAIL_ADDRESS>&Justin Solomon
MIT CSAIL
<EMAIL_ADDRESS>
###### 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.

##  1 Introduction

Sampling from an unnormalized probability density is a ubiquitous task in
statistics, mathematical physics, and machine learning. While Markov chain
Monte Carlo (MCMC) methods (Brooks et al., 2011) provide a way to obtain
unbiased samples at the price of potentially long mixing times, variational
inference (VI) methods (Blei et al., 2017) approximate the target measure with
simpler (e.g., parametric) distributions at a lower computational cost. In
this work, we focus on a particular class of VI methods that approximate the
target measure using a collection of interacting particles. A primary example
is 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.

While MCMC and VI methods have found great success in sampling from
unconstrained distributions, they often break down for distributions supported
in a constrained domain. Constrained sampling is needed when the target
density is undefined outside a given domain (e.g., the Dirichlet
distribution), when the target density is not integrable in the entire
Euclidean 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)). A few recent approaches (Brubaker et al., 2012;
Byrne & Girolami, 2013; Liu & Zhu, 2018; 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 on a case-by-case basis from the constraints.

We present an optimization-based 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 a double
integral of the quotient of a _mollifier_ —smooth approximation of
$\delta_{0}$, 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 $\delta_{0}$, the MIE converges to the
$\chi^{2}$ divergence to the target measure up to an additive constant and the
gradient flow of the MIE agrees with that of the $\chi^{2}$ divergence,
providing theoretical basis for using MIED as a sampling method.

While mollifiers can be interpreted as kernels with diminishing bandwidths,
our analysis is fundamentally different from that of SVGD where a fixed-
bandwidth kernel is used to define a reproducing kernel Hilbert space (RKHS)
on which the Stein discrepancy has a closed-form (Gorham & Mackey, 2017).
Deriving a 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 the same form for constrained and unconstrained
domains and approximates the $\chi^{2}$ divergence so long as the bandwidth is
sufficiently small so that short-range interaction dominates the energy.

Algorithmically, we use first-order optimization to minimize MIEs discretized
using particles. We introduce a log-sum-exp trick to 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 reparameterization using a differentiable (not necessarily bijective)
map or the dynamic barrier method by Gong & Liu (2021) to handle a generic
differentiable inequality constraint. Our method is effective as it only uses
first-order derivatives of both the target density and the inequality
constraint (or the reparameterization 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 gradient flow of the Kullback-Leibler (KL) divergence has been
extensively studied, and many popular sampling algorithms can be viewed as
discretizations of the KL-divergence gradient flow. Two primary examples are
Langevin Monte Carlo (LMC) and Stein variational gradient descent (SVGD). LMC
simulates Langevin diffusion and can be viewed as a forward-flow splitting
scheme for the KL-divergence gradient flow (Wibisono, 2018). At each iteration
of LMC, particles are pulled along $-\nabla\log p$ where $p$ is the target
density, while random Gaussian noise is injected. In contrast with LMC, SVGD
is a deterministic algorithm that updates a collection of particles using a
combination of an attractive force involving $-\nabla\log p$ 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 $\chi^{2}$ divergence (Chewi et 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 et al., 2022; Korba et al.,
2020; Salim et al., 2022).

Sampling in constrained domains. Sampling in constrained domains is more
challenging compared 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-Jacobian term that can be costly to evaluate.
Brubaker et al. (2012) propose a constrained version of HMC for sampling on
implicit submanifolds, but their algorithm is expensive as they need to solve
a nonlinear system 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 mirror-Langevin
diffusion for constrained sampling provided that a 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 incorporate 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 $E(\\{x_{i}\\}_{i=1}^{N})=\sum_{i\neq
j}W(x_{i},x_{j})$ for a kernel $W(x,y)$. The gradient flow of $E$ on $N$
particles can give rise to phenomena such as flocking and swarming, and a long
line of mathematical research focuses on mean-field convergence of the
particles following the gradient flow of $E$ to continuous solutions of
certain PDEs as $N\to\infty$ (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)=\norm{x-y}^{-s}_{2}$ with sufficiently big $s$, minimizing $E$ over a
compact set yields uniform samples as $N\to\infty$. 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 methods. In contrast,
we take a different approach by using integrable kernels in MIEs through
mollifiers.This enables us to apply variational analysis to show our
convergence results. We compare our formulation with theirs further in
Appendix B. Recently, Korba et al. (2021) propose sampling via descending on
kernel Stein discrepancy which can also be viewed as a type of interaction
energy, but like SVGD, it is limited to unconstrained sampling and can be slow
as higher-order derivatives of $\log p$ are required. Current work Craig et
al. (2022) propose approximating the $\chi^{2}$ 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 at each step which can be costly. Consequently their method is only
applied to special target densities for which the integral has a closed-form.
In comparison, our method can handle general target densities with cheap
update in each step.

##  3 Mollified interaction energy

Notation. Let $X\subset\mathbf{R}^{n}$ be the domain we work in. We use
$\mathcal{L}_{n}$ to denote the Lebesgue measure on $\mathbf{R}^{n}$ and write
$\int f(x)\mathop{}\\!\mathrm{d}x$ to denote integration with respect to
$\mathcal{L}_{n}$. We will assume all measures are Borel. We use
$\mathcal{H}_{d}$ to denote the $d$-dimensional Hausdorff measure. We use
$\omega_{N}$ to denote a set of points
$\\{x_{1},\ldots,x_{N}\\}\subset\mathbf{R}^{n}$. For $x\in\mathbf{R}^{n}$, let
$\delta_{x}$ be the Dirac delta measure at $x$ and $\delta_{\omega_{N}}$ be
the empirical measure $\frac{1}{N}\sum_{k=1}^{N}\delta_{x_{k}}$. We denote
$L^{p}(X):=\\{f:X\to\mathbf{R}\text{ Lebesgue measurable with
}\int_{X}\absolutevalue{f(x)}^{p}\mathop{}\\!\mathrm{d}x<\infty\\}$, and for
$f\in L^{p}(X)$ we let
$\norm{f}_{p}:=\left(\int_{X}\absolutevalue{f(x)}^{p}\mathop{}\\!\mathrm{d}x\right)^{1/p}$.
We use $\norm{f}_{\infty}$ to denote the essential supremum of
$\absolutevalue{f}$. For $f,g:\mathbf{R}^{n}\to\mathbf{R}$, we define
convolution $(f*g)(x):=\int f(y)g(x-y)\mathop{}\\!\mathrm{d}y$ provided this
integral exists. We use $C^{k}(X)$ to denote continuously $k$-differentiable
functions and $C^{\infty}(X)$ to indicate the smooth functions. The space of
continuous $k$-differentiable functions with compact support on $X$ is
$C_{c}^{k}(X)$. We denote by $\mathcal{P}(X)$ the set of probability measures,
and $\mathcal{P}_{2}(X)$ the set of probability measures with bounded second
moments. We denote by $\mathcal{M}_{\textup{sign}}(X)$ the set of finite
signed measures. Given a Lebesgue measurable map $\bm{T}:X\to X$ and
$\mu\in\mathcal{P}_{2}(X)$, ${\bm{T}}_{\\#}\mu$ is the pushforward measure of
$\mu$ by $\bm{T}$.

#### Setup.

The domain $X\subset\mathbf{R}^{n}$ where we constrain samples in is assumed
to be closed and full-dimensional, i.e., $\mathcal{L}_{n}(X)>0$; the
unconstrained case corresponds to $X=\mathbf{R}^{n}$. The target density,
always denoted as $p$, is assumed to satisfy $p\in C^{1}(X)$, $p(x)>0$ for all
$x\in X$, and $0<\int_{X}p(x)\mathop{}\\!\mathrm{d}x<\infty$. Let $\mu^{*}$ be
the target probability measure
$\mu^{*}(B):=\nicefrac{{\int_{B}p(x)\mathop{}\\!\mathrm{d}x}}{{\int_{X}p(x)\mathop{}\\!\mathrm{d}x}}$
for Borel $B\subset X$. Our goal is to sample from $\mu^{*}$.

###  3.1 Mollifiers

Our framework relies on the notion of _mollifiers_ originally introduced in
Friedrichs (1944).

######  Definition 3.1 (family of mollifiers).

We say $\\{\phi_{\epsilon}\\}_{\epsilon>0}\subset C^{1}(\mathbf{R}^{n})$ is a
family of _mollifiers_ if it satisfies

  1. (a)

For any $\epsilon>0$, $x\in\mathbf{R}^{n}$, $\phi_{\epsilon}(x)\geq 0$ and
$\phi_{\epsilon}(x)=\phi_{\epsilon}(-x)$.

  2. (b)

For any $\epsilon>0,$ $\norm{\phi_{\epsilon}}_{1}=1$ and
$\sup_{x\in\mathbf{R}^{n}}\phi_{\epsilon}(x)<\infty$.

  3. (c)

For any $\delta>0$, $p\in\\{1,\infty\\}$, $\lim_{\epsilon\to
0}\norm{\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}}_{p}=0$.

As $\epsilon\to 0$, the distribution with density $\phi_{\epsilon}$ converges
to $\delta_{0}$, and if $f\in L^{p}(\mathbf{R}^{n})$ and is continuous, then
convolution $\phi_{\epsilon}*f$ converges to $f$ both pointwise (Proposition
A.4) and in $L^{p}$ (Proposition A.7).

Our definition is different from the one used in PDE theory (Hörmander, 2015),
where some $\phi\in C_{c}^{\infty}(\mathbf{R}^{n})$ is fixed and the family of
mollifiers is generated by
$\phi_{\epsilon}(x):=\epsilon^{-n}\phi(x/\epsilon)$. While in PDE theory
mollifiers are typically used to improve the regularity of non-smooth
functions, in our framework they are used to construct interaction energies
that approximate the $\chi^{2}$ divergence.

We will use the following mollifiers. In following definitions we include
normalizing constants $Z_{\epsilon}$ that ensure $\norm{\phi_{\epsilon}}=1$.
We do not write out $Z_{\epsilon}$ explicitly as they might not admit clean
forms and they are only relevant in theory but not in our practical algorithm.
For $s>n$, the _$s$ -Riesz family of mollifiers_ is defined as
$\phi^{s}_{\epsilon}(x):=(\norm{x}^{2}_{2}+\epsilon^{2})^{-s/2}/Z^{s}_{\epsilon}$.
The _Gaussian family of mollifiers_ is defined as
$\phi^{g}_{\epsilon}(x):=\exp\left(-\frac{\norm{x}_{2}^{2}}{2\epsilon^{2}}\right)/Z_{\epsilon}^{g}$.
The _Laplace family of mollifiers_ is defined as
$\phi^{l}_{\epsilon}(x):=\exp\left(-\frac{\norm{x}_{2}}{\epsilon}\right)/Z_{\epsilon}^{l}$.
Since $\\{\phi_{\epsilon}^{g}\\}$ and $\\{\phi_{\epsilon}^{l}\\}$ correspond
to the densities of centered Gaussian and Laplace random variables, they
satisfy Definition 3.1. Proposition A.1 proves that
$\\{\phi_{\epsilon}^{s}\\}$ also satisfies Definition 3.1.

###  3.2 Mollified Interaction Energies

######  Definition 3.2 (mollified interaction energy).

Given a family of mollifiers $\\{\phi_{\epsilon}\\}_{\epsilon>0}$ satisfying
Definition 3.1, for each $\epsilon>0$, define a symmetric kernel
$W_{\epsilon}:\mathbf{R}^{n}\times\mathbf{R}^{n}\to[0,\infty]$ by

$W_{\epsilon}(x,y):=\left\\{\begin{array}[]{ll}\phi_{\epsilon}(x-y)(p(x)p(y))^{-1/2}&\text{
if }x,y\in X\\\ \infty&\text{ otherwise.}\end{array}\right.$ (1)

Define the _mollified interaction energy_ (MIE),
$\mathcal{E}_{\epsilon}:\mathcal{P}(\mathbf{R}^{n})\to[0,\infty]$, to be

$\displaystyle\mathcal{E}_{\epsilon}(\mu)$ $\displaystyle:=\iint
W_{\epsilon}(x,y)\mathop{}\\!\mathrm{d}\mu(x)\mathop{}\\!\mathrm{d}\mu(y).$
(2)

Intuitively, minimizing $\mathcal{E}_{\epsilon}(\mu)$ balances two forces: a
repulsive force from $\phi_{\epsilon}(x-y)$ makes sure $\mu$ has a good
spread, and an attractive force from $(p(x)p(y))^{-1/2}$ helps $\mu$
concentrate on high-density regions. The exponent $-1/2$ is chosen to balance
the two forces so that, as we will show in Theorem 3.3,
$\mathcal{E}_{\epsilon}(\mu)$ approximates the $\chi^{2}$ divergence with
respect to $\mu^{*}$ for small $\epsilon$.

###  3.3 Convergence to $\chi^{2}$-divergence

Recall that the $\chi^{2}$-divergence between probability measures $P,Q$ is
defined by

$\displaystyle\begin{split}\chi^{2}(Q\;\|\;P):=\left\\{\begin{array}[]{ll}\int\left(\frac{\mathop{}\\!\mathrm{d}Q}{\mathop{}\\!\mathrm{d}P}-1\right)^{2}\mathop{}\\!\mathrm{d}P&\text{
if }Q\ll P\\\ \infty&\text{ otherwise,}\end{array}\right.\end{split}$ (3)

where $Q\ll P$ denotes absolute continuity of $Q$ with respect to $P$ with
Radon-Nikodym density
$\nicefrac{{\mathop{}\\!\mathrm{d}Q}}{{\mathop{}\\!\mathrm{d}P}}$.

######  Theorem 3.3.

Suppose $\mu\in\mathcal{P}(\mathbf{R}^{n})$ satisfies
$\chi^{2}(\mu\;\|\;\mu^{*})<\infty$. Then,
$\mathcal{E}_{\epsilon}(\mu)<\infty$ for any $\epsilon>0$. Furthermore,

$\displaystyle\lim_{\epsilon\to
0}\mathcal{E}_{\epsilon}(\mu)=\chi^{2}(\mu\;\|\;\mu^{*})+1.$

We provide a succinct proof of Theorem 3.3 in Section A.2 using the theory of
mollifiers developed in Section A.1. In Remark A.9 we discuss extension of
Theorem 3.3 to cases when $X$ is not full-dimensional; in particular Theorem
3.3 still holds when $X$ is “flat”, i.e., has Hausdorff dimension $d<n$ and is
contained in a $d$-dimensional linear subspace.

###  3.4 Convexity and $\Gamma$-convergence

We next study properties of the minimizers of MIEs: Does
$\min_{\mu\in\mathcal{P}(X)}\mathcal{E}_{\epsilon}(\mu)$ admit a unique
minimum? If so, do the minima of $\\{\mathcal{E}_{\epsilon}\\}_{\epsilon>0}$
converge to $\mu^{*}$ as $\epsilon\to 0$? In order to answer affirmatively to
these questions, we will need the associated kernel
$k_{\epsilon}(x,y):=\phi_{\epsilon}(x-y)$ to satisfy the following property.

######  Definition 3.4 (i.s.p.d. kernel).

A symmetric lower semicontinuous (l.s.c.) kernel $K$ on $X\times X$ is
_integrally strictly positive definite_ (i.s.p.d.) if for every finite signed
Borel measure $\nu$ on $X$, the energy $\mathcal{E}_{K}(\nu):=\iint
K(x,y)\mathop{}\\!\mathrm{d}(\nu\times\nu)(x,y)$ is well-defined (i.e., the
integrals over the negative and positive parts of $\nu$ are not both
infinite), and $\mathcal{E}_{K}(\nu)\geq 0$ where the equality holds only if
$\nu=0$ on all Borel sets of $X$.

For the mollifiers we consider, the associated kernel
$k_{\epsilon}(x,y)=\phi_{\epsilon}(x-y)$ is i.s.p.d. on any compact set
(Pronzato & Zhigljavsky (2021, Example 1.2)).

######  Proposition 3.5.

Suppose $X$ is compact and $k_{\epsilon}(x,y):=\phi_{\epsilon}(x-y)$ is
i.s.p.d. on $X$. Then,

  1. (a)

The kernel $W_{\epsilon}(x,y)$ defined in (1) also i.s.p.d on $X$.

  2. (b)

The functional $\mathcal{E}_{\epsilon}$ is strictly convex on
$\mathcal{M}_{\textup{sign}}(X)$ and attains a unique minimum on
$\mathcal{P}(X)$.

We next show the convergence of minima of $\\{\mathcal{E}_{\epsilon}\\}$ to
$\mu^{*}$ as a consequence of $\Gamma$-convergence of the sequence
$\\{\mathcal{E}_{\epsilon}\\}$.

######  Theorem 3.6.

Suppose $k_{\epsilon}(x,y):=\phi_{\epsilon}(x-y)$ is i.s.p.d. for every
$\epsilon>0$. Then we have $\Gamma$-convergence (Definition A.12)
$\mathcal{E}_{\epsilon}\stackrel{{\scriptstyle\Gamma}}{{\to}}\chi^{2}(\cdot\;\|\;\mu^{*})+1$
as $\epsilon\to 0$. In particular, if $X$ is compact, if we denote
$\mu_{\epsilon}^{*}:=\operatorname*{arg\,min}_{\mu\in\mathcal{P}(X)}\mathcal{E}_{\epsilon}(\mu)$,
then $\mu^{*}_{\epsilon}\to\mu^{*}$ weakly and $\lim_{\epsilon\to
0}\mathcal{E}_{\epsilon}(\mu_{\epsilon}^{*})=1$.

We prove Theorem 3.6 in Section A.3 using Fourier transforms and Bochner’s
theorem. Theorem 3.6 provides basis for minimizing $\mathcal{E}_{\epsilon}$
for a small $\epsilon$, since its unique minimum will be a good approximation
of $\mu^{*}$.

###  3.5 Differential calculus of $\mathcal{E}_{\epsilon}$ in
$\mathcal{P}_{2}(\mathbf{R}^{n})$

We next study the gradient flow of $\mathcal{E}_{\epsilon}$ in Wasserstein
space $\mathcal{P}_{2}(\mathbf{R}^{n})$ for $X=\mathbf{R}^{n}$. Understanding
the gradient flow of a functional often provides insights on 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 Section A.4.

Let $\mathcal{F}:\mathcal{P}_{2}(\mathbf{R}^{n})\to\mathbf{R}$ be a
functional. The _Wasserstein gradient flow_ of $\mathcal{F}$ (Ambrosio et al.,
2005, Definition 11.1.1) is defined as the solution $\\{\mu_{t}\\}_{t\geq 0}$
of the PDE: $\frac{\partial\mu_{t}}{\partial
t}=\divergence(\mu_{t}\bm{w}_{\mathcal{F},\mu_{t}})$ where
$\bm{w}_{\mathcal{F},\mu}\in L^{2}(\mu;\mathbf{R}^{n})$ is a _Frechét
subdifferential_ of $\mathcal{F}$ at $\mu$. Intuitively, gradient flows
capture the evolution of the variable being optimized if we were to do
gradient descent111The more precise notion is of minimizing movements
(Ambrosio et al., 2005, Chapter 2). on $\mathcal{F}$ when the step sizes go to
$0$. We next show that the gradient flow of $\mathcal{E}_{\epsilon}$ agrees
with that of $\chi^{2}(\cdot\;\|\;\mu^{*})$ in the sense that their
subdifferentials coincide as $\epsilon\to 0$.

######  Proposition 3.7.

Assume $\mu\in\mathcal{P}_{2}(\mathbf{R}^{n})$ has density $q\in
C_{c}^{1}(\mathbf{R}^{n})$222We assume that $\mu$ has compact support because
$W_{\epsilon}$ can be unbounded and cause integrability issues due to the
presence of $(p(x)p(y))^{-1/2}$ term.. Then any strong Frechét subdifferential
(19) of $\mathcal{E}_{\epsilon}$ at $\mu$ takes the form

$\displaystyle\bm{w}_{\epsilon,\mu}(x)=2\nabla\left(p(x)^{-1/2}(\phi_{\epsilon}*q/\sqrt{p})(x)\right),\text{
for }\mu\text{-a.e. }x\in\mathbf{R}^{n}.$ (4)

On the other hand, any strong Frechét subdifferential of
$\chi^{2}(\cdot\;\|\;\mu^{*})$ at $\mu$ takes the form

$\displaystyle\bm{w}_{\chi^{2},\mu}(x)=2\nabla\frac{q(x)}{p(x)},\text{ for
}\mu\text{-a.e. }x\in\mathbf{R}^{n}.$ (5)

Moreover, if for sufficiently small $\epsilon$, $\phi_{\epsilon}$ has compact
support independent of $\epsilon$, then

$\displaystyle\lim_{\epsilon\to
0}\bm{w}_{\epsilon,\mu}(x)=\bm{w}_{\chi^{2},\mu}(x),\text{ for }\mu\text{-a.e.
}x\in\mathbf{R}^{n}.$ (6)

While simulating $\chi^{2}$ divergence gradient flow is often intractable
(Trillos & Sanz-Alonso, 2020, Section 3.3), our MIE admits a straightforward
practical algorithm (Section 4). Chewi et al. (2020) show that SVGD can be
viewed as a kernelized gradient flow of $\chi^{2}$ divergence. However it is
unclear how to choose the kernel to recover precisely the
$\chi^{2}$-divergence gradient flow. On the other hand, they show the
exponential decay of KL divergence of $\chi^{2}$-divergence gradient flow if
the target measure satisfies a Poincaré inequality.

We next show that $\mathcal{E}_{\epsilon}$ is _displacement convex_ at
$\mu^{*}\in\mathcal{P}_{2}(\mathbf{R}^{n})$ as $\epsilon\to 0$, obtaining a
result similar to Korba et al. (2021, Corollary 4). This hints that gradient
flow initialized near $\mu^{*}$ will have fast convergence.

######  Proposition 3.8.

Suppose $p\in C^{2}_{c}(\mathbf{R}^{n})$. Suppose that for sufficiently small
$\epsilon$, $\phi_{\epsilon}$ has compact support independent of $\epsilon$.
Assume $k_{\epsilon}(x,y):=\phi_{\epsilon}(x-y)$ is i.s.p.d. Let $\bm{\xi}\in
C_{c}^{\infty}(\mathbf{R}^{n};\mathbf{R}^{n})$ and
$\mu_{t}:=(\bm{I}+t\bm{\xi})_{\\#}\mu^{*}$. Then $\lim_{\epsilon\to
0}\derivative[2]{t}\bigg{|}_{t=0}\mathcal{E}_{\epsilon}(\mu_{t})\geq 0.$

##  4 A practical sampling algorithm

We now present a first-order particle-based algorithm for constrained and
unconstrained sampling by minimizing a discrete version of (2). Substituting
an empirical distribution for $\mu$ in (2) gives the discrete mollified
interaction energy

$\displaystyle
E_{\epsilon}(\omega_{N}):=\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\phi_{\epsilon}(x_{i}-x_{j})(p(x_{i})p(x_{j}))^{-1/2}$
(7)

for $N$ particles $\omega_{N}=\\{x_{1},\ldots,x_{N}\\}$.

Denote $\omega^{*}_{N,\epsilon}\in\operatorname*{arg\,min}_{\omega_{N}\subset
X}E_{\epsilon}(\omega_{N})$ and
$\mu^{*}_{\epsilon}=\operatorname*{arg\,min}_{\mu\in\mathcal{P}_{2}(X)}\mathcal{E}_{\epsilon}(\mu)$,
If $X$ is compact, by Borodachov et al. (2019, Corollary 4.2.9), we have weak
convergence $\delta_{\omega^{*}_{N,\epsilon}}\to\mu^{*}_{\epsilon}$ as
$N\to\infty$. If in addition $\phi_{\epsilon}$ is i.s.p.d., then by Theorem
3.6, we have weak convergence $\mu^{*}_{\epsilon}\to\mu^{*}$. This shows that
minimizing (7) with a large $N$ and a small $\epsilon$ will result in an
empirical distribution of particles that approximates $\mu^{*}$. Our sampling
method, mollified interaction energy descent (MIED), is simply running
gradient descent on (7) over particles. 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 $\phi_{\epsilon}$ can also be hard to compute (e.g., for the
Riesz family). While these normalization constants do not affect the minima of
in (7), they can still affect gradient step sizes during optimization.
Moreover, $\phi_{\epsilon}$ can be very large when $\epsilon$ is small, and in
many Bayesian applications $p$ can be tiny and only $\log p$ is numerically
significant. To address these issues, we optimize the logarithm of (7) using
the log-sum-exp trick (Blanchard et al., 2021) to improve numerical stability
and to get rid of the arbitrary scaling of the normalizing constants:

$\displaystyle\log
E_{\epsilon}(\omega_{N}):=\log\sum_{i=1}^{N}\sum_{j=1}^{N}\exp(\log\phi_{\epsilon}(x_{i}-x_{j})-\frac{1}{2}(\log
p(x_{i})+\log p(x_{j})))-2\log N.$ (8)

Special treatment of the diagonal terms. Since $\phi_{\epsilon}$ goes to the
Dirac delta as $\epsilon\to 0$, the discretization of (2) on a neighborhood of
the diagonal $\\{(x,y):x=y\\}$ needs to be handled with extra care. In (7) the
diagonal appears as $\sum_{i=1}^{N}\phi_{\epsilon}(0)p(x_{i})^{-1}$ which can
dominant the summation when $\epsilon$ is small and then $\phi_{\epsilon}(0)$
becomes too large. For mollifiers we consider, we use a different diagonal
term $\sum_{i=1}^{N}\phi_{\epsilon}(h_{i}/\kappa_{n})p(x_{i})^{-1}$ where
$h_{i}:=\min_{j\neq i}\norm{x_{i}-x_{j}}$ and $\kappa_{n}\geq 1$ is a constant
depending only on the dimension $n$. Since
$0\leq\phi_{\epsilon}(h_{i}/\kappa)\leq\phi_{\epsilon}(0)$, the energy
obtained will be bounded between (7) and the version of (7) without the
diagonal terms. Hence by the proof of Theorem 4.2.2 of Borodachov et al.
(2019), the discrete minimizers of $E_{\epsilon}$ still converge to
$\mu^{*}_{\epsilon}$ as $N\to\infty$ and $\epsilon\to 0$. Empirically we found
the choice $\kappa_{n}=(1.3n)^{1/n}$ work well for the Riesz family of
mollifiers and we use the Riesz family primarily for our experiments.

Constraint handling. If $X\neq\mathbf{R}^{n}$, we need to minimize (8) subject
to constraints $\omega_{N}=\\{x_{1},\ldots,x_{N}\\}\subset X$. Since the
constraint is the same for each particle $x_{i}$, we want our algorithm to
remain parallelizable across particles.

We consider two types of constraints.

  1. (a)

There exists a differentiable map $f:\mathbf{R}^{n}\to X$ that so that
$\mathcal{L}_{n}(X\setminus f(\mathbf{R}^{n}))=0$.

  2. (b)

The set $X$ is given by $\\{x\in\mathbf{R}^{n}:g(x)\leq 0\\}$ for a
differentiable $g:\mathbf{R}^{n}\to\mathbf{R}$.

For (a), we reduce the problem to unconstrained optimization in
$\mathbf{R}^{n}$ using objective $\log E_{\epsilon}(f(\omega_{N}))$ with
$\omega_{N}=\\{x_{1},\ldots,x_{N}\\}\subset\mathbf{R}^{n}$ and
$f(\omega_{N}):=\\{f(x_{1}),\ldots,f(x_{N})\\}$. For (b), we apply the dynamic
barrier method by Gong & Liu (2021) to particles in parallel. In our
implementation, for both types of constraints, we use Adam optimizer (Kingma &
Ba, 2014) that modulates the update directions for each particle using history
of the moments of the gradients.

##  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 $\\{\phi^{s}_{\epsilon}\\}$ with $s=n+10^{-4}$ and
$\epsilon=10^{-8}$: we found minimizing the MIE with such mollifiers results
in well-separated particles so that we can take $\epsilon$ to be very small as
our theory recommends. This is not the case for the Gaussian or the Laplace
family as setting $\epsilon$ too small can cause numerical issues even when
particles are well-separated.

###  5.1 Unconstrained sampling

Gaussians in varying dimensions. We first compare MIED with SVGD (Liu & Wang,
2016) and KSDD (Korba et al., 2021) on Gaussians of varying dimensions and
with different number of particles. In Figure 1, we see that as the number of
dimensions increases, MIED results in best samples in terms of $W_{2}$
distance (Wasserstein-2 distance computed using linear programming) with
respect to $10^{4}$ 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 5) so the $W_{2}$ distance is lower, but it is biased
since $\epsilon>0$ so for energy distance SVGD performs better. More details
can be found in Section C.1.

Figure 1: Gaussian experiments results in varying dimensions and different
numbers of particles. For each method we plot the metric ($W_{2}$ distance or
energy distance) versus the number of particles, averaged over 10 trials
(shaded region indicates standard deviation).

Product of two Student’s $t$-distributions. Next, we consider a 2-dimensional
distribution constructed as the product of two independent $t$-distributions
with degree of freedom $\nu=2$ composed with a linear transform. Unlike
Gaussians, Student’s $t$-distributions have heavy tails. On the left of Figure
2, we 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 SVGD captures the distribution in $[-a,a]^{2}$
better for $a\leq 3$, MIED yields lower metrics for a bigger $a$.

Figure 2:  Left: visualization of samples from each method for the product of
two Student’s $t$-distribution (composed with a linear transform). Right:
metrics of each method when restricting to $[-a,a]^{2}$. As $t$-distributions
have heavy tail, if we draw i.i.d. 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]^{2}$ before computing the metrics for
$a\in[2,10]$.

Bayesian logistic regression. We compare MIED with SVGD and KSDD for the
Bayesian logistic regression setting in (Liu & Wang, 2016; Korba et al.,
2021). We further include a naïve baseline, independent particle descent
(IPD), which simply runs gradient descent independently on each particle to
maximize the posterior probability. In addition to using test accuracy as the
metric, we include $W_{2}$ distance and energy distance between the samples
from each method and $10^{4}$ samples from NUTS (Hoffman et al., 2014) after
sufficient burn-in steps. We summarize the results in Table 1. All methods,
including the naïve baseline IPD, are comparable in terms of test accuracy. In
other words, accuracy is not a good metric for comparing the quality of
posterior samples. Bayesian inference is typically preferred 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 $W_{2}$
distance and all methods (except IPD) are comparable in terms of the energy
distance.

Dataset | NUTS | IPD | SVGD | KSDD | MIED
---|---|---|---|---|---
banana ($d=3$) | 0.55 | -6.14/-3.58/0.55 | -7.81/-5.24/0.55 |  -8.24/-5.76/0.55 | -7.37/-5.06/0.55
breast cancer ($d=10$) | 0.64 | -1.51/-1.03/0.60 | -1.62/-2.06/0.60 | -1.71/-2.23/0.60 |  -1.99/-2.18/0.60
diabetis ($d=9$) | 0.78 | -2.18/-1.55/0.77 | -3.09/-3.42/0.77 | -2.91/-3.90/0.77 |  -3.11/-3.13/0.77
flare solar ($d=10$) | 0.59 | 3.30/2.65/0.48 | 6.91/4.09/0.52 |  1.77/-0.08/0.55 | 7.09/4.25/0.48
german ($d=21$) | 0.65 | -1.80/-1.25/0.65 | -1.89/-2.63/0.64 | -1.27/-2.83/0.65 |  -1.96/-2.80/0.65
heart ($d=14$) | 0.87 | -0.40/-0.56/0.87 | -0.41/-1.50/0.87 | -0.10/-1.76/0.87 |  -0.92/-1.67/0.87
image ($d=19$) | 0.82 | 6.53/4.31/0.83 | 7.17/4.01/0.83 | 2.16/-0.50/0.82 |  1.14/-1.88/0.82
ringnorm ($d=21$) | 0.77 | -3.82/-2.45/0.77 |  -4.11/-5.98/0.77 | 1.07/-2.21/0.76 | -4.03/-5.70/0.77
splice ($d=61$) | 0.85 |  -1.47/-1.18/0.85 | -1.22/-2.65/0.85 | 2.04/-0.05/0.84 | 1.45/0.70/0.82
thyroid ($d=6$) | 0.84 | 1.95/0.53/0.84 | 1.17/-0.00/0.84 | 2.42/1.57/0.74 |  0.84/-0.37/0.84
titanic ($d=4$) | 0.40 |  -1.59/-0.16/0.40 | -0.46/-0.31/0.40 | -0.63/-0.39/0.40 | -1.00/-0.45/0.40
twonorm ($d=21$) | 0.97 | -1.21/-1.13/0.97 | -1.32/-2.78/0.97 | 1.55/-0.62/0.97 |  -1.44/-3.21/0.97
waveform ($d=22$) | 0.77 | -2.67/-1.87/0.78 | -2.98/-5.23/0.78 | -2.60/-4.18/0.77 |  -3.09/-3.17/0.78
Table 1: 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 columnes are logarithmic $W_{2}$ distance, logarithmic energy
distance, and test accuracy. Bold indicates the best numbers. We use
$80\%/20\%$ training/test split. All methods are run with identical
initialization and learning rate $0.01$. Results are reported after $10^{4}$
iterations.

###  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 various choices of
mollifiers—we always choose the smallest $\epsilon$ while the optimization
remains stable. We compare our method with SVMD and MSVGD proposed for
constrained sampling by Shi et al. (2021) using mirror map
$\phi(\theta)=\sum_{i=1}^{n}\left(\log\frac{1}{1-\theta_{i}}+\log\frac{1}{1+\theta_{i}}\right)$
as in Ahn & Chewi (2021). 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}$ for all methods.

Figure 3: Convergence of metrics and visualization of samples for each method
for uniform sampling from a 2D box.

The left plot of Figure 3 shows that quantitatively our method achieves much
lower energy distance and $W_{2}$ distance (compared to a $5000$ i.i.d.
samples uniformly drawn in $[-1,1]^{2}$.). On the right of Figure 3, we
visualize the samples from each method. We see that samples from SVGD collapse
and samples from KSDD are non-uniform (we choose kernel size $0.1$ for KSDD
that produces the best result; otherwise it will collapse like SVGD). SVMD and
MSVGD fail to draw samples near the boundary because they use mapping
$\nabla\phi^{*}(\eta)=\nicefrac{{(\sqrt{1+\eta^{2}}-1)}}{{\eta}}$ (coordinate-
wise arithmetic) to push unconstrained samples to $[-1,1]^{2}$, yet
coordinates of $\eta$ need to go to $\infty$ to make $\nabla\phi^{*}(\eta)$
land close to the boundary. We verify this phenomenon by using
$\nabla\phi^{*}(\eta)$ as the reparametrization map in MIED: indeed as shown
the the bottom right of Figure 3, with such reparameterization MIED also fails
to capture samples near the boundary. When using $\tanh$ as the
reparameterization map, MIED produces decent results for all three families of
mollifiers (Riesz, Gaussian, Laplace). This highlight the flexibility of our
method with constraint handling, whereas for SVMD and MSVGD, the mirror map
has to chosen carefully to be the gradient of a strongly convex function while
capturing the constraints.

In Section C.2, we conduct a 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 the mirror map is available, SVMD
and MSVGD can obtain better performance compared to MIED.

Fairness Bayesian neural networks. We train fair bayesian neural networks to
predict whether the annual income of a person is at least $\$50,000$ with
gender as the protected attribute following the same setup as in Liu et al.
(2021) where we are given a dataset
$\mathcal{D}=\\{x^{(i)},y^{(i)},z^{(i)}\\}_{i=1}^{|\mathcal{D}|}$ consisting
of feature vectors $x^{(i)}$, labels $y^{i}$ (whether the income is
$\geq\$50,000$), and genders $z^{(i)}$ (protected attribute). The target
density is taken to be the posterior of logistic regression with a two-layer
Bayesian neural network with weights $\theta$, and we put a standard Gaussian
prior on each entry of $\theta$ independently. Given $t>0$, the fairness
constraint is
$g(\theta)=(\textup{cov}_{(x,y,z)\sim\mathcal{D}}[z,\hat{y}(x;\theta)])^{2}-t\leq
0$. We use $80\%/20\%$ training/test split as in Liu et al. (2021). 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
$t\in\\{10^{-5},10^{-4},0.0001,0.001,0.002,0.005,0.01\\}$. We use the source
code provided by Liu et al. (2021) with default hyperparameters.333The source
code provided by the authors of Liu et al. (2021) does not implement the
calculation of $g(\theta)$ faithfully as written in the formula, so we
modified it. All methods use 2000 iterations for training, and for our method
we use learning rate $0.001$. One of their four methods (Control+SVGD) got
stuck at initialization (with accuracy around $0.75$), so we omit its result
from the plot. Our method recovers a much larger Pareto front compared to the
alternatives. On the right of Figure 4, we visualize the curves of the energy
and the covariance versus the number of training iterations: as expected we
see a smaller $t$ results in bigger MIE (lower log-likelihood) and smaller
covariance between the prediction and the protected attribute.

Figure 4: Left: trade-off curves of demographic parity versus test accuracy
for MIED and methods from Liu et al. (2021). Right: MIEs and
$(\textup{cov}_{(x,y,z)\sim\mathcal{D}}[z,\hat{y}(x;\theta)])^{2}$ (measured
on the training data) versus the number of training iterations, for various
$t$.

##  6 Conclusion

We present a new sampling method by minimizing MIEs discretized as particles
for unconstrained and constrained sampling. This is motivated by the insight
that MIEs converge to $\chi^{2}$ divergence with respect to the target measure
as the mollifier parameter goes to 0. The proposed method achieves promising
results on the sampling problems we consider.

Below we highlight three limitations. First, as discussed in Remark A.9, our
theory only applies when the domain is full-dimensional or flat. Extending our
theory to handle cases where the domain is an arbitrary $d$-rectifiable set is
an important next step as it allows handling of more complicated constraints
such as nonlinear equality constraints. Secondly, when $\epsilon>0$, the
minimizer of MIE can be different from the target mesaure. Finding a way to
debias MIE (e.g., like how Sinkhorn distances are debiased (Feydy et al.,
2019)) is an interesting direction. Lastly, the connection between the
minimizers of the discretized MIE (7) and those of the continuous MIE (2) is
only established in the limit as $N\to\infty$. We hope to investigate how well
the empirical distribution of particles minimizing (7) approximates the target
measure when $N$ is finite.

#### Reproducibility statement.

We provide self-contained proofs in Appendix A for the theoretical results
stated in the main text.

We have included the complete source code in the supplementary necessary to
reproduce all our experiments.

## References

  * Ahn & Chewi (2021) Kwangjun Ahn and Sinho Chewi.  Efficient constrained sampling via the mirror-langevin algorithm.  _Advances in Neural Information Processing Systems_ , 34:28405–28418, 2021.
  * Ambrosio et al. (2005) Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré.  _Gradient flows: in metric spaces and in the space of probability measures_.  Springer Science & Business Media, 2005.
  * Balasubramanian et al. (2022) Krishna Balasubramanian, Sinho Chewi, Murat A Erdogdu, Adil Salim, and Shunshi Zhang.  Towards a theory of non-log-concave sampling: first-order stationarity guarantees for langevin monte carlo.  In _Conference on Learning Theory_ , pp. 2896–2923. PMLR, 2022\.
  * Blanchard et al. (2021) Pierre Blanchard, Desmond J Higham, and Nicholas J Higham.  Accurately computing the log-sum-exp and softmax functions.  _IMA Journal of Numerical Analysis_ , 41(4):2311–2330, 2021.
  * Blei et al. (2017) David M Blei, Alp Kucukelbir, and Jon D McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Borodachov et al. (2019) Sergiy V Borodachov, Douglas P Hardin, and Edward B Saff.  _Discrete energy on rectifiable sets_.  Springer, 2019.
  * Brooks et al. (2011) Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng.  _Handbook of markov chain monte carlo_.  CRC press, 2011.
  * Brubaker et al. (2012) Marcus Brubaker, Mathieu Salzmann, and Raquel Urtasun.  A family of mcmc methods on implicitly defined manifolds.  In _Artificial intelligence and statistics_ , pp. 161–172. PMLR, 2012.
  * Byrne & Girolami (2013) Simon Byrne and Mark Girolami.  Geodesic monte carlo on embedded manifolds.  _Scandinavian Journal of Statistics_ , 40(4):825–845, 2013.
  * Carrillo et al. (2014) 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.
  * Chewi et al. (2020) Sinho Chewi, Thibaut Le Gouic, Chen Lu, Tyler Maunu, and Philippe Rigollet.  Svgd as a kernelized wasserstein gradient flow of the chi-squared divergence.  _Advances in Neural Information Processing Systems_ , 33:2098–2109, 2020.
  * Chizat (2022) Lenaic Chizat.  Sparse optimization on measures with over-parameterized gradient descent.  _Mathematical Programming_ , 194(1):487–532, 2022\.
  * Craig et al. (2022) Katy Craig, Karthik Elamvazhuthi, Matt Haberland, and Olga Turanova.  A blob method method for inhomogeneous diffusion with applications to multi-agent control and sampling.  _arXiv preprint arXiv:2202.12927_ , 2022.
  * Durmus et al. (2019) Alain Durmus, Szymon Majewski, and Błażej Miasojedow.  Analysis of langevin monte carlo via convex optimization.  _The Journal of Machine Learning Research_ , 20(1):2666–2711, 2019.
  * Feydy et al. (2019) Jean Feydy, Thibault Séjourné, François-Xavier 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 Statistics_ , pp. 2681–2690. PMLR, 2019.
  * Friedrichs (1944) Kurt Otto Friedrichs.  The identity of weak and strong extensions of differential operators.  _Transactions of the American Mathematical Society_ , 55(1):132–151, 1944.
  * Gong & Liu (2021) Chengyue Gong and Xingchao Liu.  Bi-objective trade-off with dynamic barrier gradient descent.  _NeurIPS 2021_ , 2021.
  * Gorham & Mackey (2017) Jackson Gorham and Lester Mackey.  Measuring sample quality with kernels.  In _International Conference on Machine Learning_ , pp. 1292–1301. PMLR, 2017.
  * Hoffman et al. (2014) Matthew D Hoffman, Andrew Gelman, et al.  The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo.  _J. Mach. Learn. Res._ , 15(1):1593–1623, 2014\.
  * Hörmander (2015) Lars Hörmander.  _The analysis of linear partial differential operators I: Distribution theory and Fourier analysis_.  Springer, 2015.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Klenke (2013) Achim Klenke.  _Probability theory: a comprehensive course_.  Springer Science & Business Media, 2013.
  * Korba et al. (2020) Anna Korba, Adil Salim, Michael Arbel, Giulia Luise, and Arthur Gretton.  A non-asymptotic analysis for stein variational gradient descent.  _Advances in Neural Information Processing Systems_ , 33:4672–4682, 2020.
  * Korba et al. (2021) Anna Korba, Pierre-Cyril Aubin-Frankowski, Szymon Majewski, and Pierre Ablin.  Kernel stein discrepancy descent.  In _International Conference on Machine Learning_ , pp. 5719–5730. PMLR, 2021.
  * Kühn (2016) Franziska Kühn.  Probability and heat kernel estimates for lévy (-type) processes.  2016\.
  * Liu & Zhu (2018) 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.
  * Liu (2017) Qiang Liu.  Stein variational gradient descent as gradient flow.  _Advances in neural information processing systems_ , 30, 2017.
  * Liu & Wang (2016) Qiang Liu and Dilin Wang.  Stein variational gradient descent: A general purpose bayesian inference algorithm.  _Advances in neural information processing systems_ , 29, 2016.
  * Liu et al. (2021) Xingchao Liu, Xin Tong, and Qiang Liu.  Sampling with trusthworthy constraints: A variational gradient framework.  _Advances in Neural Information Processing Systems_ , 34:23557–23568, 2021.
  * Pronzato & Zhigljavsky (2021) Luc Pronzato and Anatoly Zhigljavsky.  Minimum-energy measures for singular kernels.  _Journal of Computational and Applied Mathematics_ , 382:113089, 2021.
  * Salim et al. (2022) Adil Salim, Lukang Sun, and Peter Richtarik.  A convergence theory for svgd in the population limit under talagrand’s inequality t1.  In _International Conference on Machine Learning_ , pp. 19139–19152. PMLR, 2022.
  * Serfaty (2020) Sylvia Serfaty.  Mean field limit for coulomb-type flows.  _Duke Mathematical Journal_ , 169(15):2887–2935, 2020.
  * Shi et al. (2021) Jiaxin Shi, Chang Liu, and Lester Mackey.  Sampling with mirrored stein operators.  _arXiv preprint arXiv:2106.12506_ , 2021.
  * Székely & Rizzo (2013) 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–1272, 2013.
  * Trillos & Sanz-Alonso (2020) Nicolas Garcia Trillos and Daniel Sanz-Alonso.  The bayesian update: variational formulations and gradient flows.  _Bayesian Analysis_ , 15(1):29–56, 2020.
  * Varadhan (2001) SR Srinivasa Varadhan.  _Probability theory_.  Number 7. American Mathematical Soc., 2001.
  * Wibisono (2018) Andre Wibisono.  Sampling as optimization in the space of measures: The langevin dynamics as a composite optimization problem.  In _Conference on Learning Theory_ , pp. 2093–3027. PMLR, 2018\.
  * Xu (2021) Wenkai Xu.  Generalised kernel stein discrepancy (gksd): A unifying approach for non-parametric goodness-of-fit testing.  _arXiv preprint arXiv:2106.12105_ , 2021.
  * Zhang et al. (2020) 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\.

##  Appendix A Detailed Analysis

###  A.1 Preliminaries on mollifiers

######  Proposition A.1.

For $s>n$, the $s$-Riesz family of mollifiers, defined as
$\phi^{s}_{\epsilon}(x):=(\norm{x}^{2}_{2}+\epsilon^{2})^{-s/2}/Z^{s}_{\epsilon}$,
satisfies Definition 3.1.

In order to prove Proposition A.1, we first prove a lemma.

######  Lemma A.2.

Let $b\geq 0$. Then for any $\epsilon>0$,

$\displaystyle
Z_{\epsilon}^{s}\int_{B_{\epsilon}(0)}\phi^{s}_{\epsilon}(y)\norm{y}_{2}^{b}\mathop{}\\!\mathrm{d}y=\mathcal{H}_{n-1}(S^{n-1})\left(\int_{0}^{1}\frac{t^{n+b-1}}{(t^{2}+1)^{s/2}}\right)\epsilon^{n+b-s},$
(9)

where $\mathcal{H}_{n-1}(S^{n-1})$ is the volume of the $(n-1)$-dimensional
sphere.

Furthermore, assuming $s>b+n$, then for any $\epsilon>0,\delta>0$,

$\displaystyle Z_{\epsilon}^{s}\int_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi^{s}_{\epsilon}(y)\norm{y}_{2}^{b}\mathop{}\\!\mathrm{d}y$
$\displaystyle\leq\mathcal{H}_{n-1}(S^{n-1})\frac{\delta^{n+b-s}}{s-(n+b)}.$
(10)

###### Proof.

For (9), we compute

$\displaystyle
Z_{\epsilon}^{s}\int_{B_{\epsilon}(0)}\phi^{s}_{\epsilon}(y)\norm{y}_{2}^{b}\mathop{}\\!\mathrm{d}y$
$\displaystyle=\int_{B_{\epsilon}(0)}\frac{\norm{y}_{2}^{b}}{(\norm{y}_{2}^{2}+\epsilon^{2})^{s/2}}=\mathcal{H}_{n-1}(S^{n-1})\int_{0}^{\epsilon}\frac{r^{n+b-1}}{(r^{2}+\epsilon^{2})^{s/2}}\mathop{}\\!\mathrm{d}r$
$\displaystyle=\mathcal{H}_{n-1}(S^{n-1})\int_{0}^{1}\frac{(\epsilon
t)^{n+b-1}}{(\epsilon^{2}t^{2}+\epsilon^{2})^{s/2}}\epsilon\mathop{}\\!\mathrm{d}t$
$\displaystyle=\mathcal{H}_{n-1}(S^{n-1})\left(\int_{0}^{1}\frac{t^{n+b-1}}{(t^{2}+1)^{s/2}}\mathop{}\\!\mathrm{d}t\right)\epsilon^{n+b-s},$

where we use substitution $r=\epsilon t$ on the second line.

If $s>b+n$, then

$\displaystyle Z_{\epsilon}^{s}\int_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi^{s}_{\epsilon}(0,y)\norm{y}_{2}^{b}\mathop{}\\!\mathrm{d}y$
$\displaystyle=\int_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\frac{\norm{y}_{2}^{b}}{(\norm{y}_{2}^{2}+\epsilon^{2})^{s/2}}$
$\displaystyle\leq\int_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\frac{\norm{y}_{2}^{b}}{\norm{y}_{2}^{s}}=\mathcal{H}_{n-1}(S^{n-1})\int_{\delta}^{\infty}r^{n-1+b-s}$
$\displaystyle=\mathcal{H}_{n-1}(S^{n-1})\frac{\delta^{n+b-s}}{s-(n+b)},$

where the last integral uses $s>b+n$. ∎

###### Proof of Proposition A.1.

It is clear that (a) of Proposition A.1 holds for $\phi_{\epsilon}^{s}(x)$.
Since $\phi_{\epsilon}^{s}(x)$ is bounded we have $\phi_{\epsilon}^{s}\in
L^{\infty}(\mathbf{R}^{n})$. For any $\delta>0$, we have, for
$\epsilon\leq\delta$,

$\displaystyle
Z_{\epsilon}^{s}\int_{B_{\delta}(0)}\phi^{s}_{\epsilon}(x)\mathop{}\\!\mathrm{d}x\geq
Z_{\epsilon}^{s}\int_{B_{\epsilon}(0)}\phi^{s}_{\epsilon}(x)\mathop{}\\!\mathrm{d}x=C\epsilon^{n-s},$

where we use (9) with $b=0$ and $C$ is a constant depending only on $n,s$. On
the other hand, using (10) with $b=0$,

$\displaystyle Z_{\epsilon}^{s}\int_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi^{s}_{\epsilon}(x)\mathop{}\\!\mathrm{d}x\leq
C^{\prime}\delta^{n-s},$

where $C^{\prime}$ depends only on $n,s$. With $\delta=\epsilon$, we see that
$Z_{\epsilon}^{s}\phi_{\epsilon}^{s}\in L^{1}(\mathbf{R}^{n})$ so (b) is
satisfied.

Since $\delta$ is fixed and $s>n$, we see that

$\lim_{\epsilon\to
0}\frac{\int_{B_{\delta}(0)}\phi_{\epsilon}^{s}(x)\mathop{}\\!\mathrm{d}x}{\int_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}^{s}(x)\mathop{}\\!\mathrm{d}x}\geq\lim_{\epsilon\to
0}\frac{C\epsilon^{n-s}}{C^{\prime}\delta^{n-s}}=\infty.$

Since $\int_{B_{\delta}(0)}\phi_{\epsilon}^{s}(x)\mathop{}\\!\mathrm{d}x$
$+\int_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}^{s}(x)\mathop{}\\!\mathrm{d}x=1$, we have shown
(c). ∎

In the rest of this section, we assume $\\{\phi_{\epsilon}\\}_{\epsilon>0}$ is
a family of mollifiers satisfying Definition 3.1.

######  Lemma A.3.

For any $p\in[1,\infty]$, for any $\delta>0$, $\phi_{\epsilon}\in
L^{p}(\mathbf{R}^{n})$ and $\lim_{\epsilon\to
0}\norm{\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}}_{p}=0$. holds.

###### Proof.

Assume $p\notin\\{1,\infty\\}$ since both cases are covered in the
assumptions. Then for any $\delta>0$, by Hölder’s inequality,

$\displaystyle\norm{\phi_{\epsilon}}_{p}^{p}=\norm{\phi_{\epsilon}\cdot\phi_{\epsilon}^{p-1}}_{1}\leq\norm{\phi_{\epsilon}}_{1}\norm{\phi_{\epsilon}^{p-1}}_{\infty}=\norm{\phi_{\epsilon}}_{1}\norm{\phi_{\epsilon}}_{\infty}^{p-1}<\infty.$

Similarly,

$\displaystyle\norm{\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}}_{p}^{p}$
$\displaystyle=\norm{\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}^{p}}_{1}=\norm{\phi_{\epsilon}\cdot\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}^{p-1}}_{1}$
$\displaystyle\leq\norm{\phi_{\epsilon}}_{1}\norm{\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}^{p-1}}_{\infty}=\norm{\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}}_{\infty}^{p-1}.$

Letting $\epsilon\to 0$ and applying (c) gives $\lim_{\epsilon\to
0}\norm{\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}}_{p}=0$. ∎

######  Proposition A.4.

Let $f\in L^{p}(\mathbf{R}^{n})$, $p\in[1,\infty]$. Then for every
$\epsilon>0$, the integral

$\int f(x-y)\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y$

exists, so that $(f*\phi_{\epsilon})(x)$ is finite. Moreover, if $f$ is
continous at $x$, then

$\displaystyle\lim_{\epsilon\to 0}(f*\phi_{\epsilon})(x)=f(x).$ (11)

###### Proof.

Observe that by Hölder’s inequality, for $q$ such that $1/p+1/q=1$ (allowing
infinity),

$\displaystyle\int\absolutevalue{f(x-y)\phi_{\epsilon}(y)}\mathop{}\\!\mathrm{d}y$
$\displaystyle\leq\norm{f}_{p}\norm{\phi_{\epsilon}}_{q}<\infty.$

Hence $(f*\phi_{\epsilon})(x)$ is integrable and finite.

For any $\epsilon>0$, note that

$\displaystyle\absolutevalue{(f*\phi_{\epsilon})(x)-f(x)}$
$\displaystyle=\absolutevalue{\int
f(x-y)\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y-f(x)}.$

Since $\int\phi_{\epsilon}(x)\mathop{}\\!\mathrm{d}x=1$, we have

$\displaystyle\absolutevalue{(f*\phi_{\epsilon})(x)-f(x)}$
$\displaystyle=\absolutevalue{\int\left(f(x-y)-f(x)\right)\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y}\leq\int\absolutevalue{f(x-y)-f(x)}\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y.$

Fix $t>0$. Continuity of $f$ at $x$ implies there exists $\delta>0$ such that
$\absolutevalue{f(x-y)-f(x)}<t$ for all $y\in B_{\delta}(0)$. Then

$\displaystyle\int_{B_{\delta}(0)}\absolutevalue{f(x-y)-f(x)}\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y$
$\displaystyle\leq
t\int_{B_{\delta}(0)}\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y\leq t.$

On the other hand, by Hölder’s inequality,

$\displaystyle\int_{\mathbf{R}^{n}\setminus
B_{\delta}(x)}\absolutevalue{f(x-y)-f(x)}\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y\leq\norm{\tau_{x}f-f}_{p}\norm{\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}}_{q}\leq
2\norm{f}_{p}\norm{\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}}_{q}.$

Hence

$\displaystyle\absolutevalue{(f*\phi_{\epsilon})(x)-f(x)}$ $\displaystyle\leq
t+2\norm{f}_{p}\norm{\mathbf{1}_{\mathbf{R}^{n}\setminus
B_{\delta}(0)}\phi_{\epsilon}}_{q}.$

By Lemma A.3, since $\norm{f}_{p}<\infty$, taking $\epsilon\to 0$ we get, for
any $t>0$,

$\displaystyle\limsup_{\epsilon\to
0}\absolutevalue{(f*\phi_{\epsilon})(x)-f(x)}$ $\displaystyle\leq t.$

Now let $t\to 0$ we get (11). ∎

######  Corollary A.5.

Let $f\in L^{p}(\mathbf{R}^{n})$, $p\in[1,\infty]$, be a continuous function.
If $\\{f_{\epsilon}\\}$ is a sequence of functions that converge uniformly to
$f$, then for every $x$,

$\displaystyle\lim_{\epsilon\to 0}(f_{\epsilon}*\phi_{\epsilon})(x)=f(x).$
(12)

###### Proof.

Note that

$\displaystyle\absolutevalue{(\phi_{\epsilon}*f_{\epsilon})(x)-f(x)}$
$\displaystyle\leq\absolutevalue{(\phi_{\epsilon}*f_{\epsilon})(x)-(\phi_{\epsilon}*f)(x)}+\absolutevalue{(\phi_{\epsilon}*f)(x)-f(x)}.$

Proposition A.4 shows the second term goes to $0$ as $\epsilon\to 0$. For the
first term, we have

$\displaystyle\absolutevalue{(\phi_{\epsilon}*f_{\epsilon})(x)-(\phi_{\epsilon}*f)(x)}$
$\displaystyle=\absolutevalue{\int\left(f_{\epsilon}(x-y)-f(x-y)\right)\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y}$
$\displaystyle\leq\sup_{x}\absolutevalue{f_{\epsilon}(x)-f(x)}\to 0$

by uniform convergence. ∎

######  Lemma A.6.

For $f\in L^{p}(\mathbf{R}^{n})$, $p\in[1,\infty)$, we have

$\displaystyle\lim_{y\to 0}\norm{\tau_{y}f-f}_{p}=0,$

where we use $\tau_{a}f$ to denote the translated function
$\tau_{a}f(x):=f(x-a)$.

###### Proof.

Fix $\epsilon>0$. It is a standard fact that $C_{c}(\mathbf{R}^{n})$ is dense
in $L^{p}(\mathbf{R}^{n})$. Hence there exists $g\in C_{c}(\mathbf{R}^{n})$
such that $\norm{f-g}_{p}<\epsilon$. Since $g$ is continuous with compact
support, it is uniform continuous. Then there exists $\delta>0$ with
$\absolutevalue{g(x-y)-g(x)}<\epsilon^{1/p}/\mathcal{L}_{n}(K)$ if $y\in
B_{\delta}(x)$. Hence for such $y$ we have
$\norm{\tau_{y}g-g}^{p}_{p}<\epsilon$ with $\mathcal{L}_{n}(K)<\infty$. Thus

$\displaystyle\norm{\tau_{y}f-f}_{p}$
$\displaystyle\leq\norm{\tau_{y}f-\tau_{y}g}_{p}+\norm{\tau_{y}g-g}_{p}+\norm{g-f}_{p}\leq
3\epsilon.$

∎

######  Proposition A.7.

Let $f\in L^{p}(\mathbf{R}^{n})$, $p\in(1,\infty)$. Then for every
$\epsilon>0$,

$\displaystyle\norm{f*\phi_{\epsilon}}_{p}\leq 3\norm{f}_{p}.$

In particular, $f*\phi_{\epsilon}\in L^{p}(\mathbf{R}^{n})$. Moreover,

$\displaystyle\lim_{\epsilon\to 0}\norm{f*\phi_{\epsilon}-f}_{p}=0.$

###### Proof.

For every $\epsilon>0$, we have

$\displaystyle\norm{f*\phi_{\epsilon}-f}^{p}_{p}$
$\displaystyle=\int\absolutevalue{\int(f(x-y)-f(x))\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y}^{p}\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq\int\int\absolutevalue{f(x-y)-f(x)}^{p}\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y\mathop{}\\!\mathrm{d}x$
$\displaystyle=\int\left(\int\absolutevalue{f(x-y)-f(x)}^{p}\mathop{}\\!\mathrm{d}x\right)\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y$
$\displaystyle=\int\norm{\tau_{y}f-f}_{p}^{p}\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y,$

where we use Jensen’s inequality with the observation that
$\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y$ is a probability measure and
Tonelli’s theorem to exchange the order of integration. To show the first
claim, note that,

$\displaystyle\norm{f*\phi_{\epsilon}-f}^{p}_{p}\leq\int\norm{\tau_{y}f-f}_{p}^{p}\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y$
$\displaystyle\leq\int(2\norm{f}_{p})^{p}\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y=(2\norm{f}_{p})^{p}.$

Hence
$\norm{f*\phi_{\epsilon}}_{p}\leq\norm{f}_{p}+\norm{f*\phi_{\epsilon}-f}_{p}\leq
3\norm{f}_{p}$.

For the second claim, by Lemma A.6, the function
$y\mapsto\norm{\tau_{y}f-f}_{p}^{p}$ is continuous at $y=0$ with limit $0$.
Hence by Proposition A.4, we are done by taking $\epsilon\to 0$. ∎

###  A.2 Convergence to $\chi^{2}$-divergence

We start by giving an alternative formula for $\chi^{2}$ divergence with
respect to the target measure.

######  Lemma A.8.

Define a functional $\mathcal{E}:\mathcal{P}(\mathbf{R}^{n})\to[0,\infty]$ by

$\displaystyle\begin{split}\mathcal{E}(\mu):=\left\\{\begin{array}[]{ll}\int_{X}\frac{q(x)^{2}}{p(x)}\mathop{}\\!\mathrm{d}x&\text{
if
}\mathop{}\\!\mathrm{d}\mu(x)=q(x)\mathop{}\\!\mathrm{d}\mathcal{L}_{n}(x)\text{
and $\mu\ll\mu^{*}$}\\\ \infty&\text{
otherwise.}\end{array}\right.\end{split}$ (13)

Then for every $\mu\in\mathcal{P}(\mathbf{R}^{n})$,

$\mathcal{E}(\mu)=\chi^{2}(\mu\;\|\;\mu^{*})+1.$

###### Proof.

If $\mu$ is not absolutely continuous with respect to $\mu^{*}$, then both
sides are infinity. Otherwise, $\mu(X)=1$ and $\mu$ has some density $q$. We
then compute

$\displaystyle\chi^{2}(\mu\;\|\;\mu^{*})+1$
$\displaystyle=\int\left(\frac{\mathop{}\\!\mathrm{d}\mu}{\mathop{}\\!\mathrm{d}\mu^{*}}(x)-1\right)^{2}\mathop{}\\!\mathrm{d}\mu^{*}(x)+1=\int_{X}\left(\frac{q(x)}{p(x)}-1\right)^{2}p(x)\mathop{}\\!\mathrm{d}x+1$
$\displaystyle=\int_{X}\frac{q(x)^{2}}{p(x)}\mathop{}\\!\mathrm{d}x-2\int_{X}q(x)\mathop{}\\!\mathrm{d}x+2=\int_{X}\frac{q(x)^{2}}{p(x)}\mathop{}\\!\mathrm{d}x=\mathcal{E}(\mu).$

∎

###### Proof of Theorem 3.3.

Since $\chi^{2}(\mu\;\|\;\mu^{*})<\infty$, we know $\mu\ll\mu^{*}$, so
$\mu(X)=1$. We let $q(x)$ be the density of $\mu$ which satisfies $q(x)=0$ for
$x\in\mathbf{R}^{n}\setminus X$. We compute using Tonelli’s theorem (since our
integrand is positive):

$\displaystyle\mathcal{E}_{\epsilon}(\mu)$
$\displaystyle=\iint\phi_{\epsilon}(x-y)(p(x)p(y))^{-1/2}q(x)q(y)\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y$
$\displaystyle=\int\left(\phi_{\epsilon}*\frac{q}{\sqrt{p}}\right)(x)\frac{q}{\sqrt{p}}(x)\mathop{}\\!\mathrm{d}x=\int(\phi_{\epsilon}*g)(x)g(x)\mathop{}\\!\mathrm{d}x,$

where we define $g(x):=q(x)/\sqrt{p(x)}$ on $X$ and $g(x):=0$ for $x\notin X$.
Moreover, by Lemma A.8,

$\displaystyle\chi^{2}(\mu\;\|\;\mu^{*})+1=\int
g(x)^{2}\mathop{}\\!\mathrm{d}x.$

The assumption $\chi^{2}(\mu\;\|\;\mu^{*})<\infty$ then implies $g\in
L^{2}(\mathbf{R}^{n})$. By Proposition A.7, we have $\phi_{\epsilon}*g\in
L^{2}(\mathbf{R}^{n})$. Next notice that

$\displaystyle\absolutevalue{\mathcal{E}_{\epsilon}(\mu)-\left(\chi^{2}(\mu\;\|\;\mu^{*})+1\right)}$
$\displaystyle\leq\int\absolutevalue{(\phi_{\epsilon}*g)(x)-g(x)}g(x)\mathop{}\\!\mathrm{d}x\leq\norm{\phi_{\epsilon}*g-g}_{2}\norm{g}_{2}<\infty.$

This shows the first claim. Finally, the last expression goes to $0$ as
$\epsilon\to 0$ since $\norm{\phi_{\epsilon}*g-g}_{2}\to 0$ by Proposition
A.7. ∎

######  Remark A.9.

It is possible to prove versions of Theorem 3.3 when $X$ has Hausdorff
dimension $d<n$: in such cases $\chi^{2}$-divergence 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 $\mathbf{R}^{n}$, 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$, a similar calculation
yields
$\mathcal{E}_{\epsilon}(\mu)=\int_{X}\left(\int_{X}\phi_{\epsilon}(x-y)\left(\frac{q}{\sqrt{p}}\right)(y)\mathop{}\\!\mathrm{d}\mathcal{H}_{d}(y)\right)\left(\frac{q}{\sqrt{p}}\right)(x)\mathop{}\\!\mathrm{d}\mathcal{H}_{d}(x).$
For a similar argument to go through, we will need
$\int_{X}\phi_{\epsilon}(x-y)\mathop{}\\!\mathrm{d}\mathcal{H}_{d}(y)$ to be
independent of $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 a
version of Theorem 3.3 for $X$ with Hausdorff dimension $d<n$ is an
interesting future direction.

###  A.3 Convexity and $\Gamma$-convergence

We start by recalling a few definitions regarding functionals in
$\mathcal{P}(\mathbf{R}^{n})$.

######  Definition A.10.

We say a functional
$\mathcal{F}:\mathcal{P}(\mathbf{R}^{n})\to(-\infty,\infty]$ is _proper_ if
there exists $\mu\in\mathcal{P}(\mathbf{R}^{n})$ such that
$\mathcal{F}(\mu)<\infty$, and is _lower semicontinuous_ (l.s.c.) if for any
weakly convergence sequence $\mu_{k}\to\mu$, we have
$\liminf_{k\to\infty}\mathcal{F}(\mu_{k})\geq\mathcal{F}(\mu)$.

######  Lemma A.11.

For any $\epsilon>0$, the functional
$\mathcal{E}_{\epsilon}:\mathcal{P}(\mathbf{R}^{n})\to(-\infty,\infty]$ is
proper and l.s.c. Moreover, if $X$ is compact, the minimum
$\min_{\mu\in\mathcal{P}(X)}\mathcal{E}_{\epsilon}(\mu)$ is attained by some
measure in $\mathcal{P}(X)$.

###### Proof.

Taking any $x\in X$, since $\phi_{\epsilon}$ is bounded and $p(x)>0$, we see
that $\mathcal{E}_{\epsilon}(\delta_{x})<\infty$ so $\mathcal{E}_{\epsilon}$
is proper. Moreover, given weakly convergence $\mu_{k}\to\mu$, by Portmanteau
theorem and the fact that $W_{\epsilon}$ is nonnegative and l.s.c., we
conclude that $\mathcal{E}_{\epsilon}$ is also l.s.c.

The set of probability distributions
$\mathcal{P}(X)\subset\mathcal{P}(\mathbf{R}^{n})$ is tight by the compactness
of $X$. It is closed since if $\\{\mu_{k}\\}\subset\mathcal{P}(X)$ weakly
converges to $\mu$, then by Portmanteau theorem, $\mu(X)\geq\limsup\mu(X)=1$
so that $\mu\in\mathcal{P}(X)$. Hence by Prokhorov’s theorem (Theorem 5.1.3
(Ambrosio et al., 2005)), $\mathcal{P}(X)$ is (sequentially) compact. It is
then an elementary result that any l.s.c. function attains its minimum on a
compact set. ∎

We next prove Proposition 3.5 regarding the convexity of
$\mathcal{E}_{\epsilon}$ and the uniqueness of its minimum.

###### Proof of Proposition 3.5.

Let $\mu$ be any finite signed Borel measure. The compactness assumption and
the fact that $p\in C^{1}(X)$ imply $p(x)>\delta$ for any $x\in X$ for some
$\delta>0$. Hence $p(x)^{-1/2}\leq\delta^{-1/2}$, so the weighted measure
$\tilde{\mu}$ defined by
$\mathop{}\\!\mathrm{d}\tilde{\mu}(x):=p^{-1/2}(x)\mathop{}\\!\mathrm{d}\mu(x)$
is also finite. By the definition of i.s.p.d. kernels, we have
$\mathcal{E}_{W_{\epsilon}}(\mu)=\mathcal{E}_{k_{\epsilon}}(\tilde{\mu})>0$ if
$\tilde{\mu}$ is not the zero measure, which is equivalent to $\mu$ not being
zero since $p>0$. Thus $W_{\epsilon}$ is i.s.p.d. on $X$. Also note that
$(p(x)p(y))^{-1/2}<\delta^{-1}$ for all $x,y\in X$, so
$\mathcal{E}_{W_{\epsilon}}$ is always finite. By Lemma 1.1 of Pronzato &
Zhigljavsky (2021), we conclude that $\mathcal{E}_{W_{\epsilon}}$ is strictly
convex in $\mathcal{M}_{\textup{sign}}$, the space of finite signed measures,
and in particular it is convex on $\mathcal{P}(X)$. Hence combined with the
existence result from Lemma A.11 we conclude $\mathcal{E}_{\epsilon}$ attains
a unique minimum in $\mathcal{P}(X)$. ∎

Our goal in the rest of this section is to prove the $\Gamma$-convergence of
$\mathcal{E}_{\epsilon}$ to
$\mathcal{E}:\mathcal{P}(\mathbf{R}^{n})\to[0,\infty]$ defined in (13). By
Lemma A.8, we have $\mathcal{E}(\mu)=\chi^{2}(\mu\;\|\;\mu^{*})+1$. Hence
together we will have proved Theorem 3.6.

######  Definition A.12 ($\Gamma$-convergence).

A sequence of functionals
$\mathcal{F}_{\epsilon}:\mathcal{P}(\mathbf{R}^{n})\to(-\infty,\infty]$ is
said to _$\Gamma$ -converge_ to
$\mathcal{F}:\mathcal{P}(\mathbf{R}^{n})\to(-\infty,\infty]$, denoted as
$\mathcal{F}_{\epsilon}\stackrel{{\scriptstyle\Gamma}}{{\to}}\mathcal{F}$, if:

  1. (a)

For any sequence $\mu_{\epsilon}\in\mathcal{P}(\mathbf{R}^{n})$ converging
weakly to $\mu\in\mathcal{P}(\mathbf{R}^{n})$, $\liminf_{\epsilon\to
0}\mathcal{F}_{\epsilon}(\mu_{\epsilon})\geq\mathcal{F}(\mu)$;

  2. (b)

For any $\mu\in\mathcal{P}(\mathbf{R}^{n})$, there exists a sequence
$\mu_{\epsilon}\in\mathcal{P}(\mathbf{R}^{n})$ converging weakly to $\mu$ with
$\limsup_{\epsilon\to
0}\mathcal{F}_{\epsilon}(\mu_{\epsilon})\leq\mathcal{F}(\mu)$.

We will show
$\mathcal{E}_{\epsilon}\stackrel{{\scriptstyle\Gamma}}{{\to}}\mathcal{E}$
using Fourier transforms and Bochner’s theorem.

######  Definition A.13 (Fourier transform).

For $f\in L^{1}(\mathbf{R}^{n})$, its _Fourier transform_ $\hat{f}$ is the
complex-valued function defined via

$\hat{f}(\xi):=\int e^{-2\pi i\xi\cdot x}f(x)\mathop{}\\!\mathrm{d}x.$

More generally, for a signed finite measure
$\mu\in\mathcal{M}_{\textup{sign}}(\mathbf{R}^{n})$, its _Fourier transform_
$\hat{\mu}$ is the complex-valued function defined via

$\hat{\mu}(\xi):=\int e^{-2\pi i\xi\cdot x}\mathop{}\\!\mathrm{d}\mu(x).$

This integral is always well-defined and moreover $\hat{\mu}$ is uniformly
continuous; see Borodachov et al. (2019, Section 1.10).

We will prove the following weak version (under the additional assumption that
a mollifier $\phi$ is integrable) of Bochner’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, we cannot
directly use the Fourier inversion formula since it is not obvious how to
check the integrability of $\hat{\phi}$ when $\phi$ is a mollifier.

######  Lemma A.14.

Suppose $\phi\in L^{1}(\mathbf{R}^{n})$ is even, bounded, continuous, and
$k(x,y):=\phi(x-y)$ is i.s.p.d. Then its Fourier transform $\hat{\phi}$ is
real, nonnegative, and the following inversion formula holds:

$\phi(x)=\int e^{2\pi
ix\cdot\xi}\hat{\phi}(\xi)\mathop{}\\!\mathrm{d}\xi\quad\text{ for all
$x\in\mathbf{R}^{n}$}.$

###### Proof.

The proof is adapted from that of Varadhan (2001, Theorem 2.7) and is extended
to the multi-dimensional case.

Since $\hat{\phi}(\xi):=\int e^{-2\pi
ix\cdot\xi}\phi(x)\mathop{}\\!\mathrm{d}x$ and $\phi(x)=\phi(-x)$, with a
change of variable $x^{\prime}=-x$, we obtain
$\overline{\hat{\phi}(\xi)}=\hat{\phi}(\xi)$, so $\hat{\phi}$ is real.

Next we show that $\hat{\phi}$ is nonnegative. For $T>0$, we compute, for a
fixed $\xi\in\mathbf{R}^{n}$,

$\displaystyle\frac{1}{T^{n}}\int_{[0,T]^{n}}\int_{[0,T]^{n}}e^{-2\pi
i(t-s)\cdot\xi}\phi(t-s)\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}s$
$\displaystyle=$
$\displaystyle\frac{1}{T^{n}}\int_{[-T,T]^{n}}\left(\int_{\prod_{i}[\absolutevalue{u_{i}},2T-\absolutevalue{u_{i}}]}2^{-n}\mathop{}\\!\mathrm{d}v\right)e^{-2\pi
iu\cdot\xi}\phi(u)\mathop{}\\!\mathrm{d}u$ $\displaystyle=$
$\displaystyle\frac{1}{T^{n}}\int_{[-T,T]^{n}}\left(\prod_{i=1}^{n}\frac{2T-2\absolutevalue{u_{i}}}{2}\right)e^{-2\pi
iu\cdot\xi}\phi(u)\mathop{}\\!\mathrm{d}u$ $\displaystyle=$
$\displaystyle\int_{[-T,T]^{n}}\left(\prod_{i=1}^{n}\left(1-\frac{\absolutevalue{u_{i}}}{T}\right)\right)e^{-2\pi
iu\cdot\xi}\phi(u)\mathop{}\\!\mathrm{d}u,$

where we have used change of variable $u=t-s$, $v=t+s$. Since $\phi\in
L^{1}(\mathbf{R}^{n})$, by dominated convergence theorem, we have as
$T\to\infty$

$\displaystyle\hat{\phi}(\xi)=\int e^{-2\pi i\xi\cdot
x}\phi(x)\mathop{}\\!\mathrm{d}x=\lim_{T\to\infty}\frac{1}{T^{n}}\int_{[0,T]^{n}}\int_{[0,T]^{n}}e^{-2\pi
i(t-s)\cdot\xi}\phi(t-s)\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}s.$

For $t,s\in\mathbf{R}^{n}$, we have

$\displaystyle\real\left(e^{-2\pi i(t-s)\cdot\xi}\phi(t-s))\right)$
$\displaystyle=$ $\displaystyle\cos(2\pi t\cdot\xi)k(t,s)\cos(2\pi
s\cdot\xi)+\sin(2\pi t\cdot\xi)k(t,s)\sin(2\pi s\cdot\xi).$

For a fixed $T$, if we define a finite measure $\mu$ as
$\mathop{}\\!\mathrm{d}\mu=\mathbf{1}_{t\in[0,T]^{n}}\cos(2\pi
t\cdot\xi)\mathop{}\\!\mathrm{d}t$, then $k$ being i.s.p.d. implies

$\iint\phi(t-s)\mathop{}\\!\mathrm{d}\mu\mathop{}\\!\mathrm{d}\mu=\int_{[0,T]^{n}}\int_{[0,T]^{n}}\cos(2\pi
t\cdot\xi)\phi(t-s)\cos(2\pi
s\cdot\xi)\mathop{}\\!\mathrm{d}t\mathop{}\\!\mathrm{d}s\geq 0,$

and similarly for $\sin$. Since $\hat{\phi}$ is real, we conclude that
$\hat{\phi}$ is nonnegative.

Finally we prove Lemma A.14. For $\sigma>0$, define
$\hat{\phi}_{\sigma}(\xi):=\hat{\phi}(\xi)e^{-\sigma^{2}\norm{\xi}^{2}_{2}}$.
Since $\hat{\phi}$ is bounded (because $\phi\in L^{1}(\mathbf{R}^{n})$), we
see that $\hat{\phi}_{\sigma}\in L^{1}(\mathbf{R}^{n})$. We compute, for
$x\in\mathbf{R}^{n}$, using Fubini’s theorem,

$\displaystyle\int e^{2\pi
ix\cdot\xi}\hat{\phi}_{\sigma}(\xi)\mathop{}\\!\mathrm{d}\xi$
$\displaystyle=\int\hat{\phi}(\xi)e^{-\sigma^{2}\norm{\xi}^{2}_{2}}e^{2\pi
ix\cdot\xi}\mathop{}\\!\mathrm{d}\xi$ $\displaystyle=\int\left(\int e^{-2\pi
iy\cdot\xi}\phi(y)\mathop{}\\!\mathrm{d}y\right)e^{-\sigma^{2}\norm{\xi}^{2}_{2}}e^{2\pi
ix\cdot\xi}\mathop{}\\!\mathrm{d}\xi$ $\displaystyle=\int\left(\int e^{-2\pi
i(y-x)\cdot\xi}e^{-\sigma^{2}\norm{\xi}_{2}^{2}}\mathop{}\\!\mathrm{d}\xi\right)\phi(y)\mathop{}\\!\mathrm{d}y$
$\displaystyle=\int(\pi/\sigma^{2})^{n/2}e^{-\pi^{2}\norm{x-y}^{2}_{2}/\sigma^{2}}\phi(y)\mathop{}\\!\mathrm{d}y,$

where we use the Fourier transform formula (Borodachov et al., 2019, (4.4.1))
of the Gaussian distribution. Notice that
$p_{\sigma}(y):=(\pi/\sigma^{2})^{n/2}e^{-\pi^{2}\norm{x-y}^{2}_{2}/\sigma^{2}}$
is the density of a multivariate Gaussian centered at $x$ with covariance
$\sigma^{2}/(2\pi^{2})\cdot\bm{I}$. Hence $\int e^{2\pi
ix\cdot\xi}\hat{\phi}_{\sigma}(\xi)\mathop{}\\!\mathrm{d}\xi=(p_{\sigma}*\phi)(0)$.
Since $\phi$ is bounded by assumption, with $x=0$ we find
$\int\hat{\phi}_{\sigma}(\xi)\mathop{}\\!\mathrm{d}\xi\leq\norm{\phi}_{\infty}$.
Taking $\sigma\to 0$, by monotone convergence theorem, we have
$\int\hat{\phi}(\xi)\mathop{}\\!\mathrm{d}\xi\leq\norm{\phi}_{\infty}$, so
together with the fact that $\hat{\phi}\geq 0$ we have $\hat{\phi}\in
L^{1}(\mathbf{R}^{n})$. Finally, since $\hat{\phi}_{\sigma}\leq\hat{\phi}$, by
dominated convergence theorem and Proposition A.4 (note $p_{\sigma}$ is
centered at $x$), we have

$\int e^{2\pi
ix\cdot\xi}\hat{\phi}(\xi)\mathop{}\\!\mathrm{d}\xi=\lim_{\sigma\to 0}\int
e^{2\pi
ix\cdot\xi}\hat{\phi}_{\sigma}(\xi)\mathop{}\\!\mathrm{d}\xi=\lim_{\sigma\to
0}(p_{\sigma}*\phi)(0)=\phi(x).$

∎

######  Proposition A.15.

Given $\phi:\mathbf{R}^{n}\to\mathbf{R}$ satisfying the assumptions of Lemma
A.14, for any $\nu\in\mathcal{M}_{\textup{sign}}(\mathbf{R}^{n})$, it holds
that

$\displaystyle\iint\phi(x-y)\mathop{}\\!\mathrm{d}\nu(x)\mathop{}\\!\mathrm{d}\nu(y)=\int\absolutevalue{\hat{\nu}(\xi)}^{2}\hat{\phi}(\xi)\mathop{}\\!\mathrm{d}\xi,$

###### Proof.

By Lemma A.14,

$\displaystyle\iint\phi(x-y)\mathop{}\\!\mathrm{d}\nu(x)\mathop{}\\!\mathrm{d}\nu(y)$
$\displaystyle=\iint\left(\int e^{-2\pi
i(x-y)\cdot\xi}\hat{\phi}(\xi)\mathop{}\\!\mathrm{d}\xi\right)\mathop{}\\!\mathrm{d}\nu(x)\mathop{}\\!\mathrm{d}\nu(y)$
$\displaystyle=\int\left(\int e^{-2\pi
ix\cdot\xi}\mathop{}\\!\mathrm{d}\nu(x)\right)\left(\int e^{2\pi
iy\cdot\xi}\mathop{}\\!\mathrm{d}\nu(y)\right)\hat{\phi}(\xi)\mathop{}\\!\mathrm{d}\xi$
$\displaystyle=\int\hat{\nu}(\xi)\overline{\hat{\nu}(\xi)}\hat{\phi}(\xi)\mathop{}\\!\mathrm{d}\xi=\int\absolutevalue{\hat{\nu}(\xi)}^{2}\hat{\phi}(\xi)\mathop{}\\!\mathrm{d}\xi.$

where we use Fubini’s theorem (all measures are finite and $e^{-2\pi i\cdot}$
is bounded) to exchange the order of integration. ∎

###### Proof of Theorem 3.6.

First note that Definition A.12(b) is immediate from Theorem 3.3 with
$\mu_{\epsilon}=\mu$: if $\mathcal{E}(\mu)=\infty$, then trivially
$\limsup_{\epsilon\to 0}\mathcal{E}_{\epsilon}(\mu)\leq\mathcal{E}(\mu)$;
otherwise we apply Theorem 3.3. So we focus on proving criterion Definition
A.12(a).

Fix a sequence $\mu_{\epsilon}\in\mathcal{P}(\mathbf{R}^{n})$ that converges
weakly to $\mu\in\mathcal{P}(\mathbf{R}^{n})$, and our goal is to show
$\liminf_{\epsilon\to
0}\mathcal{E}_{\epsilon}(\mu_{\epsilon})\geq\mathcal{E}(\mu)$. Without loss of
generality, we may assume $\mathcal{E}_{\epsilon}(\mu_{\epsilon})<\infty$ for
all $\epsilon>0$ (these terms have no effect in $\liminf$), which implies
$\mu_{\epsilon}(X)=1$. By Portmanteau’s theorem and the assumption that $X$ is
closed, we have $\mu(X)\geq\limsup_{\epsilon\to 0}\mu_{\epsilon}(X)=1$. So all
of $\mu_{\epsilon}$ and $\mu$ will have support in $X$.

For $m>0,\epsilon>0$, define a nonnegative measure $\nu_{\epsilon,m}$ by

$\mathop{}\\!\mathrm{d}\nu_{\epsilon,m}(x):=h_{m}(x)p^{-1/2}(x)\mathop{}\\!\mathrm{d}\mu_{\epsilon}(x),$

where $h_{m}:\mathbf{R}^{n}\to\mathbf{R}$ is a continuous monotonically
decreasing cutoff function satisfying $h_{m}(x)=1$ if $\norm{x}_{2}^{2}<m$ and
$h_{m}(x)=0$ if $\norm{x}_{2}^{2}>m+1$, and that $h_{m}(x)\leq
h_{m^{\prime}}(x)$ for $m<m^{\prime}$. Then since $p>0$ is continuous, it is
bounded below on any compact set in $X$, and hence $\nu_{\epsilon,m}$ is
finite. Also define, for $m>0$, a nonnegative measure $\nu_{m}$ by

$\mathop{}\\!\mathrm{d}\nu_{m}(x):=h_{m}(x)p^{-1/2}(x)\mathop{}\\!\mathrm{d}\mu(x),$

which is again finite following the same reasoning.

Then for any $m>0$,

$\displaystyle\mathcal{E}_{\epsilon}(\mu_{\epsilon})$
$\displaystyle=\iint\phi_{\epsilon}(x-y)\mathop{}\\!\mathrm{d}\nu_{\epsilon}(x)\mathop{}\\!\mathrm{d}\nu_{\epsilon}(y)$
$\displaystyle\geq\iint\phi_{\epsilon}(x-y)\mathop{}\\!\mathrm{d}\nu_{\epsilon,m}(x)\mathop{}\\!\mathrm{d}\nu_{\epsilon,m}(y)=\int\absolutevalue{\hat{\nu}_{\epsilon,m}(\xi)}^{2}\hat{\phi}_{\epsilon}(\xi)\mathop{}\\!\mathrm{d}\xi,$
(15)

where we apply Proposition A.15 for the last equality. On the other hand, note
that

$\displaystyle\hat{\nu}_{\epsilon,m}(\xi)=\int e^{-2\pi i\xi\cdot
x}\mathop{}\\!\mathrm{d}\nu_{\epsilon,m}(x)=\int e^{-2\pi i\xi\cdot
x}h_{m}(x)p^{-1/2}(x)\mathop{}\\!\mathrm{d}\mu_{\epsilon}(x).$

Since $\mu_{\epsilon}\to\mu$ weakly and the last integrand is a continuous
bounded function, we have

$\lim_{\epsilon\to 0}\hat{\nu}_{\epsilon,m}(\xi)=\int e^{-2\pi i\xi\cdot
x}h_{m}(x)p^{-1/2}(x)\mathop{}\\!\mathrm{d}\mu(x)=\hat{\nu}_{m}(\xi).$

On the other hand, $\hat{\phi}_{\epsilon}(\xi)=\int e^{-2\pi i\xi\cdot
x}\phi_{\epsilon}(x)\mathop{}\\!\mathrm{d}x$. Since by Definition 3.1,
$\phi_{\epsilon}\mathop{}\\!\mathrm{d}x$ converges to $\delta_{0}$ in
probability (and in particular weakly), we have $\lim_{\epsilon\to
0}\hat{\phi}_{\epsilon}(\xi)=1$.

Applying Fatou’s lemma to (15), we obtain, for any $m>0$,

$\liminf_{\epsilon\to
0}\mathcal{E}_{\epsilon}(\mu_{\epsilon})\geq\int\absolutevalue{\hat{\nu}_{m}(\xi)}^{2}\mathop{}\\!\mathrm{d}\xi.$

If
$\int\absolutevalue{\hat{\nu}_{m}(\xi)}^{2}\mathop{}\\!\mathrm{d}\xi=\infty$
for any $m>0$, then we are done since $\liminf_{\epsilon\to
0}\mathcal{E}_{\epsilon}(\mu_{\epsilon})=\infty$. Otherwise, by Kühn (2016,
Lemma 2.11), for every $m>0$, $\nu_{m}$ has density in
$L^{2}(\mathbf{R}^{n})$. This imples $\mu$ has density everywhere. Suppose
$\mathop{}\\!\mathrm{d}\mu(x)=q(x)\mathop{}\\!\mathrm{d}\mathcal{L}_{n}(x)$.
By Plancherel’s theorem and the monotone convergence theorem, we have

$\displaystyle\liminf_{\epsilon\to 0}\mathcal{E}_{\epsilon}(\mu_{\epsilon})$
$\displaystyle\geq\lim_{m\to\infty}\int\absolutevalue{\hat{\nu}_{m}(\xi)}^{2}\mathop{}\\!\mathrm{d}\xi$
$\displaystyle=\lim_{m\to\infty}\int_{X}\left(h_{m}(x)p^{-1/2}(x)q(x)\right)^{2}\mathop{}\\!\mathrm{d}x$
$\displaystyle=\int_{X}p^{-1}(x)q(x)^{2}\mathop{}\\!\mathrm{d}x=\mathcal{E}(\mu).$

This completes the proof of
$\mathcal{E}_{\epsilon}\stackrel{{\scriptstyle\Gamma}}{{\to}}\mathcal{E}$.

Now suppose $X$ is compact, and let
$\mu_{\epsilon}^{*}:=\arg\min_{\mu\in\mathcal{P}(X)}\mathcal{E}_{\epsilon}(\mu)$.
To establish $\mu^{*}_{\epsilon}\to\mu^{*}$ weakly (resp. $\lim_{\epsilon\to
0}\mathcal{E}_{\epsilon}(\mu^{*}_{\epsilon})=\mathcal{E}(\mu^{*})=1$), it
suffices to show that every subsequence of $\\{\mu^{*}_{\epsilon}\\}$ (resp.
$\\{\mathcal{E}_{\epsilon}(\mu_{\epsilon}^{*})\\}$) has a further convergence
subsequence converging to $\mu^{*}$ (resp. $\mathcal{E}(\mu^{*})$). With a
slight abuse of notation, we assume the sequence of $\epsilon$ is chosen so
that $\\{\mu_{\epsilon}^{*}\\}$ (resp.
$\\{\mathcal{E}_{\epsilon}(\mu_{\epsilon}^{*})\\}$) is already some
subsequence of the original sequence. As argued in the proof of Lemma A.11,
$\mathcal{P}(X)$ is compact with respect to weak convergence. Hence
$\\{\mu_{\epsilon}^{*}\\}$ has a weakly convergence subsequence
$\mu_{\epsilon_{k}}^{*}\to\nu$ for some $\nu\in\mathcal{P}(X)$. The
$\Gamma$-convergence
$\mathcal{E}_{\epsilon}\stackrel{{\scriptstyle\Gamma}}{{\to}}\mathcal{E}$
implies

$\displaystyle\liminf_{k\to\infty}\mathcal{E}_{\epsilon_{k}}(\mu^{*})\geq\liminf_{k\to\infty}\mathcal{E}_{\epsilon_{k}}(\mu_{\epsilon_{k}}^{*})\geq\mathcal{E}(\nu)\geq\limsup_{k\to\infty}\mathcal{E}_{\epsilon_{k}}(\nu)\geq\limsup_{k\to\infty}\mathcal{E}_{\epsilon_{k}}(\mu_{\epsilon_{k}}^{*}),$

where we have used, for each inequality,
$\mu_{\epsilon_{k}}^{*}=\operatorname*{arg\,min}_{\mathcal{P}(X)}\mathcal{E}_{\epsilon_{k}}$,
Definition A.12((a)), the first paragraph of this proof, and again the fact
that
$\mu_{\epsilon_{k}}^{*}=\operatorname*{arg\,min}_{\mathcal{P}(X)}\mathcal{E}_{\epsilon_{k}}$.
Since
$\liminf_{k\to\infty}\mathcal{E}_{\epsilon_{k}}(\mu^{*})=\mathcal{E}(\mu^{*})$
by Theorem 3.3, we have

$\mathcal{E}(\mu^{*})\geq\lim_{k\to\infty}\mathcal{E}_{\epsilon_{k}}(\mu_{\epsilon_{k}}^{*})=\mathcal{E}(\nu)\geq\mathcal{E}(\mu^{*}),$

where the last inequality follows because $\mu^{*}$ is the minimizer of
$\mathcal{E}$. Then
$\lim_{k\to\infty}\mathcal{E}_{\epsilon_{k}}(\mu_{\epsilon_{k}}^{*})=\mathcal{E}(\mu^{*})=1$.
Moreover, $\chi^{2}(\nu\;\|\;\mu^{*})=0$. This can only happen if $\nu$ and
$\mu^{*}$ agree on all Borel sets, so $\nu=\mu^{*}$. ∎

###  A.4 Differential calculus of $\mathcal{E}_{\epsilon}$ in
$\mathcal{P}_{2}(\mathbf{R}^{n})$

######  Lemma A.16.

Let $x_{0}\in\mathbf{R}$, $h>0$, and $\Omega\subset\mathbf{R}^{m}$ be a
compact set. Suppose $f:(x_{0}-h,x_{0}+h)\times\Omega\to\mathbf{R}$ is jointly
continuous and the derivative
$\partialderivative{x}f:(x_{0}-h,x_{0}+h)\times\Omega\to\mathbf{R}$ exists and
is jointly continuous. Then
$\int_{\Omega}f(x,\omega)\mathop{}\\!\mathrm{d}\omega$ is differentiable for
$x\in(x_{0}-h,x_{0}+h)$, and

$\derivative{x}\int_{\Omega}f(x,\omega)\mathop{}\\!\mathrm{d}\omega=\int_{\Omega}\partialderivative{x}f(x,\omega)\mathop{}\\!\mathrm{d}\omega$

where the integration is with respect to the Lebesgue measure
$\mathcal{L}_{m}$ on $\Omega$.

###### Proof.

Fix $x\in(x_{0}-h,x_{0}+h)$ and let $t>0$ be small enough such that
$[x-t,x+t]\subset(x_{0}-h,x_{0}+h)$. Note that
$\int_{\Omega}f(x,\omega)\mathop{}\\!\mathrm{d}\omega$ is well-defined by the
dominated convergence theorem since
$\sup_{x\in[x-t,x+t],\omega\in\Omega}f(x,\omega)<\infty$ and
$\mathcal{L}_{m}(\Omega)<\infty$. Define
$\theta(\omega):=\sup_{x\in[x-t,x+t]}\absolutevalue{\partialderivative{x}f(x,\omega)}$
which is finite since $\partialderivative{x}f(x,\omega)$ is continuous, so
that $\theta(\omega)\geq\partialderivative{x}\absolutevalue{f(x,\omega)}$ for
all $x\in[x-t,x+t]$ and $\theta$ is integrable since
$\mathcal{L}_{m}(\Omega)<\infty$. Hence by the differentiation lemma (Klenke,
2013, Theorem 6.28) we are done. ∎

The following lemma is similar to Korba et al. (2021, Proposition 3) but with
different assumptions: we do not put any integrability assumptions on
$W_{\epsilon}$, but we do restrict measures to have compact support. For a
differentiable $f:\mathbf{R}^{n}\times\mathbf{R}^{n}\to\mathbf{R}$, we use
$\nabla_{1}f(x,y)$ to denote the gradient with respect to $x$. We also use
$\mathrm{H}_{1}f(x,y)$ to denote the Hessian of $f$ with respect to $x$. We
similarly denote $\nabla_{2}f(x,y)$ and $\mathrm{H}_{2}f(x,y)$.

######  Lemma A.17.

Assume $\mu\in\mathcal{P}_{2}(\mathbf{R}^{n})$ has density $q\in
C_{c}^{1}(\mathbf{R}^{n})$. Let $\bm{\xi}\in
C_{c}^{\infty}(\mathbf{R}^{n};\mathbf{R}^{n})$. Denote
$\bm{s}(x,t):=x+t\bm{\xi}(x)$. Then for all $t>0$ and all $\epsilon>0$,

$\displaystyle\derivative{t}\mathcal{E}_{\epsilon}(\mu_{t})$
$\displaystyle=2\iint\left(\bm{\xi}(x)^{\top}\nabla_{1}W_{\epsilon}(\bm{s}(x,t),\bm{s}(y,t))\right)\mathop{}\\!\mathrm{d}\mu(x)\mathop{}\\!\mathrm{d}\mu(y).$
(16)

In particularly,

$\displaystyle\derivative{t}\bigg{|}_{t=0}\mathcal{E}_{\epsilon}(\mu_{t})$
$\displaystyle=2\iint\left(\bm{\xi}(x)^{\top}\nabla_{1}W_{\epsilon}(x,y)\right)\mathop{}\\!\mathrm{d}\mu(x)\mathop{}\\!\mathrm{d}\mu(y).$
(17)

Moreover,

$\displaystyle\derivative[2]{t}\bigg{|}_{t=0}\mathcal{E}_{\epsilon}(\mu_{t})$
$\displaystyle=2\iint\left(\bm{\xi}(x)^{\top}\nabla_{1}\nabla_{2}W_{\epsilon}(x,y)\bm{\xi}(y)+\bm{\xi}(x)^{\top}\mathrm{H}_{1}W_{\epsilon}(x,y)\bm{\xi}(x)\right)\mathop{}\\!\mathrm{d}\mu(x)\mathop{}\\!\mathrm{d}\mu(y).$
(18)

###### Proof.

We compute

$\displaystyle\derivative{t}\mathcal{E}_{\epsilon}(\mu_{t})$
$\displaystyle=\derivative{t}\left(2\iint
W_{\epsilon}(\bm{s}(x,t),\bm{s}(y,t))\mathop{}\\!\mathrm{d}\mu(x)\mathop{}\\!\mathrm{d}\mu(y)\right).$

Since $\mu$ has compact support and $(x,y,t)\mapsto
W_{\epsilon}(\bm{s}(x,t),\bm{s}(y,t))q(x)q(y)$ is jointly continuous and its
derivative with respect to $t$ is also jointly continuous, by Lemma A.16, we
can push $\derivative{t}$ inside the double integral and we obtain (16).
Another application of Lemma A.16 shows that if we take derivative with
respect to $t$ again on (16) and evaluate at $0$ we obtain (18). ∎

######  Lemma A.18.

Let $f\in C^{1}(\mathbf{R}^{n})$, $p\in[1,\infty]$. Assume $\phi_{\epsilon}$
has compact support for some $\epsilon>0$. Then for all $i=1,\ldots,n$,
$\phi_{\epsilon}*f$ is differentiable and

$\partialderivative{x_{i}}(\phi_{\epsilon}*f)(x)=\left(\phi_{\epsilon}*\partialderivative{x_{i}}f\right)(x).$

###### Proof.

Since $\operatorname*{\textup{supp}}(\phi_{\epsilon})$ is compact and $f\in
C^{1}(\mathbf{R}^{n})$ is bounded on any compact set, we know
$\phi_{\epsilon}*f$ is well-defined at every $x\in\mathbf{R}^{n}$. Note that

$\displaystyle\partialderivative{x_{i}}(\phi_{\epsilon}*f)(x)$
$\displaystyle=\partialderivative{x_{i}}\left(\int
f(x-y)\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y\right)$
$\displaystyle\stackrel{{\scriptstyle(?)}}{{=}}\int\partialderivative{f}{x_{i}}\/(x-y)\phi_{\epsilon}(y)\mathop{}\\!\mathrm{d}y=\left(\phi_{\epsilon}*\partialderivative{x_{i}}f\right)(x).$

Since $\operatorname*{\textup{supp}}(\phi_{\epsilon})$ is compact and
$(x,y)\mapsto f(x-y)\phi(y)$ is $C^{1}$, by Lemma A.16 we justify the
existence of the derivative and the exchange of differentiation and
integration (?). ∎

Recall that a vector field $\bm{w}\in L^{2}(\mu;\mathbf{R}^{n})$ is a _strong
Fréchet subdifferential_ of a functional
$\mathcal{F}:\mathcal{P}_{2}(\mathbf{R}^{n})\to(-\infty,+\infty]$ if for all
$\bm{T}\in L^{2}(\mu;\mathbf{R}^{n})$, the following holds (Definition 10.1.1
of Ambrosio et al. (2005)):

$\displaystyle\mathcal{F}(\bm{T}_{\\#}\mu)-\mathcal{F}(\mu)\geq\int\bm{w}(x)^{\top}(\bm{T}(x)-x)\mathop{}\\!\mathrm{d}\mu(x)+o\left(\norm{\bm{T}-\bm{I}}_{L^{2}(\mu;\mathbf{R}^{n})}\right).$
(19)

Note that we cannot apply Ambrosio et al. (2005, Lemma 10.4.1) directly to
prove (4) because interaction energies cannot be written as (10.4.1) in their
setup.

###### Proof of Proposition 3.7.

Let $\bm{\xi}\in C^{\infty}_{c}(\mathbf{R}^{n};\mathbf{R}^{n})$. By Lemma
A.17, we have

$\displaystyle\derivative{t}\bigg{|}_{t=0}\mathcal{E}_{\epsilon}(\mu_{t})$
$\displaystyle=2\iint\left(\bm{\xi}(x)^{\top}\nabla_{1}W_{\epsilon}(x,y)\right)\mathop{}\\!\mathrm{d}\mu(x)\mathop{}\\!\mathrm{d}\mu(y)$
$\displaystyle=2\iint\left(\bm{\xi}(x)^{\top}\nabla_{1}\left(\phi_{\epsilon}(x-y)(p(x)p(y))^{-1/2}\right)\right)q(y)\mathop{}\\!\mathrm{d}y\mathop{}\\!\mathrm{d}\mu(x)$
$\displaystyle=2\int\bm{\xi}(x)^{\top}\nabla\left(p(x)^{-1/2}(\phi_{\epsilon}*q/\sqrt{p})(x)\right)\mathop{}\\!\mathrm{d}\mu(x),$

where the last step follows from applying Lemma A.16 since $q$ has compact
support. Now suppose $\bm{w}\in L^{2}(\mu;\mathbf{R}^{n})$ is a strong Fréchet
subdifferential satisfying (19). For the sequence $\\{\bm{T}_{t}\\}$, we have
by definition

$\displaystyle\mathcal{E}_{\epsilon}(\mu_{t})-\mathcal{E}_{\epsilon}(\mu)$
$\displaystyle\geq\int\bm{w}(x)^{\top}(\bm{T}_{t}(x)-x)\mathop{}\\!\mathrm{d}\mu(x)+o\left(\norm{\bm{T}_{t}-\bm{I}}_{L^{2}(\mu;\mathbf{R}^{n})}\right)$
$\displaystyle=\int\bm{w}(x)^{\top}(t\bm{\xi}(x))\mathop{}\\!\mathrm{d}\mu(x)+o(t).$

Hence

$\liminf_{t\downarrow
0}\frac{\mathcal{E}_{\epsilon}(\mu_{t})-\mathcal{E}_{\epsilon}(\mu)}{t}\geq\int\bm{w}(x)^{\top}\bm{\xi}(x)\mathop{}\\!\mathrm{d}\mu(x)\geq\liminf_{t\uparrow
0}\frac{\mathcal{E}_{\epsilon}(\mu_{t})-\mathcal{E}_{\epsilon}(\mu)}{t}.$

The previous calculation shows
$\derivative{t}|_{t=0}\mathcal{E}_{\epsilon}(\mu_{t})$ exists, and hence it is
equal to $\int\bm{w}(x)^{\top}\bm{\xi}(x)\mathop{}\\!\mathrm{d}\mu(x)$. This
proves for any $\bm{\xi}\in C_{c}^{\infty}(\mathbf{R}^{n};\mathbf{R}^{n})$,

$\displaystyle\int\bm{w}(x)^{\top}\bm{\xi}(x)\mathop{}\\!\mathrm{d}\mu(x)$
$\displaystyle=\int\bm{\xi}(x)^{\top}\nabla\left(2p(x)^{-1/2}(\phi_{\epsilon}*q/\sqrt{p})(x)\right)\mathop{}\\!\mathrm{d}\mu(x).$

Hence we have shown (4).

Finally we show $\lim_{\epsilon\to
0}\bm{w}_{\epsilon}(x)=\bm{w}_{\chi^{2}}(x)$ for $\mu$-a.e. $x$ under the
additional assumption that $\phi_{\epsilon}$ has compact support. By Ambrosio
et al. (2005, Lemma 10.4.1) with
$F(x,\rho(x))=\left(\frac{\rho(x)}{p(x)}-1\right)^{2}p(x)$, we find the strong
subdifferential of $\chi^{2}(\cdot\;\|\;\mu^{*})$ is the one in (5). To show
(6), we compute for $\mu$-a.e. $x$,

$\displaystyle\bm{w}_{\epsilon}(x)$
$\displaystyle=\nabla\left(p(x)^{-1/2}(\phi_{\epsilon}*q/\sqrt{p})(x)\right)$
$\displaystyle=\nabla(p(x)^{-1/2})(\phi_{\epsilon}*q/\sqrt{p})(x)+p(x)^{-1/2}\nabla(\phi_{\epsilon}*q/\sqrt{p})(x)$
$\displaystyle=\nabla(p(x)^{-1/2})(\phi_{\epsilon}*q/\sqrt{p})(x)+p(x)^{-1/2}(\phi_{\epsilon}*\nabla(q/\sqrt{p}))(x),$

where for the last equality we have applied Lemma A.18. Now taking
$\epsilon\to 0$, by Proposition A.4 using the fact that
$\operatorname*{\textup{supp}}(\phi_{\epsilon})$ is compact (so that
$q/\sqrt{p}$ and $\nabla(q/\sqrt{p})$ are bounded on
$x+\operatorname*{\textup{supp}}(\phi_{\epsilon})$), we obtain (6). ∎

The statement of Proposition 3.8 is similar in form as Korba et al. (2021,
Corollary 4) but in our case we do not have the second term in Equation 18
vanishing and we need to take the limit $\epsilon\to 0$. We also do not resort
to RKHS theory in the proof.

###### Proof of Proposition 3.8.

By (18), we have

$\displaystyle\derivative[2]{t}\bigg{|}_{t=0}\mathcal{E}_{\epsilon}(\mu_{t})$
$\displaystyle=2(\mathcal{F}_{\epsilon}+\mathcal{G}_{\epsilon}),$

where

$\displaystyle F_{\epsilon}$
$\displaystyle=\iint\left(\bm{\xi}(x)^{\top}\nabla_{1}\nabla_{2}W_{\epsilon}(x,y)\bm{\xi}(y)\right)\mathop{}\\!\mathrm{d}\mu^{*}(x)\mathop{}\\!\mathrm{d}\mu^{*}(y)$
$\displaystyle G_{\epsilon}$
$\displaystyle=\iint\left(\bm{\xi}(x)^{\top}\mathrm{H}_{1}W_{\epsilon}(x,y)\bm{\xi}(x)\right)\mathop{}\\!\mathrm{d}\mu^{*}(x)\mathop{}\\!\mathrm{d}\mu^{*}(y).$

We tackle $F_{\epsilon}$ first. Observe that successive application of
integration by parts using the fact that $\bm{\xi}$ has compact support gives

$\displaystyle F_{\epsilon}$
$\displaystyle=\sum_{i,j=1}^{n}\iint\bm{\xi}_{i}(x)\partialderivative{}{x_{i}}{y_{j}}W_{\epsilon}(x,y)\bm{\xi}_{j}(y)p(x)p(y)\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y$
$\displaystyle=-\sum_{i,j=1}^{n}\iint\partialderivative{x_{i}}\/(\bm{\xi}_{i}(x)p(x))\partialderivative{y_{j}}W_{\epsilon}(x,y)\bm{\xi}_{j}(y)p(y)\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y$
$\displaystyle=\sum_{i,j=1}^{n}\iint\partialderivative{x_{i}}\/(\bm{\xi}_{i}(x)p(x))W_{\epsilon}(x,y)\partialderivative{y_{j}}\/(\bm{\xi}_{j}(y)p(y))\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y$
$\displaystyle=\iint\left(\sum_{i=1}^{n}\partialderivative{x_{i}}\/(\bm{\xi}_{i}(x)p(x))\right)W_{\epsilon}(x,y)\left(\sum_{j=1}^{n}\partialderivative{y_{j}}\/(\bm{\xi}_{j}(y)p(y))\right)\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y.$

If we view $\sum_{i=1}^{n}\partialderivative{x_{i}}\/(\bm{\xi}_{i}(x)p(x))$ as
the density of a signed measure (it is integrable since it has compact
support), and since $W_{\epsilon}$ is i.s.p.d. on the support of $\bm{\xi}$ by
Proposition 3.5(a), we see that each double integral in the last expression is
non-negative. Hence $F_{\epsilon}\geq 0$.

Next we show $\lim_{\epsilon\to 0}G_{\epsilon}=0$. Since $\mu^{*}$ has compact
support, by Fubini’s theorem,

$\displaystyle G_{\epsilon}$
$\displaystyle=\int\bm{\xi}(x)^{\top}\mathrm{H}_{1}\left(\int
W_{\epsilon}(x,y)p(y)\mathop{}\\!\mathrm{d}y\right)\bm{\xi}(x)p(x)\mathop{}\\!\mathrm{d}x.$

To expand the integral inside the Hessian operator, we have

$\displaystyle\mathrm{H}_{1}\left(\int
W_{\epsilon}(x,y)p(y)\mathop{}\\!\mathrm{d}y\right)$
$\displaystyle=\mathrm{H}\left(p(x)^{-1/2}(\phi_{\epsilon}*\sqrt{p})(x)\right).$

First by the chain rule and Lemma A.18, we have

$\displaystyle\derivative{x}\left(p(x)^{-1/2}(\phi_{\epsilon}*\sqrt{p})(x)\right)$
$\displaystyle=\derivative{x}\left(p(x)^{-1/2}\right)(\phi_{\epsilon}*\sqrt{p})+p(x)^{-1/2}\left(\phi_{\epsilon}*\derivative{x}\sqrt{p}(x)\right).$

Differentiating again while applying Lemma A.18, we obtain after rearranging
terms,

$\displaystyle\mathrm{H}_{1}\left(\int
W_{\epsilon}(x,y)p(y)\mathop{}\\!\mathrm{d}y\right)$
$\displaystyle=\left(\derivative[2]{x}p(x)^{-1/2}\right)(\phi_{\epsilon}*\sqrt{p})(x)$
$\displaystyle+2\left(\derivative{x}p(x)^{-1/2}\right)\left(\phi_{\epsilon}*\derivative{x}\sqrt{p}\right)(x)$
$\displaystyle+p(x)^{-1/2}\left(\phi_{\epsilon}*\derivative[2]{x}\sqrt{p}\right)(x).$
(20)

By Proposition A.4 and the fact that $p\in C^{2}_{c}(\mathbf{R}^{n})$, we have

$\displaystyle\lim_{\epsilon\to 0}\mathrm{H}_{1}\left(\int
W_{\epsilon}(x,y)p(y)\mathop{}\\!\mathrm{d}y\right)$ $\displaystyle=$
$\displaystyle\left(\derivative[2]{x}p(x)^{-1/2}\right)\sqrt{p}(x)+2\left(\derivative{x}p(x)^{-1/2}\right)\derivative{x}\sqrt{p}(x)+p(x)^{-1/2}\derivative[2]{x}\sqrt{p}(x)$
$\displaystyle=$
$\displaystyle\derivative[2]{x}\left(p(x)^{-1/2}\sqrt{p(x)}\right)=0.$

Finally, we have

$\displaystyle\lim_{\epsilon\to 0}G_{\epsilon}$
$\displaystyle=\lim_{\epsilon\to
0}\int\bm{\xi}(x)^{\top}\mathrm{H}_{1}\left(\int
W_{\epsilon}(x,y)p(y)\mathop{}\\!\mathrm{d}y\right)\bm{\xi}(x)p(x)\mathop{}\\!\mathrm{d}x$
$\displaystyle=\int\bm{\xi}(x)^{\top}\left(\lim_{\epsilon\to
0}\mathrm{H}_{1}\left(\int
W_{\epsilon}(x,y)p(y)\mathop{}\\!\mathrm{d}y\right)\right)\bm{\xi}(x)p(x)\mathop{}\\!\mathrm{d}x$
$\displaystyle=0,$

where interchanging the limit and the integral is jusfitied by the dominated
convergence theorem and the fact that $p$ and $\phi_{\epsilon}$ have compact
support (we need compact support assumption of $\phi_{\epsilon}$ to make sure
convolutions appearing in (20) are uniformly bounded when $\epsilon$ is
sufficiently small).

∎

##  Appendix B Weighted hypersingular Riesz energy

We recall results most relevant to us from Borodachov et al. (2019). Suppose
$X\subset\mathbf{R}^{n}$ is compact, of Hausdorff dimension $d$, and
$d$-rectifiable, i.e., the image of a bounded set in $\mathbf{R}^{d}$ under a
Lipschitz mapping. Given a non-vanishing continuous probability density $p$ on
$X$, define a measure
$\mathcal{H}_{d}^{p}(B):=\int_{B}p(x)\mathop{}\\!\mathrm{d}\mathcal{H}_{d}(x)$
for any Borel set $B\subset X$. The target measure $h_{d}^{p}$ is defined to
be
$h_{d}^{p}(B):=\nicefrac{{\mathcal{H}_{d}^{p}(B)}}{{\mathcal{H}_{d}^{p}(X)}}$.
The _weighted $s$-Riesz energy_ is defined as the interaction energy

$\displaystyle E_{s}(\omega_{N}):=\sum_{i\neq
j}\frac{(p(x_{i})p(x_{j}))^{-s/2d}}{\norm{x_{i}-x_{j}}^{s}}.$ (21)

The following result states that minimizers of the weighted $s$-Riesz energy
approximate the target measure when the $s$ is sufficiently large.

######  Theorem B.1 (Theorem 11.1.2, Borodachov et al. (2019)).

Suppose $s>d$. For $\omega_{N}^{*}\in\operatorname*{arg\,min}E_{s}$ with
$\omega_{N}^{*}=\\{x_{1}^{N},\ldots,x_{N}^{N}\\}$, we have weak convergence
$\delta_{\omega_{N}^{*}}\to h_{d}^{p}$ as $N\to\infty$.

A similar result with a slightly different assumption holds for $s=d$ (Theorem
11.1.3 of Borodachov et al. (2019)).

#### Comparison of (21) with discretized MIE (7).

Our theory is only valid for the full-dimensional case, i.e., $n=d$ (see a few
exceptions discussed in Remark A.9). When this is the case and if the
mollifier is taken to be from the $s$-Riesz family, (7) becomes, for $s>n$,

$\displaystyle
E_{\epsilon}(\omega_{N}):=\sum_{i,j=1}^{N}\frac{(p(x_{i})p(x_{j}))^{-1/2}}{\left(\norm{x_{i}-x_{j}}^{2}_{2}+\epsilon^{2}\right)^{s/2}}.$

Compared with (21), in our case there is an $\epsilon$ in the denominator, so
that the continuous version of the energy $\mathcal{E}_{\epsilon}$ does not
blow up all the time (this is in contrast with the continuous version of
(21)—see Borodachov et al. (2019, Theorem 4.3.1)). Moreover, the exponential
scaling on $p$ is different: in our case we use $-1/2$ whereas in (21) it is
$-s/2n$. The diagonal $i=j$ is included in our energy, but this is
inconsequential as another valid discretization is to discard the diagonal
term (Borodachov et al., 2019, Theorem 4.2.2). On the other hand, our theory
allows a bigger class of mollifiers that are not necessarily of Riesz
families. 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 $\mathbf{R}^{n}$.

##  Appendix C Experiment details and additional results

###  C.1 Gaussians in varying dimensions

We generate the $n$-dimensional Gaussians used to produce Figure 1 as follows.
We generate a matrix $\sqrt{A}\in\mathbf{R}^{n\times n}$ with i.i.d. entries
uniformly in $[-1,1]$. Then we set $A=\sqrt{A}\sqrt{A}^{\top}/\det(\sqrt{A})$.
This way $\det(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.

We visualize the samples from each method for $n=2$ in Figure 5.

Figure 5: Visualization of samples from each method for a 2D skewed Gaussian.

We notice that MIED is capable of generating well-separated samples, while for
SVGD there is a 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.

###  C.2 20-dimensional Dirichlet distribution

We sample from the 20-dimensional Dirichlet distribution in the same setup as
in Shi et al. (2021) with 50 particles. Results and visualization are shown in
Figure 6. We see that unlike sampling from a 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 a well-tested choice for simplex
constraint, yet obtaining a good mirror map for a generic constraint, even if
it is linear, can be challenging. Our method does not have such limitation as
it can easily incorporate exisiting constrained optimization tools.

Figure 6: Top: Metrics vs. the number of iteration for sampling the
20-dimensional Dirichlet distribution. Bottom: visualization of samples from
each method.