# Sampling from the Mean-Field Stationary Distribution

 Yunbum Kook Matthew S. Zhang
Sinho Chewi Murat A. Erdogdu Mufan (Bill) Li  School of Computer Science at
Georgia Institute of Technology<EMAIL_ADDRESS>Department of Computer
Science at University of Toronto, and Vector Institute,
<EMAIL_ADDRESS>School of Mathematics at Institute for Advanced
Study<EMAIL_ADDRESS>Department of Computer Science at University of Toronto,
and Vector Institute<EMAIL_ADDRESS>Department of Operations Research
and Financial Engineering at Princeton University<EMAIL_ADDRESS>
###### Abstract

We study the complexity of sampling from the stationary distribution of a
mean-field SDE, or equivalently, the complexity of minimizing a functional
over the space of probability measures which includes an interaction term. Our
main insight is to _decouple_ the two key aspects of this problem: (1)
approximation of the mean-field SDE via a finite-particle system, via uniform-
in-time propagation of chaos, and (2) sampling from the finite-particle
stationary distribution, via standard log-concave samplers. Our approach is
conceptually simpler and its flexibility allows for incorporating the state-
of-the-art for both algorithms and theory. This leads to improved guarantees
in numerous settings, including better guarantees for optimizing certain two-
layer neural networks in the mean-field regime.

##  1 Introduction

The minimization of energy functionals $\mathcal{E}$ over the Wasserstein
space $\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$ of probability measures has
attracted substantial research activity in recent years, encompassing numerous
application domains, including distributionally robust optimization [Kuh+19,
YKW22], sampling [JKO98, Wib18, Che23], and variational inference [LW16,
Lam+22, Dia+23, JCP23, Lac23, YY23].

A canonical example of such a functional is $\mathcal{E}(\mu)=\int
V\,\mathrm{d}\mu+\int\log\mu\,\mathrm{d}\mu$, where
$V:\mathbb{R}^{d}\to\mathbb{R}$ is called the potential. Up to an additive
constant, which is irrelevant for the optimization, this energy functional
equals the KL divergence $\operatorname{\mathsf{KL}}(\mu\mathbin{\|}\pi)$ with
respect to the density $\pi\propto\exp(-V)$, and the celebrated result of
[JKO98] identifies the Wasserstein gradient flow of $\mathcal{E}$ with the
Langevin diffusion. This link has inspired a well-developed theory for log-
concave sampling, with applications to Bayesian inference and randomized
algorithms; see [Che23] for an exposition.

The energy functional above contains two terms, corresponding to two of the
fundamental examples of functionals considered in Villani’s well-known
treatise on optimal transport [Vil03]. Namely, they are the “potential energy”
and the entropy, the latter being a special case of the “internal energy.”
However, Villani identifies a third fundamental functional—the “interaction
energy”—with the _pairwise_ form given by

$\displaystyle\mathcal{E}(\mu)\coloneqq\int V(x)\,\mu(\mathrm{d}x)+\iint
W(x-y)\,\mu(\mathrm{d}x)\,\mu(\mathrm{d}y)+\frac{\sigma^{2}}{2}\int\log\mu(x)\,\mu(\mathrm{d}x)\,.$
($\mathsf{pE}$)

More generally, in this work we consider minimizing the _generic_ entropy-
regularized energy

$\displaystyle\mathcal{E}(\mu)\coloneqq\mathcal{F}(\mu)+\frac{\sigma^{2}}{2}\int\log\mu\,\mathrm{d}\mu\,$
($\mathsf{gE}$)

where $\mathcal{F}:\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})\to\mathbb{R}$ is a
known functional. The minimization of the energy ($\mathsf{gE}$) has recently
been of interest due to its role in analysing neural network training dynamics
in the mean-field regime, including with [SNW22] and without [CB18, MMN18]
entropic regularization, as well as with Fisher regularization [Cla+23].

For the sake of exposition, let us first focus on minimizing the pairwise
energy ($\mathsf{pE}$). A priori, this question is more difficult than log-
concave sampling; for instance, $\pi$ does not admit a closed form but rather
is the solution to a non-linear equation

$\displaystyle\pi(x)\propto\exp\Bigl{(}-\frac{2}{\sigma^{2}}\,V(x)-\frac{2}{\sigma^{2}}\int
W(x-\cdot)\,\mathrm{d}\pi\Bigr{)}\,.$ (1.1)

However, here too there is a well-developed mathematical theory which suggests
a principled algorithmic approach. Just as the Wasserstein gradient flow of
($\mathsf{pE}$) in the case when $W=0$ can be identified with the Langevin
diffusion, the Wasserstein gradient flow of ($\mathsf{pE}$) in the case when
$W\neq 0$ corresponds to a (pairwise) McKean–Vlasov SDE, i.e., an SDE whose
coefficients depend on the marginal law of the process, given as

$\displaystyle\mathrm{d}X_{t}=-\Bigl{(}\nabla V(X_{t})+\int\nabla
W(X_{t}-\cdot)\,\mathrm{d}\pi_{t}\Bigr{)}\,\mathrm{d}t+\sigma\,\mathrm{d}B_{t}\,,$
($\mathsf{pMV}$)

where $\pi_{t}=\operatorname{\mathsf{law}}(X_{t})$, $W$ is even, and
$\\{B_{t}\\}_{t\geq 0}$ is a standard Brownian motion on $\mathbb{R}^{d}$.
Since the McKean–Vlasov SDE is the so-called _mean-field limit_ of interacting
particle systems, we can approximately sample from the minimizer $\pi$ by
numerically discretizing a system of SDEs, which describe the evolution of $N$
particles$\\{X_{t}^{1:N}\\}_{t\geq
0}\coloneqq\\{(X^{1}_{t},\dotsc,X^{N}_{t})\\}_{t\geq 0}$ as:

$\displaystyle\mathrm{d}X_{t}^{i}=-\Bigl{(}\nabla
V(X_{t}^{i})+\frac{1}{N-1}\sum_{j\in[N]\backslash i}\nabla
W(X_{t}^{i}-X_{t}^{j})\Bigr{)}\,\mathrm{d}t+\sigma\,\mathrm{d}B_{t}^{i}\,,\quad\forall\,i\in[N]\,,$
($\mathsf{pMV}_{N}$)

where $\\{B^{i}:i\in[N]\\}$ is a collection of independent Brownian motions.
Moreover, the error from approximating the mean-field limit via this finite
particle system has been studied in the literature on _propagation of chaos_
[Szn91]. Similarly, the Wasserstein gradient flow for ($\mathsf{gE}$)
corresponds to the _mean-field Langevin dynamics_ and admits an analogous
particle approximation.

The bounds for propagation of chaos have been refined over time, with [LL23]
recently establishing a tight error dependence
$\mathcal{O}(\nicefrac{{1}}{{N}})$ on the total number of particles $N$. These
bounds, however, do not translate immediately into algorithmic guarantees.
Existing sampling analyses study the propagation of chaos and discretization
as a single entangled problem, which thus far have only been able to use
weaker ${\mathcal{O}}(\sqrt{\nicefrac{{1}}{{N}}})$ rates for the former.
Furthermore, there has been recent interest in using more sophisticated
particle-based algorithms, e.g., “non-linear” Hamiltonian Monte Carlo [BS23]
and the mean-field underdamped Langevin dynamics [FW23] to reduce the
discretization error. Currently, this requires repeatedly carrying out the
propagation of chaos and time discretization analyses from the ground up for
each instance.

This motivates us to pose the following questions: (1) Can we incorporate
improvements in the propagation of chaos literature, such as the
$\mathcal{O}(\nicefrac{{1}}{{N}})$ error dependence shown in [LL23], to
improve existing theoretical guarantees? (2) Can we leverage recent advances
in the theory of log-concave sampling to design better algorithms?

Our main proposal in this work is to decouple the error into two terms,
representing the propagation of chaos and discretization errors respectively.
This simple and _modular_ approach immediately allows us to answer both
questions in the affirmative. Namely, we show how to combine established
propagation of chaos bounds in various settings [[, including the sharp rate
of]]lacker2023sharp with a large class of sophisticated off-the-shelf log-
concave samplers, such as interacting versions of the randomized midpoint
discretization of the underdamped Langevin dynamics [SL19, HBE20], Metropolis-
adjusted algorithms [Che+21, WSC22, AC23], and the proximal sampler [LST21,
Che+22, FYC23]. Our framework yields improvements upon prior state-of-the-art,
such as [BS23, FW23], and provides a clear path for future ones.

###  1.1 Contributions and Organization

##### Propagation of chaos at stationarity.

We provide three propagation of chaos results which hold in the
$\mathcal{W}_{2}$, $\sqrt{\operatorname{\mathsf{KL}}}$, and
$\sqrt{\operatorname{\mathsf{FI}}}$ “metrics”; the rates reflect the distance
of the $k$-particle marginal of the finite-particle system from $\pi^{\otimes
k}$: (1) In the setting of ($\mathsf{pE}$), under strong displacement
convexity, we obtain a $\mathcal{O}(\sqrt{\nicefrac{{k}}{{N}}})$ rate by
adapting techniques from [Szn91, Mal01]; (2) without assuming displacement
convexity, but assuming a weaker interaction, we obtain the sharp rate of
$\widetilde{\mathcal{O}}(\nicefrac{{k}}{{N}})$ following [LL23]; (3) finally,
in the general setting of ($\mathsf{gE}$), and assuming $\mathcal{F}$ is
convex along linear interpolations, we obtain a
$\mathcal{O}(\sqrt{\nicefrac{{k}}{{N}}})$ rate following [CRW22].

Unlike prior works, our proofs are carried out at stationarity; thus, our
proofs are _self-contained_ , streamlined, and include various improvements
(e.g., weaker assumptions and explicit bounds). As a result, our work also
serves as a helpful exposition to the mathematics of propagation of chaos.

##### Discretization.

Once the error due to particle approximation is controlled, we then obtain
improved complexity guarantees by applying recent advances in the theory of
log-concave sampling to the finite-particle stationary distribution. See Table
1 for a summary of our results, and the discussion in §4 for comparisons with
prior works and an application to neural network training.

Once again, the importance of our framework is its _modularity_ , which allows
for any combination of uniform-in-time propagation of chaos bounds and log-
concave sampler, provided that the finite-particle stationary distribution
satisfies certain isoperimetric properties needed for the sampling guarantees.
Toward this end, we also provide tools for verifying these isoperimetric
properties with constants that hold independently of the number of particles
(see §3.2.1).

###  1.2 Related Work

##### Mean-field equations.

The McKean–Vlasov SDE was first formulated in the works [McK66, Fun84, Mél96],
with origins dating to much earlier [Bol72]. It has applications in many
domains, from fluid dynamics [Vil02] to game theory [LL07, CD18]; see [CD22,
CD22a] for a comprehensive survey. The kinetic version of this equation is
known as the Boltzmann equation, and propagation of chaos has similarly been
studied under a variety of assumptions [BGM10, Mon17, GM21, GLM22]. One
prominent application within machine learning is the study of infinitely wide
two-layer neural networks in the mean-field regime (see §4.2).

##### Propagation of chaos and sampling for ($\mathsf{pE}$).

The original propagation of chaos arguments of [Szn91] were first made uniform
in time in [Mal01, Mal03] in both entropy and $\mathcal{W}_{2}$. The
aforementioned works all achieve an error of order
$\widetilde{\mathcal{O}}(\sqrt{\nicefrac{{k}}{{N}}})$, and require a strong
convexity assumption on $V$ and $W$. These were later adapted for non-smooth
potentials [JW17, JW18, BJW23]. Finally, [CRW22] obtained an entropic
propagation of chaos bound under a higher-order smoothness assumption. See
[CD22] for a more complete bibliography.

The breakthrough result of [Lac21, Lac22] obtained the sharp bound of
$\widetilde{\mathcal{O}}(\nicefrac{{k}}{{N}})$ when the interaction is
sufficiently weak, and this bound was made uniform in time in [LL23]. Their
approach differs significantly from previous proofs by considering a local
analysis based on the recursive BBGKY hierarchy. These results have been
extended to other divergences, e.g., the $\chi^{2}$ divergence, but without a
uniform-in-time guarantee [HR23]. In addition, [MRW24] showed an extension of
this result under a “convexity at infinity” assumption.

The question of sampling from minimizers of ($\mathsf{pE}$) was first studied
in [Tal96, BT97, AK02]. These works focused on the Euler–Maruyama
discretization of the finite-particle system ($\mathsf{pMV}_{N}$), under
$L^{\infty}$-boundedness of the gradients. Subsequently, the convergence of
the Euler–Maruyama scheme has been studied in many works, including but not
limited to [BH22, RES22, Li+23]. Finally, [BS23] considered a non-linear
version of Hamiltonian Monte Carlo; we give a detailed comparison with their
work in §4.

##### Propagation of chaos and sampling for ($\mathsf{gE}$).

The mean-field (underdamped) Langevin algorithm for minimizing ($\mathsf{gE}$)
was proposed and studied in [CRW22, Che+23]. Under alternative assumptions
(see §3.1.2), they established propagation of chaos with a
$\mathcal{O}(\sqrt{\nicefrac{{k}}{{N}}})$ rate, for both the overdamped and
the underdamped finite-particle approximations. Recent works from the machine
learning community [NWS22, SNW22, FW23, SWN23] studied the application of
these algorithms for optimizing two-layer neural networks and obtained
sampling guarantees. We provide a detailed comparison with their works in
§4.2.

##  2 Background and Notation

Let $\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$ be the set of probability
measures on $\mathbb{R}^{d}$ that admit a density with respect to the Lebesgue
measure and have finite second moment. We will also abuse notation and use the
same symbol for a measure and its density when there is no confusion. We use
superscripts for the particle index, and subscripts for the time variable. We
will use $\mathcal{O},\widetilde{\mathcal{O}}$ to signify upper bounds up to
numeric constants and polylogarithms respectively. We recall the definitions
of convexity and smoothness:

######  Definition 1.

A function $U:\mathbb{R}^{d}\to\mathbb{R}$ is $\alpha$-uniformly convex
(allowing for $\alpha\leq 0$) and $\beta$-smooth if the following hold
respectively

$\displaystyle\langle\nabla U(x)-\nabla U(y),x-y\rangle$
$\displaystyle\geq\alpha\,\lVert x-y\rVert^{2}$ $\displaystyle\text{for
all}~{}x,y\in\mathbb{R}^{d}\,,$ $\displaystyle\lVert\nabla U(x)-\nabla
U(y)\rVert$ $\displaystyle\leq\beta\,\lVert x-y\rVert$ $\displaystyle\text{for
all}~{}x,y\in\mathbb{R}^{d}\,.$

For two probability measures $\mu,\nu\in\mathcal{P}_{2,\rm
ac}(\mathbb{R}^{d})$, we define the $\mathsf{KL}$ divergence and the
(relative) Fisher information by

$\mathsf{KL}(\mu\mathbin{\|}\nu)\coloneqq\operatorname{\mathbb{E}}_{\mu}\bigl{[}\log\frac{\mu}{\nu}\bigr{]}\qquad\text{and}\qquad\mathsf{FI}(\mu\mathbin{\|}\nu)\coloneqq\operatorname{\mathbb{E}}_{\mu}\bigl{[}\bigl{\lVert}\nabla\log\frac{\mu}{\nu}\bigr{\rVert}^{2}\bigr{]}\,,$

with the convention
$\mathsf{KL}(\mu\mathbin{\|}\nu)=\mathsf{FI}(\mu\mathbin{\|}\nu)=\infty$
whenever $\mu\not\ll\nu$.

We recall the definition of the log-Sobolev inequality, which is used both for
propagation of chaos arguments as well as mixing time bounds.

######  Definition 2 (Log-Sobolev Inequality).

A measure $\pi$ satisfies a log-Sobolev inequality with parameter
$C_{\mathsf{LSI}}$ if for all $\mu\in\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$,

$\displaystyle\mathsf{KL}(\mu\mathbin{\|}\pi)\leq\frac{C_{\mathsf{LSI}}}{2}\,\mathsf{FI}(\mu\mathbin{\|}\pi)\,.$
(LSI)

When $\log(1/\pi)$ is $\alpha$-uniformly convex for $\alpha>0$, it follows
from the Bakry–Émery condition that $\pi$ satisfies (LSI) with constant
$C_{\mathsf{LSI}}\leq 1/\alpha$ [BGL14, Proposition 5.7.1].

We can also define the $p$-Wasserstein distance $\mathcal{W}_{p}(\mu,\pi)$,
$p\geq 1$, between $\mu,\pi$ as

$\displaystyle\mathcal{W}_{p}^{p}(\mu,\pi)=\inf_{\gamma\in\Gamma(\mu,\pi)}\int\lVert
x-y\rVert^{p}\,\gamma(\mathrm{d}x,\mathrm{d}y),$

where $\Gamma(\mu,\pi)$ is the set of all joint probability measures on
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $\mu,\pi$ respectively.

Lastly, we recall that the celebrated Otto calculus interprets the space
$\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$, equipped with the $\mathcal{W}_{2}$
metric, as a formal Riemannian manifold [Ott01]. In particular, the
Wasserstein gradient of a functional $\mathcal{L}:\mathcal{P}_{2,\rm
ac}(\mathbb{R}^{d})\to\mathbb{R}\cup\\{\infty\\}$ is given as
$\nabla_{\mathcal{W}_{2}}\mathcal{L}=\nabla\delta\mathcal{L}$. Here,
$\delta\mathcal{L}$ is the first variation defined as follows: for all
$\nu_{0},\nu_{1}\in\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$,
$\delta\mathcal{L}(\nu_{0}):\mathbb{R}^{d}\to\mathbb{R}$ satisfies

$\displaystyle\lim_{t\searrow
0}\frac{\mathcal{L}((1-t)\,\nu_{0}+t\,\nu_{1})-\mathcal{L}(\nu_{0})}{t}=\langle\delta\mathcal{L}(\nu_{0}),\nu_{1}-\nu_{0}\rangle\coloneqq\int\delta\mathcal{L}(\nu_{0})\,\mathrm{d}(\nu_{1}-\nu_{0})\,.$

The first variation is defined up to an additive constant, but the Wasserstein
gradient is unambiguous. See [AGS08] for a rigorous development. As a
shorthand, we will write
$\delta\mathcal{L}(\nu_{0},x)\coloneqq\delta\mathcal{L}(\nu_{0})(x)$ and
similarly
$\nabla_{\mathcal{W}_{2}}\mathcal{L}(\nu_{0},x)\coloneqq\nabla_{\mathcal{W}_{2}}\mathcal{L}(\nu_{0})(x)$.

###  2.1 SDE Systems and Their Stationary Distributions

####  2.1.1 The Pairwise McKean–Vlasov Setting

In the formalism introduced in the previous section, we note that
($\mathsf{pMV}$) can be interpreted as Wasserstein gradient flow for
($\mathsf{pE}$). In this paper, we refer to ($\mathsf{pMV}$) as the _pairwise
McKean–Vlasov_ process. As noted in the introduction, it has the stationary
distribution (1.1) which minimizes ($\mathsf{pE}$). Recall also that the
equation ($\mathsf{pMV}$) is the mean-field limit of the finite-particle
system ($\mathsf{pMV}_{N}$). This $N$-particle system has the following
stationary distribution: for
$x^{1:N}=[x^{1},\dotsc,x^{N}]\in\mathbb{R}^{d\times N}$,

$\displaystyle\mu^{1:N}(x^{1:N})\propto\exp\Bigl{(}-\frac{2}{\sigma^{2}}\sum_{i\in[N]}V(x^{i})-\frac{1}{\sigma^{2}\,(N-1)}\sum_{i\in[N]}\sum_{j\in[N]\setminus
i}W(x^{i}-x^{j})\Bigr{)}\,.$ (2.1)

The system ($\mathsf{pMV}_{N}$) can be viewed as an approximation to
($\mathsf{pMV}$), with the expectation term in the drift replaced by an
empirical average. Note that the measure $\mu^{1:N}$ is
exchangeable.111Exchangeability refers to the property that the law of
$[x^{1},\ldots,x^{N}]$ equals the law of
$[x^{\sigma(1)},\ldots,x^{\sigma(N)}]$ for any permutation $\sigma$ of
$\\{1,\ldots,N\\}$. While the standard approach is to apply an Euler–Maruyama
discretization to ($\mathsf{pMV}_{N}$) in order to sample from
($\mathsf{pMV}$), our perspective is to write more sophisticated samplers for
$\mu^{1:N}$ directly. Indeed, unlike (1.1), the finite-particle stationary
distribution (2.1) is explicit and amenable to sampling methods.

####  2.1.2 The General McKean–Vlasov Setting

More generally, we consider the functional ($\mathsf{gE}$) where $\mathcal{F}$
is of the form
$\mathcal{F}(\mu)=\mathcal{F}_{0}(\mu)+\frac{\lambda}{2}\int\lVert\cdot\rVert^{2}\,\mathrm{d}\mu$
with $\lambda\geq 0$. The second term acts as regularization and is common in
the literature [FW23, SWN23]. We can describe its Wasserstein gradient flow as
the marginal law of a particle trajectory satisfying the following SDE, which
we call the _general McKean–Vlasov equation_ :

$\displaystyle\mathrm{d}X_{t}=\\{-\nabla_{\mathcal{W}_{2}}\mathcal{F}_{0}(\pi_{t},X_{t})-\lambda
X_{t}\\}\,\mathrm{d}t+\sigma\,\mathrm{d}B_{t}\,,$ ($\mathsf{gMV}$)

where $\pi_{t}=\operatorname{\mathsf{law}}(X_{t})$, and $\\{B_{t}\\}_{t\geq
0}$ is a standard Brownian motion on $\mathbb{R}^{d}$. The stationary
distribution $\pi$ of ($\mathsf{gMV}$), and its linearization $\pi_{\mu}$
around a measure $\mu\in\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$, satisfy the
following equations:

$\displaystyle\\!\\!\pi(x)\propto\exp\Bigl{(}-\frac{2}{\sigma^{2}}\,\delta\mathcal{F}_{0}(\pi,x)-\frac{\lambda\,\lVert
x\rVert^{2}}{\sigma^{2}}\Bigr{)}\>\;\text{and}\>\;\pi_{\mu}(x)\propto\exp\Bigl{(}-\frac{2}{\sigma^{2}}\,\delta\mathcal{F}_{0}(\mu,x)-\frac{\lambda\,\lVert
x\rVert^{2}}{\sigma^{2}}\Bigr{)}\,.$ (2.2)

The latter is called the proximal Gibbs distribution with respect to $\mu$.
The general dynamics corresponds to the mean-field limit of the following
finite-particle system described by an $N$-tuple of stochastic processes
$\\{X_{t}^{1:N}\\}_{t\geq 0}\coloneqq\\{(X^{1}_{t},\dotsc,X^{N}_{t})\\}_{t\geq
0}$:

$\displaystyle\mathrm{d}X_{t}^{i}=\\{-\nabla_{\mathcal{W}_{2}}\mathcal{F}_{0}(\rho_{X_{t}^{1:N}},X_{t}^{i})-\lambda
X_{t}^{i}\\}\,\mathrm{d}t+\sigma\,\mathrm{d}B_{t}^{i}\,,$ ($\mathsf{gMV}_{N}$)

and $\rho_{x^{1:N}}=\frac{1}{N}\sum_{i=1}^{N}\delta_{x^{i}}$ is the empirical
measure of the particle system. The stationary distribution for
($\mathsf{gMV}_{N}$) is given as follows [CRW22, (2.16)]: for
$x^{1:N}=[x^{1},\dotsc,x^{N}]\in\mathbb{R}^{d\times N}$,

$\displaystyle\mu^{1:N}(x^{1:N})\propto\exp\Bigl{(}-\frac{2N}{\sigma^{2}}\,\mathcal{F}_{0}(\rho_{x^{1:N}})-\frac{\lambda}{\sigma^{2}}\,\lVert
x^{1:N}\rVert^{2}\Bigr{)}\,.$ (2.3)

One can show that
$\nabla_{x^{i}}\mathcal{F}_{0}(\rho_{x^{1:N}})=\frac{1}{N}\,\nabla_{\mathcal{W}_{2}}\mathcal{F}_{0}(\rho_{x^{1:N}},x^{i})$,
and hence ($\mathsf{gMV}_{N}$) is simply the Langevin diffusion corresponding
to stationary measure (2.3). Moreover, when $\lambda=0$ and choosing
$\mathcal{F}_{0}(\mu)=\int V(x)\,\mu(\mathrm{d}x)+\iint
W(x-y)\,\mu(\mathrm{d}x)\,\mu(\mathrm{d}y)$, then the equations
($\mathsf{gMV}$), (2.2), ($\mathsf{gMV}_{N}$), and (2.3) reduce to
($\mathsf{pMV}$), (1.1), ($\mathsf{pMV}_{N}$), and (2.1), respectively.

##  3 Technical Ingredients

Our general approach for sampling from the stationary distribution $\pi$ in
either (1.1) or (2.2) is to directly apply an off-the-shelf sampler for the
finite-particle stationary distribution $\mu^{1:N}$. The theoretical
guarantees for this procedure require two main ingredients: (1) control of the
“bias”—i.e., the error incurred by approximating $\pi$ by the $1$-particle
marginal of $\mu^{1:N}$—and (2) verification of isoperimetric properties which
allow for fast sampling from the measure $\mu^{1:N}$.

###  3.1 Bias Control via Uniform-in-Time Propagation of Chaos

In this section, we focus on the first ingredient, namely, obtaining control
of the bias via uniform-in-time propagation of chaos results. Proofs for this
section are given in §A.

####  3.1.1 Pairwise McKean–Vlasov Setting

We first consider the pairwise McKean–Vlasov setting described in §2.1.1. Our
first propagation of chaos result uses the following three assumptions.

######  Assumption 1.

The potentials $V,W$ are $\beta_{V},\beta_{W}$-smooth respectively.

######  Assumption 2.

The distribution $\pi$ satisfies (LSI) with parameter $C_{\mathsf{LSI}}(\pi)$.

######  Assumption 3.

The ratio
$\rho\coloneqq\nicefrac{{\sigma^{4}}}{{8\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)}}$
is at least $3$.

Remark. Note that from (1.1), we typically would expect
$C_{\mathsf{LSI}}^{2}(\pi)$ to also scale as $\sigma^{4}$ (e.g., in the case
when $V$ and $W$ are $\alpha$-uniformly convex for $\alpha>0$). Therefore,
Assumption 3 is typically invariant to the scaling of $\sigma$ and can be
satisfied even for $\sigma\searrow 0$. Under these assumptions, we obtain a
sharp propagation of chaos result via a similar argument as [Lac22, LL23]. We
note that the former is more permissive regarding the constant in Assumption 3
as compared to this work.

######  Theorem 3 (Sharp Propagation of Chaos).

Under Assumptions 1, 2 and 3, for any $N\geq 100$ and $k\in[N]$, it holds that
$\mathsf{KL}(\mu^{1:k}\mathbin{\|}\pi^{\otimes
k})=\widetilde{\mathcal{O}}(dk^{2}/N^{2})$. Thus,
$\mathsf{KL}(\mu^{1:k}\mathbin{\|}\pi^{\otimes k})<\varepsilon^{2}$ if

$\displaystyle N\geq
100\vee\widetilde{\Omega}\bigl{(}k\sqrt{d}\,\varepsilon^{-1}\bigr{)}\,.$ (3.1)

We note that the rate in Theorem 3 is sharp; see the Gaussian case in Example
10. A condition such as Assumption 3 is in general necessary, since otherwise
the minimizer of ($\mathsf{pE}$) may not even be unique [[, see the example
and discussion in]]lacker2023sharp. However, it can be restrictive, as it
requires the interaction to be sufficiently weak. With the following convexity
assumption, we can obtain a propagation of chaos result without Assumption 3.

######  Assumption 4.

The potentials $V,W$ are $\alpha_{V},\alpha_{W}$-uniformly convex with
$\alpha_{V}+{\alpha_{W}^{-}}>0$. Here,
${\alpha_{W}^{-}}\coloneqq\alpha_{W}\wedge 0$ denotes the negative part of
$\alpha_{W}$.

The following weaker result consists of two parts. The first, a Wasserstein
propagation of chaos result, is based on [Szn91]. The second, building on the
first, is a uniform-in-time entropic propagation of chaos bound following from
a Fisher information bound. The arguments are similar to those in [Mal01,
Mal03], albeit simplified (since we work at stationarity) and presented here
with explicit constants.

######  Theorem 4 (Weak Propagation of Chaos).

Under Assumptions 1 and 4, for any
$N\geq\frac{\alpha_{V}-{\alpha_{W}^{-}}}{\alpha_{V}+{\alpha_{W}^{-}}}\vee 2,$
if we denote $\alpha\coloneqq\alpha_{V}+{\alpha_{W}^{-}}$, then

$\displaystyle\mathcal{W}_{2}^{2}(\mu^{1:k},\pi^{\otimes k})$
$\displaystyle\leq\frac{4\beta_{W}^{2}\sigma^{2}d}{\alpha^{3}}\,\frac{k}{N}\,,$
(3.2) $\displaystyle\mathsf{KL}(\mu^{1:k}\mathbin{\|}\pi^{\otimes
k})\leq\frac{\sigma^{2}}{4\alpha}\operatorname{\mathsf{FI}}(\mu^{1:k}\mathbin{\|}\pi^{\otimes
k})$
$\displaystyle\leq\frac{132\beta_{W}^{2}\,{(\beta_{V}+\beta_{W})}^{2}\,d}{\alpha^{4}}\,\frac{k}{N}\,.$
(3.3)

####  3.1.2 General McKean–Vlasov Setting

In the more general case where we aim to minimize ($\mathsf{gE}$) for a
generic functional $\mathcal{F}$ of the form
$\mathcal{F}(\mu)=\mathcal{F}_{0}(\mu)+\frac{\lambda}{2}\int\lVert\cdot\rVert^{2}\,\mathrm{d}\mu$,
we impose the following assumptions. They can be largely seen as
generalizations of the conditions for the pairwise case, and they are
inherited from [CRW22, SWN23]. There is an additional convexity condition
(Assumption 5), which in the pairwise McKean–Vlasov setting amounts to
positive semidefiniteness of the kernel $(x,y)\mapsto W(x-y)$ on
$\mathbb{R}^{d}\times\mathbb{R}^{d}$; thus, in general, the following
assumptions are incomparable with the ones in §3.1.1.

######  Assumption 5.

The functional $\mathcal{F}_{0}$ is convex in the usual sense. For all
$\nu_{0},\nu_{1}\in\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$, $t\in[0,1]$,

$\displaystyle\mathcal{F}_{0}((1-t)\,\nu_{0}+t\,\nu_{1})\leq(1-t)\,\mathcal{F}_{0}(\nu_{0})+t\,\mathcal{F}_{0}(\nu_{1})\,.$

######  Assumption 6.

The functional $\mathcal{F}$ is smooth in the sense that for all
$x,y\in\mathbb{R}^{d}$, $\nu,\nu^{\prime}\in\mathcal{P}_{2,\rm
ac}(\mathbb{R}^{d})$, there is a uniform constant $\beta$ such that

$\displaystyle\lVert\nabla_{\mathcal{W}_{2}}\mathcal{F}(\nu,x)-\nabla_{\mathcal{W}_{2}}\mathcal{F}(\nu^{\prime},y)\rVert\leq\beta\,(\lVert
x-y\rVert+\mathcal{W}_{1}(\nu,\nu^{\prime}))\,.$

######  Assumption 7.

The proximal Gibbs measures satisfy (LSI) with a uniform constant: namely, it
holds that $C_{\mathsf{LSI}}(\pi)\vee\sup_{\mu\in\mathcal{P}_{2,\rm
ac}(\mathbb{R}^{d})}C_{\mathsf{LSI}}(\pi_{\mu})\leq\overline{C}_{\mathsf{LSI}}$.

Remark. These assumptions taken together cover settings not covered in the
preceding sections, including optimization of two-layer neural networks. See
[CRW22, Remark 3.1] and §4.2.

Under these assumptions, we can derive an entropic propagation of chaos bound
by following the proof of [CRW22]. Through a tighter analysis, we are able to
reduce the dependence on the condition number
$\kappa\coloneqq\overline{C}_{\mathsf{LSI}}\beta/\sigma^{2}$ from $\kappa^{2}$
to $\kappa$.

######  Theorem 5 (Propagation of Chaos for General Functionals).

Under Assumptions 5, 6, and 7, for
$N\geq\nicefrac{{160\beta\overline{C}_{\mathsf{LSI}}}}{{\sigma^{2}}},$ we have

$\frac{1}{2\overline{C}_{\mathsf{LSI}}}\mathcal{W}_{2}^{2}(\mu^{1:k},\pi^{\otimes
k})\leq\operatorname{\mathsf{KL}}(\mu^{1:k}\mathbin{\|}\pi^{\otimes
k})\leq\frac{33\beta\overline{C}_{\mathsf{LSI}}dk}{\sigma^{2}N}\,.$

Among these assumptions, the hardest to verify is the uniform LSI of
Assumption 7. Following [SWN23], we introduce the following sufficient
condition for the validity of Assumption 7; see Lemma 21 for a more precise
statement.

######  Assumption 8.

There exists a uniform bound on the Wasserstein gradient of the interaction
term $\mathcal{F}_{0}$: for some constant $B<\infty$ and all
$\mu\in\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$, $x\in\mathbb{R}^{d}$,
$\lVert\nabla_{\mathcal{W}_{2}}\mathcal{F}_{0}(\mu,x)\rVert\leq B\,.$

######  Lemma 6 (Informal).

Assumptions 6 and 8 imply Assumption 7 with an explicit constant
$\overline{C}_{\mathsf{LSI}}$, given in terms of $d$, $B$, $\beta$, $\lambda$,
and $\sigma$.

###  3.2 Isoperimetric Properties of the Stationary Distributions

In this section, we verify the isoperimetric properties of $\pi,\mu^{1:N}$,
with proofs provided in §B.

####  3.2.1 Pairwise McKean–Vlasov Setting

If $V$, $W$ satisfy Assumptions 1 and 4 (i.e., $V$ and $W$ have bounded
Hessians), then the potential for (1.1), i.e., $\log(1/\pi)$, is
$\frac{2}{\sigma^{2}}\,(\alpha_{V}+\alpha_{W})$-convex and
$\frac{2}{\sigma^{2}}\,(\beta_{V}+\beta_{W})$-smooth. By the Bakry–Émery
condition, $\pi$ satisfies (LSI) with parameter
$C_{\mathsf{LSI}}(\pi)\leq\nicefrac{{\sigma^{2}}}{{2\,(\alpha_{V}+\alpha_{W})}}$.

Similarly, for the invariant measure $\mu^{1:N}$ in (2.1), we can prove the
following.

######  Lemma 7.

If $V$ and $W$ satisfy Assumption 1, then $\log(1/\mu^{1:N})$ is
$\frac{2}{\sigma^{2}}\,(\beta_{V}+\frac{N}{N-1}\,\beta_{W})$-smooth.

If $V$ and $W$ satisfy Assumption 4, then $\log(1/\mu^{1:N})$ is
$\frac{2}{\sigma^{2}}\,(\alpha_{V}+\frac{N}{N-1}\,{\alpha_{W}^{-}})$-convex.222Only
the negative part of $\alpha_{W}$ contributes to the strong log-concavity of
$\mu^{1:N}$. This is consistent with [Vil03, Theorem 5.15], which asserts that
when $\alpha_{W}>0$, the interaction energy $\mu\mapsto\iint
W(x-y)\,\mu(\mathrm{d}x)\,\mu(\mathrm{d}y)$ is $\alpha_{W}$-strongly
displacement convex over the subspace of probability measures with fixed mean,
but only weakly convex over the full Wasserstein space.

We now consider the non-log-concave case. It is standard in the sampling
literature that the assumption of (LSI) for the stationary distribution yields
mixing time guarantees. Since our strategy is to sample from (2.1), we
therefore seek an LSI for $\mu^{1:N}$, formalized as the following assumption.

######  Assumption 9.

The distribution $\mu^{1:N}$ satisfies (LSI) with parameter
$C_{\mathsf{LSI}}(\mu^{1:N})$.

In this section, we provide an easily verifiable condition, based on the
Holley–Stroock condition [HS87], for this assumption to hold with an
$N$-independent constant.

######  Assumption 10.

The potentials $V$ and $W$ can be decomposed as $V=V_{0}+V_{1}$ and
$W=W_{0}+W_{1}$ such $V_{0}$, $W_{0}$ satisfy Assumption 4 and
$\operatorname{\mathsf{osc}}(V_{1}),\operatorname{\mathsf{osc}}(W_{1})<\infty$,
where for a function $U:\mathbb{R}^{d}\to\mathbb{R}$ we define
$\operatorname{\mathsf{osc}}(U)\coloneqq\sup U-\inf U$.

Under this assumption, a careful application of the Holley–Stroock
perturbation principle yields the following lemma.

######  Lemma 8.

Under Assumption 10, $\pi,\mu^{1:N}$ satisfy (LSI) with parameters

$\displaystyle C_{\mathsf{LSI}}(\pi)$
$\displaystyle\leq\frac{\sigma^{2}}{2\,(\alpha_{V_{0}}+\alpha_{W_{0}})}\exp\Bigl{(}\frac{2}{\sigma^{2}}\,\bigl{(}\operatorname{\mathsf{osc}}(V_{1})+\operatorname{\mathsf{osc}}(W_{1})\bigr{)}\Bigr{)}\,$
$\displaystyle C_{\mathsf{LSI}}(\mu^{1:N})$
$\displaystyle\leq\frac{\sigma^{2}}{2\,(\alpha_{V_{0}}+\frac{N}{N-1}\,{\alpha_{W_{0}}^{-}})}\exp\Bigl{(}\frac{2}{\sigma^{2}}\,\bigl{(}\operatorname{\mathsf{osc}}(V_{1})+\operatorname{\mathsf{osc}}(W_{1})\bigr{)}\Bigr{)}\,.$

####  3.2.2 General McKean–Vlasov Setting

In the setting ($\mathsf{gE}$) with
$\mathcal{F}(\mu)=\mathcal{F}_{0}(\mu)+\frac{\lambda}{2}\int\lVert\cdot\rVert^{2}\,\mathrm{d}\mu$,
we verify that Assumption 8 yields (LSI) for $\mu^{1:N}$ with an
$N$-independent constant. See Lemma 22 for a more precise statement.

######  Lemma 9 (Informal).

In the mean-field Langevin setting of §2.1.2, suppose that Assumption 8 holds.
Then, Assumption 9 holds with $C_{\mathsf{LSI}}(\mu^{1:N})$ depending on $d$,
$B$, $\lambda$, and $\sigma$, but not on $N$.

##  4 Sampling from the Mean-Field Target

In this section, we present results for sampling from $\pi$. As outlined in
Algorithm 1, we use off-the-shelf log-concave samplers to sample from
$\mu^{1:N}$, during which we access the first-order333For our results
involving the proximal sampler, we also assume access to a proximal oracle for
simplicity. oracle for $\mu^{1:N}$ (i.e., an oracle for evaluation of
$\log\mu^{1:N}$ up to an additive constant, and for evaluation of
$\nabla\log\mu^{1:N}$). For $N$ sufficiently large, the first particle given
by Algorithm 1 is approximately distributed according to $\pi$: for
$\hat{\mu}^{1:N}$ the law of the output of the log-concave sampler and its
$1$-particle marginal distribution $\hat{\mu}^{1}$,

$\mathcal{W}_{2}(\hat{\mu}^{1},\pi)\leq\mathcal{W}_{2}(\hat{\mu}^{1},\mu^{1})+\mathcal{W}_{2}(\mu^{1},\pi)\leq\sqrt{\frac{1}{N}}\,\mathcal{W}_{2}(\hat{\mu}^{1:N},\mu^{1:N})+\mathcal{W}_{2}(\mu^{1},\pi)\,,$

where the inequality follows from exchangeability (Lemma 23). A similar
decomposition also holds for $\operatorname{\mathsf{KL}}$, although the
argument is more technical. We defer its presentation to §E.

Input: the number $N$ of total particles, a log-concave sampler
$\mathsf{LC{\text{-}}Sampler}$
Output: $k$ particles $\hat{X}^{1:k}$

1: Sample $\hat{X}^{1:N}\sim\hat{\mu}^{1:N}$ via
$\mathsf{LC{\text{-}}Sampler}$, so that $\hat{\mu}^{1:N}\approx\mu^{1:N}$,
e.g., in $\mathcal{W}_{2}$ or $\sqrt{\operatorname{\mathsf{KL}}}$.

2: Output the first $k$ particles $\hat{X}^{1:k}$.

Algorithm 1 Sampling from the Mean-Field Stationary Distribution

To bound the second term by $\varepsilon$, it suffices to choose $N$ according
to the propagation of chaos results in §3.1. Our results are summarized in
Table 1, in which we record the total number of oracle calls $M$ for
$\mu^{1:N}$ made by the sampler $M$ and the number of particles $N$ needed to
achieve $\varepsilon$ error in the desired metric, hiding polylogarithmic
factors. Note that in the pairwise McKean–Vlasov setting, each oracle call to
$\mu^{1:N}$ requires $N$ calls to an oracle for $V$, and $\binom{N}{2}$ calls
to an oracle for $W$.

Algorithm | “Metric” | Assumptions | $M$ | $N$
---|---|---|---|---
LMC | $\nicefrac{{\sqrt{\alpha}}}{{\sigma}}\,\mathcal{W}_{2}$ | 1, 2, 3, 9 | ${\kappa^{2}d}/{\varepsilon^{2}}$ | ${d^{1/2}}/{\varepsilon}$
MALA–PS | $\kappa d^{3/4}/\varepsilon^{1/2}$
ULMC–PS | ${\kappa^{3/2}d^{1/2}}/{\varepsilon}$
ULMC+ | $\nicefrac{{\sqrt{\alpha}}}{{\sigma}}\,\mathcal{W}_{2}$ |  1, 3, 4 | ${\kappa d^{1/3}}/{\varepsilon^{2/3}}$ | ${d^{1/2}}/{\varepsilon}$
LMC | $\sqrt{\operatorname{\mathsf{KL}}}$ | 1, 4 | ${\kappa^{2}d}/{\varepsilon^{2}}$ | ${\kappa^{4}d}/{\varepsilon^{2}}$
ULMC | ${\kappa^{3/2}d^{1/2}}/{\varepsilon}$
LMC | $\nicefrac{{\sqrt{\alpha}}}{{\sigma}}\,\mathcal{W}_{2}$ | ${\kappa d}/{\varepsilon^{2}}$ | ${\kappa^{2}d}/{\varepsilon^{2}}$
ULMC+ | ${\kappa d^{1/3}}/{\varepsilon^{2/3}}$
LMC | $\sqrt{\operatorname{\mathsf{KL}}}$ | 5, 6, 7, 9 | ${\kappa^{2}d}/{\varepsilon^{2}}$ | ${\kappa d}/{\varepsilon^{2}}$
ULMC–PS | $\kappa^{3/2}d^{1/2}/\varepsilon$
Table 1: In this table, we record $M$, the total number of oracle queries to
$\log\mu^{1:N}$ made by the log-concave sampler, and $N$, the number of
particles.

The algorithms in the table refer to: Langevin Monte Carlo (LMC); underdamped
Langevin Monte Carlo (ULMC); discretizations of the underdamped Langevin
diffusion via the randomized midpoint method [SL19] or the shifted ODE method
[FLO21] (ULMC+); and implementation of the proximal sampler [LST21, Che+22]
via the Metropolis-adjusted Langevin algorithm or via ULMC (MALA–PS and
ULMC–PS respectively). Note that LMC applied to sample from $\mu^{1:N}$ is
simply the Euler–Maruyama discretization of ($\mathsf{pMV}_{N}$), and likewise
ULMC is the algorithm considered in [FW23]. We refer to §E for proofs and
references.

To streamline the rates, we simplify the notation by defining
$\beta=\beta_{V}+\beta_{W}$ if Assumption 1 holds, otherwise we use the value
from Assumption 6. We let $\alpha=\alpha_{V}+{\alpha_{W}^{-}}$ under
Assumption 4,
$\alpha=\nicefrac{{\sigma^{2}}}{{2\max\\{C_{\mathsf{LSI}}(\mu^{1:N}),C_{\mathsf{LSI}}(\pi)\\}}}$
under Assumptions 2 and 9, and
$\alpha=\nicefrac{{\sigma^{2}}}{{2\max\\{C_{\mathsf{LSI}}(\mu^{1:N}),\overline{C}_{\mathsf{LSI}}\\}}}$
in the general McKean–Vlasov setting. Finally, we let
$\kappa\coloneqq\beta/\alpha$ denote the condition number. The additional
assumption $\kappa\leq\sqrt{d}/\varepsilon$ will be used to simplify some of
the rates.

In the following subsections, we discuss some of the results in greater
detail.

###  4.1 Pairwise McKean–Vlasov Setting

######  Example 10 (Gaussian Case).

Consider a quadratic confinement and interaction,

$V(x)=\frac{1}{2}\,x^{\mathsf{T}}Ax=\frac{1}{2}\,\lVert
x\rVert_{A}^{2}\,,\qquad W(x)=\frac{\lambda}{2}\,\lVert x\rVert^{2}\,,$

for some matrix $A\in\mathbb{R}^{d\times d}$ with $A\succ 0$, $\lambda\geq 0$.
The resulting stationary distributions can be calculated explicitly to be
Gaussians. We show in §C that for large $N$,
$\operatorname{\mathsf{KL}}(\mu^{1:k}\mathbin{\|}\pi^{\otimes
k})=\widetilde{\Theta}(\nicefrac{{dk^{2}}}{{N^{2}}}\,)$. This shows that the
rate in Theorem 3 is sharp.

######  Example 11 (Strongly Convex Case).

Consider the strongly convex case where
$\alpha=\alpha_{V}+{\alpha_{W}^{-}}>0$. The prior work [BS23] also considered
the problem of sampling from the mean-field stationary distribution $\pi$,
with $\sigma^{2}=2$. If we count the number of calls to a gradient oracle for
$V$, their complexity bound reads
$\widetilde{\mathcal{O}}(\nicefrac{{\kappa^{5/3}d^{4/3}}}{{\varepsilon^{8/3}}})$
to achieve
$\nicefrac{{\sqrt{\alpha}}}{{\sigma}}\,\mathcal{W}_{1}(\hat{\mu}^{1},\pi)\leq\varepsilon$.
We note that their assumptions are not strictly comparable to ours. They
require the interaction $W$ to be sufficiently weak, in the sense that
$\beta_{W}\lesssim\alpha$, which is similar444See eq. (2.24) therein; note
that they have a scaling factor of $\varepsilon$ in front of their interaction
term, so that our parameter $\beta_{W}$ is equivalent to their
$\varepsilon\tilde{L}$. to our Assumption 3; on the other hand, they only
assume $\alpha_{V}>0$, rather than $\alpha_{V}+{\alpha_{W}^{-}}>0$.
Nevertheless, we attempt to make some comparisons with their work below.

Without Assumption 3, ULMC+ achieves
$\nicefrac{{\sqrt{\alpha}}}{{\sigma}}\,\mathcal{W}_{2}(\hat{\mu}^{1},\pi)\leq\varepsilon$
with complexity
$\widetilde{\mathcal{O}}(\nicefrac{{\kappa^{3}d^{4/3}}}{{\varepsilon^{8/3}}})$,
which matches the guarantee of [BS23] up to the dependence on $\kappa$. We can
also obtain guarantees in $\sqrt{\operatorname{\mathsf{KL}}}$, at the cost of
an extra factor of $\kappa^{2}$.

With Assumption 3, MALA–PS has complexity
$\widetilde{\mathcal{O}}(\nicefrac{{\kappa d^{5/4}}}{{\varepsilon^{3/2}}})$
and ULMC+ has complexity $\widetilde{\mathcal{O}}(\nicefrac{{\kappa
d^{5/6}}}{{\varepsilon^{5/3}}})$, which improve substantially upon [BS23].

To summarize, in the strongly convex case, we have obtained numerous
improvements: (i) we can obtain results even without the weak interaction
condition (Assumption 3); (ii) when we assume the weak interaction condition,
we obtain improved complexities; (iii) our results hold in stronger metrics;
(iv) our approach is generic, allowing for the consideration of numerous
different samplers without needing to establish new propagation of chaos
results (by way of comparison, [BS23] developed a tailored propagation of
chaos argument for their non-linear Hamiltonian Monte Carlo algorithm).

######  Example 12 (Bounded Perturbations).

Both the results of [BS23] as well as our own allow for non-convex potentials,
albeit under different assumptions—[BS23] require strong convexity at
infinity, whereas we require (LSI) for the stationary measures $\mu^{1:N}$ and
$\pi$. In order to obtain sampling guarantees with low complexity, it is
important for the LSI constant of $\mu^{1:N}$ to be independent of $N$. We
have provided a sufficient condition for this to hold: $V$ and $W$ are bounded
perturbations of $V_{0}$ and $W_{0}$ respectively, where
$\alpha_{V_{0}}+{\alpha_{W_{0}}^{-}}>0$; see Lemma 8.

We also note that in this setting, both of our works require a weak
interaction condition. This is in general necessary in order to ensure
uniqueness of the mean-field stationary distribution, see the discussion in
§3.1.1.

###  4.2 General McKean–Vlasov Setting

######  Example 13 (General Functionals).

In the general setting, under Assumptions 5, 6, and 7, the work of [SWN23]
provided the first discretization bounds. They impose further assumptions and
their resulting complexity bound is rather complicated, but it reads roughly
$MN=\widetilde{\mathcal{O}}(\nicefrac{{\operatorname{poly}(\kappa)\,d^{2}}}{{\varepsilon^{4}}})$
for the discretization of ($\mathsf{gMV}_{N}$). Subsequently, [FW23] obtained
an improved complexity of
$MN=\widetilde{\mathcal{O}}(\nicefrac{{\kappa^{4}d^{3/2}}}{{\varepsilon^{3}}})$
via ULMC in the averaged $\mathsf{TV}$ distance. In comparison, we can improve
this complexity guarantee to
$\widetilde{\mathcal{O}}(\nicefrac{{\kappa^{5/2}d^{3/2}}}{{\varepsilon^{3}}})$,
and the guarantee even holds in $\sqrt{\operatorname{\mathsf{KL}}}$ if we
combine ULMC with the proximal sampler. It appears that we gain one factor of
$\sqrt{\kappa}$ through sharper discretization analysis (via [Zha+23], or via
the error analysis of the proximal sampler in [AC23]), and one factor of
$\kappa$ via a sharper propagation of chaos result (Theorem 5).

We also note that the result of [FW23] is based on a kinetic version of the
propagation of chaos argument from [Che+23], whereas our approach uses the
original “non-kinetic” argument from [CRW22] in the form of Theorem 5.

##### Application to Two-Layer Neural Networks.

Let us consider the problem of learning a two-layer neural network in the
mean-field regime. Let $f_{\theta}:\mathbb{R}^{d}\to\mathbb{R}$ be a function
parameterized by $\theta\in\mathbb{R}^{p}$, and for any probability measure
$\mu$ over $\mathbb{R}^{p}$, let $f_{\mu}\coloneqq\int
f_{\theta}\,\mu(\mathrm{d}\theta)$. For example, in a standard two-layer
neural network, we take $\theta=(a,w)\in\mathbb{R}\times\mathbb{R}^{d}$ and
$f_{a,w}(x)=a\,\mathsf{ReLU}(\langle w,x\rangle)$. When
$\mu=\frac{1}{m}\sum_{j=1}^{m}\delta_{(a_{j},w_{j})}$ is an empirical measure,
then $f_{\mu}$ is the function computed by a two-layer neural network with $m$
hidden neurons. In this formulation, however, we can take $\mu$ to be any
probability measure, corresponding to the _mean-field limit_ $m\to\infty$
[CB18, MMN18, Chi22, RV22, SS20].

Given a dataset $\\{(x_{i},y_{i})\\}_{i=1}^{n}$ in
$\mathbb{R}^{d}\times\mathbb{R}$ and a loss function
$\ell:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$, we can formulate neural
network training as the problem of minimizing the loss
$\mu\mapsto\sum_{i=1}^{n}\ell(f_{\mu}(x_{i}),y_{i})$. To place this within the
general McKean–Vlasov framework, we add two regularization terms: (1)
$\frac{\lambda}{2}\int\lVert\cdot\rVert^{2}\,\mathrm{d}\mu$ corresponds to
weight decay; and (2) $\frac{\sigma^{2}}{2}\int\log\mu\,\mathrm{d}\mu$ is
entropic regularization. We are now in the setting of §2.1.2, with
$\mathcal{F}_{0}(\mu)=\sum_{i=1}^{n}\ell(f_{\mu}(x_{i}),y_{i})$. To minimize
this energy, it is natural to consider the Euler–Maruyama discretization of
($\mathsf{gMV}_{N}$), which corresponds to learning the neural network via
noisy GD, and was considered in [SWN23]. Recent works [Che+23, FW23] also
considered the underdamped version of ($\mathsf{gMV}$) and its discretization.
Under the assumptions common to those works as well as our own, our results
yield improved algorithmic guarantees for this task (see Example 13).

Unfortunately, the assumptions used for the analysis of the general
McKean–Vlasov are restrictive and limit the applicability to neural network
training. For example, it suffices for $\ell$ to be convex in its first
argument (to satisfy Assumption 5), to have two bounded derivatives (w.r.t.
its first argument), and for $\theta\mapsto f_{\theta}(x_{i})$ to have two
bounded derivatives for each $x_{i}$. The last condition is satisfied, e.g.,
for $f_{\theta}(x)=\tanh(\langle\theta,x\rangle)$. For a genuinely two-layer
example, we can take $f_{\theta}(x)=\tanh(a)\tanh(\langle w,x\rangle)$ for
$\theta=(a,w)\in\mathbb{R}\times\mathbb{R}^{d}$. Under these conditions,
Assumptions 6 and 8 hold, which in turn furnish log-Sobolev inequalities via
Lemmas 6 and 9. In general, these LSI constants depend exponentially on
quantities such as the dimension $p$ of the parameter space, which is
unavoidable as global non-convex optimization is intractable in the worst
case. However, we note that these limitations are inherited from the prior
literature and not specific to our approach.

##  5 Conclusion

In this work, we propose a framework for obtaining sampling guarantees for the
minimizers of ($\mathsf{pE}$) and ($\mathsf{gE}$), based on decoupling the
problem into _(i)_ particle approximation via _propagation of chaos_ , and
_(ii)_ time-discretization via _log-concave sampling theory_. Our approach
leads to simpler proofs and improved guarantees compared to previous works,
and our results readily benefit from any improvements in either _(i)_ or
_(ii)_.

We conclude by listing some future directions of study. We believe there is
further room for improvement in the propagation of chaos results. For example,
can the sharp rate in Theorem 3 be extended to stronger metrics such as Rényi
divergences, as well as to situations when the weak interaction condition
(Assumption 3) fails, e.g., in the strongly displacement convex case or in the
setting of §3.1.2? For the sampling guarantees, the prior works [BS23, SWN23]
considered different settings, such as potentials satisfying convexity at
infinity or the use of stochastic gradients; these extensions are compatible
with our approach and could possibly lead to improvements in these cases, as
well as others. Finally, it would be interesting to develop further
applications of sampling from the mean-field distribution, e.g., applications
to neural network training under less stringent hypotheses.

## Acknowledgements

We thank Daniel Lacker for providing invaluable feedback and for pointing out
some important references.

YK was supported in part by NSF awards CCF-2007443 and CCF-2134105. MSZ was
supported by NSERC through the CGS-D program. SC was supported by the Eric and
Wendy Schmidt Fund at the Institute for Advanced Study. MAE was supported by
NSERC Grant [2019-06167] and CIFAR AI Chairs program at the Vector Institute.
MBL was supported by NSF grant DMS-2133806.

range

pages8 rangepages-1 rangepages54 rangepages29 rangepages18 rangepages88
rangepages36 rangepages-1 rangepages51 rangepages28 rangepages121
rangepages157 rangepages31 rangepages2 rangepages41 rangepages26 rangepages-1
rangepages32 rangepages-1 rangepages49 rangepages18 rangepages44 rangepages20
rangepages11 rangepages36 rangepages24 rangepages69 rangepages17 rangepages37
rangepages21 rangepages38 rangepages32 rangepages58 rangepages35 rangepages24
rangepages21 rangepages5 rangepages-1 rangepages54 rangepages17 rangepages17
rangepages74 rangepages40 rangepages49 rangepages47 rangepages28 rangepages87
rangepages49 rangepages235 rangepages-1 rangepages935 rangepages63
rangepages16 rangepages36

## References

  * [AC23] Jason M. Altschuler and Sinho Chewi  “Faster high-accuracy log-concave sampling via algorithmic warm starts”  In _2023 IEEE 64th Annual Symposium on Foundations of Computer Science (FOCS)_ , 2023, pp. 2169–2176
  * [AGS08] Luigi Ambrosio, Nicola Gigli and Giuseppe Savaré  “Gradient flows in metric spaces and in the space of probability measures”, Lectures in Mathematics ETH Zürich  Birkhäuser Verlag, Basel, 2008, pp. x+334
  * [AK02] Fabio Antonelli and Arturo Kohatsu-Higa  “Rate of convergence of a particle method to the solution of the McKean–Vlasov equation”  In _The Annals of Applied Probability_ 12.2 Institute of Mathematical Statistics, 2002, pp. 423–476
  * [BGL14] Dominique Bakry, Ivan Gentil and Michel Ledoux  “Analysis and geometry of Markov diffusion operators”  Springer, 2014
  * [BH22] Jianhai Bao and Xing Huang  “Approximations of McKean–Vlasov stochastic differential equations with irregular coefficients”  In _Journal of Theoretical Probability_ Springer, 2022, pp. 1–29
  * [BGM10] François Bolley, Arnaud Guillin and Florent Malrieu  “Trend to equilibrium and particle approximation for a weakly self consistent Vlasov–Fokker–Planck equation”  In _ESAIM: Mathematical Modelling and Numerical Analysis_ 44.5 EDP Sciences, 2010, pp. 867–884
  * [Bol72] Ludwig Boltzmann  “Further studies on the heat balance among gas molecules”  In _History of Modern Physical Sciences_ 1, 1872, pp. 262–349
  * [BT97] Mireille Bossy and Denis Talay  “A stochastic particle method for the McKean–Vlasov and the Burgers equation”  In _Mathematics of Computation_ 66.217, 1997, pp. 157–192
  * [BS23] Nawaf Bou-Rabee and Katharina Schuh  “Nonlinear Hamiltonian Monte Carlo & its particle approximation”  In _arXiv preprint arXiv:2308.11491_ , 2023
  * [BLM13] Stéphane Boucheron, Gábor Lugosi and Pascal Massart  “Concentration inequalities” A nonasymptotic theory of independence, With a foreword by Michel Ledoux  Oxford University Press, Oxford, 2013, pp. x+481
  * [BJW23] Didier Bresch, Pierre-Emmanuel Jabin and Zhenfu Wang  “Mean field limit and quantitative estimates with singular attractive kernels”  In _Duke Mathematical Journal_ 172.13 Duke University Press, 2023, pp. 2591–2641
  * [CD18] René Carmona and François Delarue  “Probabilistic theory of mean field games with applications I–II”  Springer, 2018
  * [CG22] Patrick Cattiaux and Arnaud Guillin  “Functional inequalities for perturbed measures with applications to log-concave measures and to some Bayesian problems”  In _Bernoulli_ 28.4 Bernoulli Society for Mathematical StatisticsProbability, 2022, pp. 2294–2321
  * [CD22] Louis-Pierre Chaintron and Antoine Diez  “Propagation of chaos: a review of models, methods and applications. I. Models and methods”  In _Kinet. Relat. Models_ 15.6, 2022, pp. 895–1015
  * [CD22a] Louis-Pierre Chaintron and Antoine Diez  “Propagation of chaos: a review of models, methods and applications. II. Applications”  In _Kinet. Relat. Models_ 15.6, 2022, pp. 1017–1173
  * [Che+23] Fan Chen, Yiqing Lin, Zhenjie Ren and Songbo Wang  “Uniform-in-time propagation of chaos for kinetic mean field Langevin dynamics”  In _arXiv preprint arXiv:2307.02168_ , 2023
  * [CRW22] Fan Chen, Zhenjie Ren and Songbo Wang  “Uniform-in-time propagation of chaos for mean field Langevin dynamics”  In _arXiv preprint arXiv:2212.03050_ , 2022
  * [Che+22] Yongxin Chen, Sinho Chewi, Adil Salim and Andre Wibisono  “Improved analysis for a proximal algorithm for sampling”  In _Proceedings of Thirty Fifth Conference on Learning Theory (COLT)_ 178, Proceedings of Machine Learning Research  PMLR, 2022, pp. 2984–3014
  * [Che23] Sinho Chewi  “Log-concave sampling” Book draft available at https://chewisinho.github.io., 2023
  * [Che+22a] Sinho Chewi et al.  “Analysis of Langevin Monte Carlo from Poincaré to log-Sobolev”  In _Proceedings of Thirty Fifth Conference on Learning Theory (COLT)_ 178, Proceedings of Machine Learning Research  PMLR, 2022, pp. 1–2
  * [Che+21] Sinho Chewi et al.  “Optimal dimension dependence of the Metropolis-adjusted Langevin algorithm”  In _Conference on Learning Theory (COLT)_ , 2021, pp. 1260–1300  PMLR
  * [Chi22] Lénaı̈c Chizat  “Mean-field Langevin dynamics: exponential convergence and annealing”  In _Transactions on Machine Learning Research_ , 2022
  * [CB18] Lénaı̈c Chizat and Francis Bach  “On the global convergence of gradient descent for over-parameterized models using optimal transport”  In _Advances in Neural Information Processing Systems (NeurIPS)_ 31, 2018
  * [Cla+23] Julien Claisse, Giovanni Conforti, Zhenjie Ren and Songbo Wang  “Mean field optimization problem regularized by Fisher information”  In _arXiv preprint 2302.05938_ , 2023
  * [Csi84] Imre Csiszár  “Sanov property, generalized I-projection and a conditional limit theorem”  In _The Annals of Probability_ JSTOR, 1984, pp. 768–793
  * [DKR22] Arnak S. Dalalyan, Avetik Karagulyan and Lionel Riou-Durand  “Bounding the error of discretized Langevin algorithms for non-strongly log-concave targets”  In _J. Mach. Learn. Res._ 23, 2022, pp. Paper No. 23538
  * [Dia+23] Michael Z. Diao, Krishnakumar Balasubramanian, Sinho Chewi and Adil Salim  “Forward-backward Gaussian variational inference via JKO in the Bures–Wasserstein space”  In _Proceedings of the 40th International Conference on Machine Learning_ 202, Proceedings of Machine Learning Research  PMLR, 2023, pp. 7960–7991
  * [DMM19] Alain Durmus, Szymon Majewski and Błażej Miasojedow  “Analysis of Langevin Monte Carlo via convex optimization”  In _J. Mach. Learn. Res._ 20, 2019, pp. Paper No. 7346
  * [FYC23] Jiaojiao Fan, Bo Yuan and Yongxin Chen  “Improved dimension dependence of a proximal algorithm for sampling”  In _Proceedings of Thirty Sixth Conference on Learning Theory (COLT)_ 195, Proceedings of Machine Learning Research  PMLR, 2023, pp. 1473–1521
  * [FLO21] James Foster, Terry Lyons and Harald Oberhauser  “The shifted ODE method for underdamped Langevin MCMC”  In _arXiv preprint 2101.03446_ , 2021
  * [FW23] Qiang Fu and Ashia Wilson  “Mean-field underdamped Langevin dynamics and its space-time discretization”  In _arXiv preprint arXiv:2312.16360_ , 2023
  * [Fun84] Tadahisa Funaki  “A certain class of diffusion processes associated with nonlinear parabolic equations”  In _Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete_ 67.3 Springer, 1984, pp. 331–348
  * [GLM22] Arnaud Guillin, Pierre Le Bris and Pierre Monmarché  “Convergence rates for the Vlasov–Fokker–Planck equation and uniform in time propagation of chaos in non convex cases”  In _Electronic Journal of Probability_ 27 The Institute of Mathematical Statisticsthe Bernoulli Society, 2022, pp. 1–44
  * [GM21] Arnaud Guillin and Pierre Monmarché  “Uniform long-time and propagation of chaos estimates for mean field kinetic particles in non-convex landscapes”  In _Journal of Statistical Physics_ 185 Springer, 2021, pp. 1–20
  * [HBE20] Ye He, Krishnakumar Balasubramanian and Murat A. Erdogdu  “On the ergodicity, bias and asymptotic normality of randomized midpoint sampling method”  In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 7366–7376
  * [HR23] Elias Hess-Childs and Keefer Rowan  “Higher-order propagation of chaos in ${L}^{2}$ for interacting diffusions”  In _arXiv preprint arXiv:2310.09654_ , 2023
  * [HS87] Richard Holley and Daniel Stroock  “Logarithmic Sobolev inequalities and stochastic Ising models”  In _J. Statist. Phys._ 46.5-6, 1987, pp. 1159–1194
  * [JW17] Pierre-Emmanuel Jabin and Zhenfu Wang  “Mean field limit for stochastic particle systems”  In _Active Particles, Volume 1: Advances in Theory, Models, and Applications_ Springer, 2017, pp. 379–402
  * [JW18] Pierre-Emmanuel Jabin and Zhenfu Wang  “Quantitative estimates of propagation of chaos for stochastic systems with $W^{-1,\infty}$ kernels”  In _Invent. Math._ 214.1, 2018, pp. 523–591
  * [JCP23] Yiheng Jiang, Sinho Chewi and Aram-Alexandre Pooladian  “Algorithms for mean-field variational inference via polyhedral optimization in the Wasserstein space”  In _arXiv preprint 2312.02849_ , 2023
  * [JKO98] Richard Jordan, David Kinderlehrer and Felix Otto  “The variational formulation of the Fokker–Planck equation”  In _SIAM J. Math. Anal._ 29.1, 1998, pp. 1–17
  * [Kuh+19] Daniel Kuhn, Peyman M. Esfahani, Viet Anh Nguyen and Soroosh Shafieezadeh-Abadeh  “Wasserstein distributionally robust optimization: theory and applications in machine learning”  In _Operations Research & Management Science in the Age of Analytics_, 2019, pp. 130–166
  * [Lac21] Daniel Lacker  “Hierarchies, entropy, and quantitative propagation of chaos for mean field diffusions”  In _arXiv preprint arXiv:2105.02983_ , 2021
  * [Lac22] Daniel Lacker  “Quantitative approximate independence for continuous mean field Gibbs measures”  In _Electronic Journal of Probability_ 27 The Institute of Mathematical Statisticsthe Bernoulli Society, 2022, pp. 1–21
  * [Lac23] Daniel Lacker  “Independent projections of diffusions: gradient flows for variational inference and optimal mean field approximations”  In _arXiv preprint 2309.13332_ , 2023
  * [LL23] Daniel Lacker and Luc Le Flem  “Sharp uniform-in-time propagation of chaos”  In _Probability Theory and Related Fields_ Springer, 2023, pp. 1–38
  * [Lam+22] Marc Lambert et al.  “Variational inference via Wasserstein gradient flows”  In _Advances in Neural Information Processing Systems_ , 2022
  * [LL07] Jean-Michel Lasry and Pierre-Louis Lions  “Mean field games”  In _Jpn. J. Math._ 2.1, 2007, pp. 229–260
  * [LST21] Yin Tat Lee, Ruoqi Shen and Kevin Tian  “Structured logconcave sampling with a restricted Gaussian oracle”  In _Proceedings of Thirty Fourth Conference on Learning Theory_ 134, Proceedings of Machine Learning Research  PMLR, 2021, pp. 2993–3050
  * [Li+23] Yun Li et al.  “Strong convergence of Euler–Maruyama schemes for McKean–Vlasov stochastic differential equations under local Lipschitz conditions of state variables”  In _IMA Journal of Numerical Analysis_ 43.2 Oxford University Press, 2023, pp. 1001–1035
  * [LW16] Qiang Liu and Dilin Wang  “Stein variational gradient descent: a general purpose Bayesian inference algorithm”  In _Advances in Neural Information Processing Systems_ 29 Curran Associates, Inc., 2016
  * [Mal01] Florent Malrieu  “Logarithmic Sobolev inequalities for some nonlinear PDE’s”  In _Stochastic Process. Appl._ 95.1, 2001, pp. 109–132
  * [Mal03] Florent Malrieu  “Convergence to equilibrium for granular media equations and their Euler schemes”  In _The Annals of Applied Probability_ 13.2 Institute of Mathematical Statistics, 2003, pp. 540–560
  * [McK66] Henry P. McKean  “A class of Markov processes associated with nonlinear parabolic equations”  In _Proc. Nat. Acad. Sci. U.S.A._ 56, 1966, pp. 1907–1911
  * [MMN18] Song Mei, Andrea Montanari and Phan-Minh Nguyen  “A mean field view of the landscape of two-layer neural networks”  In _Proceedings of the National Academy of Sciences_ 115.33 National Acad Sciences, 2018, pp. E7665–E7671
  * [Mél96] Sylvie Méléard  “Asymptotic behaviour of some interacting particle systems; McKean–Vlasov and Boltzmann models”  In _Probabilistic Models for Nonlinear Partial Differential Equations: Lectures given at the 1st Session of the Centro Internazionale Matematico Estivo (C.I.M.E.) held in Montecatini Terme, Italy, May 22–30, 1995_ Berlin, Heidelberg: Springer Berlin Heidelberg, 1996, pp. 42–95
  * [Mon17] Pierre Monmarché  “Long-time behaviour and propagation of chaos for mean field kinetic particles”  In _Stochastic Processes and their Applications_ 127.6 Elsevier, 2017, pp. 1721–1737
  * [MRW24] Pierre Monmarché, Zhenjie Ren and Songbo Wang  “Time-uniform log-Sobolev inequalities and applications to propagation of chaos”  In _arXiv preprint arXiv:2401.07966_ , 2024
  * [NWS22] Atsushi Nitanda, Denny Wu and Taiji Suzuki  “Convex analysis of the mean field Langevin dynamics”  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , 2022, pp. 9741–9757  PMLR
  * [Ott01] Felix Otto  “The geometry of dissipative evolution equations: the porous medium equation”  In _Comm. Partial Differential Equations_ 26.1-2, 2001, pp. 101–174
  * [OV00] Felix Otto and Cédric Villani  “Generalization of an inequality by Talagrand and links with the logarithmic Sobolev inequality”  In _Journal of Functional Analysis_ 173.2 Elsevier, 2000, pp. 361–400
  * [RES22] Gonçalo Reis, Stefan Engelhardt and Greig Smith  “Simulation of McKean–Vlasov SDEs with super-linear growth”  In _IMA Journal of Numerical Analysis_ 42.1 Oxford University Press, 2022, pp. 874–922
  * [RV22] Grant M. Rotskoff and Eric Vanden-Eijnden  “Trainability and accuracy of artificial neural networks: an interacting particle system approach”  In _Comm. Pure Appl. Math._ 75.9, 2022, pp. 1889–1935
  * [SL19] Ruoqi Shen and Yin Tat Lee  “The randomized midpoint method for log-concave sampling”  In _Advances in Neural Information Processing Systems (NeurIPS)_ 32, 2019
  * [SS20] Justin Sirignano and Konstantinos Spiliopoulos  “Mean field analysis of neural networks: A law of large numbers”  In _SIAM Journal on Applied Mathematics_ 80.2 SIAM, 2020, pp. 725–752
  * [SNW22] Taiji Suzuki, Atsushi Nitanda and Denny Wu  “Uniform-in-time propagation of chaos for the mean-field gradient Langevin dynamics”  In _The Eleventh International Conference on Learning Representations (ICLR)_ , 2022
  * [SWN23] Taiji Suzuki, Denny Wu and Atsushi Nitanda  “Convergence of mean-field Langevin dynamics: time and space discretization, stochastic gradient, and variance reduction”  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2023
  * [Szn91] Alain-Sol Sznitman  “Topics in propagation of chaos”  In _École d’Été de Probabilités de Saint-Flour XIX—1989_ 1464, Lecture Notes in Math.  Springer, Berlin, 1991, pp. 165–251
  * [Tal96] Denis Talay  “Probabilistic numerical methods for partial differential equations: elements of analysis”  In _Probabilistic models for nonlinear partial differential equations (Montecatini Terme, 1995)_ 1627, Lecture Notes in Math.  Springer, Berlin, 1996, pp. 148–196
  * [VW19] Santosh Vempala and Andre Wibisono  “Rapid convergence of the unadjusted Langevin algorithm: Isoperimetry suffices”  In _Advances in Neural Information Processing Systems (NeurIPS)_ 32, 2019
  * [Vil02] Cédric Villani  “A review of mathematical topics in collisional kinetic theory”  In _Handbook of mathematical fluid dynamics, Vol. I_ North-Holland, Amsterdam, 2002, pp. 71–305
  * [Vil03] Cédric Villani  “Topics in optimal transportation” 58, Graduate Studies in Mathematics  American Mathematical Society, Providence, RI, 2003, pp. xvi+370
  * [Wib18] Andre Wibisono  “Sampling as optimization in the space of measures: the Langevin dynamics as a composite optimization problem”  In _Proceedings of the 31st Conference on Learning Theory_ 75, Proceedings of Machine Learning Research  PMLR, 2018, pp. 2093–3027
  * [WSC22] Keru Wu, Scott Schmidler and Yuansi Chen  “Minimax mixing time of the Metropolis-adjusted Langevin algorithm for log-concave sampling”  In _The Journal of Machine Learning Research (JMLR)_ 23.1 JMLRORG, 2022, pp. 12348–12410
  * [YY23] Rentian Yao and Yun Yang  “Mean-field variational inference via Wasserstein gradient flow”  In _arXiv preprint 2207.08074_ , 2023
  * [YKW22] Man-Chung Yue, Daniel Kuhn and Wolfram Wiesemann  “On linear optimization over Wasserstein balls”  In _Math. Program._ 195.1-2, 2022, pp. 1107–1122
  * [Zha+23] Matthew S. Zhang et al.  “Improved discretization analysis for underdamped Langevin Monte Carlo”  In _Conference on Learning Theory (COLT)_ , 2023, pp. 36–71  PMLR

##  Appendix A Control of the Finite-Particle Error

In this section, we prove the results in §3.1 on the finite-particle error. We
will make extensive use of the following transport inequality, which arises as
a consequence of (LSI).

######  Lemma 14 (Talagrand’s Transport Inequality, [OV00]).

If a measure $\pi$ satisfies (LSI) with constant $C_{\mathsf{LSI}}$, then for
all measures $\mu\in\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$,

$\displaystyle\mathcal{W}_{2}^{2}(\mu,\pi)\leq
2C_{\mathsf{LSI}}\,\mathsf{KL}(\mu\mathbin{\|}\pi)\,.$ (TI)

###  A.1 LSI Case

We provide the proof of Theorem 3 under the assumption of (LSI) for the
invariant measures of $\eqref{eq:mean_field_interactive}$ and
$\eqref{eq:finite_particle_interactive}$. This relies on a BBGKY hierarchy
based on the arguments of [LL23].

Recall that $\mu^{1:k}$ is the $k$-particle distribution of the finite-
particle system. Explicitly,

$\displaystyle\log\mu^{1:k}(x^{1:k})$
$\displaystyle=\log\int\exp\Bigl{(}-\frac{2}{\sigma^{2}}\sum_{i=1}^{N}V(x^{i})-\frac{1}{\sigma^{2}\,(N-1)}\sum_{\begin{subarray}{c}i,j=1\\\
i\neq
j\end{subarray}}^{N}W(x^{i}-x^{j})\Bigr{)}\,\mathrm{d}x^{k+1:N}+\text{const.}$

Using exchangeability, we can then compute the gradient of the potential for
this measure as

$\displaystyle-\frac{\sigma^{2}}{2}\,\nabla_{x^{i}}\log\mu^{1:k}(x^{1:k})=\nabla
V(x^{i})+\frac{1}{N-1}\sum_{\begin{subarray}{c}j=1\\\ i\neq
j\end{subarray}}^{k}\nabla
W(x^{i}-x^{j})+\frac{N-k}{N-1}\operatorname{\mathbb{E}}_{\mu^{k+1|1:k}(\cdot\mid
x^{1:k})}\nabla W(x^{i}-\cdot)\,.$

Let $X^{1:k}\sim\mu^{1:k}$ and introduce the notation

$\displaystyle\mathsf{K}_{k}\coloneqq\mathsf{KL}(\mu^{1:k}\mathbin{\|}\pi^{\otimes
k})\,.$

Invoking (LSI) of the mean-field invariant measure (and tensorizing) leads to

$\displaystyle\mathsf{K}_{k}$
$\displaystyle\leq\frac{C_{\mathsf{LSI}}(\pi)}{2}\,\mathsf{FI}(\mu^{1:k}\mathbin{\|}\pi^{\otimes
k})$
$\displaystyle=\frac{2C_{\mathsf{LSI}}(\pi)}{\sigma^{4}}\sum_{i=1}^{k}\operatorname{\mathbb{E}}\Bigl{[}\Bigl{\lVert}\frac{1}{N-1}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{k}\nabla W(X^{i}-X^{j})-\int\nabla
W(X^{i}-\cdot)\,\mathrm{d}\pi$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad{}+\frac{N-k}{N-1}\int\nabla
W(X^{i}-\cdot)\,\mathrm{d}\mu^{k+1|1:k}(\cdot\mid
X^{1:k})\Bigr{\rVert}^{2}\Bigr{]}$
$\displaystyle\leq\frac{4k\,C_{\mathsf{LSI}}(\pi)}{\sigma^{4}\,{(N-1)}^{2}}\,\underset{\mathsf{A}}{\underbrace{\operatorname{\mathbb{E}}\Bigl{[}\Bigl{\lVert}\sum_{\begin{subarray}{c}j=2\end{subarray}}^{k}\Bigl{(}\nabla
W(X^{1}-X^{j})-\int\nabla
W(X^{1}-\cdot)\,\mathrm{d}\pi\Bigr{)}\Bigr{\rVert}^{2}\Bigr{]}}}$
$\displaystyle\qquad{}+\frac{4k\,C_{\mathsf{LSI}}(\pi)\,{(N-k)}^{2}}{\sigma^{4}\,{(N-1)}^{2}}\,\underset{\mathsf{B}}{\underbrace{\operatorname{\mathbb{E}}\Bigl{[}\Bigl{\lVert}\int\nabla
W(X^{1}-\cdot)\,\bigl{(}\mathrm{d}\mu^{k+1|1:k}(\cdot\mid
X^{1:k})-\mathrm{d}\pi\bigr{)}\Bigr{\rVert}^{2}\Bigr{]}}}\,,$

where the last line follows from exchangeability and $\lVert a+b\rVert^{2}\leq
2\,(\lVert a\rVert^{2}+\lVert b\rVert^{2})$ for vectors
$a,b\in\mathbb{R}^{d}$.

####  A.1.1 Bounding the Error Terms

We now handle terms $\mathsf{A},\mathsf{B}$ separately.

$\displaystyle\mathsf{A}$
$\displaystyle=\sum_{j=2}^{k}\operatorname{\mathbb{E}}[\lVert\nabla
W(X^{1}-X^{j})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(X^{1}-\cdot)\rVert^{2}]$
$\displaystyle\qquad+\sum_{\begin{subarray}{c}i,j=2\\\ i\neq
j\end{subarray}}^{k}\operatorname{\mathbb{E}}\bigl{\langle}\nabla
W(X^{1}-X^{i})-\operatorname{\mathbb{E}}_{\pi}\nabla W(X^{1}-\cdot),\,\nabla
W(X^{1}-X^{j})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(X^{1}-\cdot)\bigr{\rangle}$
$\displaystyle\underset{\text{(i)}}{=}(k-1)\operatorname{\mathbb{E}}[\lVert\nabla
W(X^{1}-X^{2})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(X^{1}-\cdot)\rVert^{2}]+$
$\displaystyle\qquad+(k-1)\,(k-2)\operatorname{\mathbb{E}}\bigl{\langle}\nabla
W(X^{1}-X^{2})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(X^{1}-\cdot),W(X^{1}-X^{3})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(X^{1}-\cdot)\bigr{\rangle}$
$\displaystyle\underset{\text{(ii)}}{\leq}(k-1)\,\beta_{W}^{2}\operatorname{\mathbb{E}}[\lVert
X-Y\rVert^{2}]$
$\displaystyle\qquad+(k-1)^{2}\operatorname{\mathbb{E}}\bigl{\langle}\nabla
W(X^{1}-X^{2})-\operatorname{\mathbb{E}}_{\pi}\nabla W(X^{1}-\cdot),\,\nabla
W(X^{1}-X^{3})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(X^{1}-\cdot)\bigr{\rangle}\,,$

where we used the exchangeability of the particles in (i) and the smoothness
of $W$ in (ii). Here, $X\sim\mu^{1}$ and $Y\sim\pi$ are independent.

Let us deal with these two terms separately. For the first term, let
$\bar{Y}\sim\pi$ be _optimally_ coupled with $X$. Then, by independence and
sub-Gaussian concentration (implied by (LSI)),

$\displaystyle\operatorname{\mathbb{E}}[\lVert X-Y\rVert^{2}]$
$\displaystyle\leq 2\operatorname{\mathbb{E}}[\lVert
X-\bar{Y}\rVert^{2}]+2\operatorname{\mathbb{E}}[\lVert
Y-\bar{Y}\rVert^{2}]=2\,\mathcal{W}_{2}^{2}(\mu^{1},\pi)+4\operatorname{\mathbb{E}}[\lVert
Y-\operatorname{\mathbb{E}}Y\rVert^{2}]$ $\displaystyle\leq
4C_{\mathsf{LSI}}(\pi)\operatorname{\mathsf{KL}}(\mu^{1}\mathbin{\|}\pi)+4dC_{\mathsf{LSI}}(\pi)\leq
4C_{\mathsf{LSI}}(\pi)\,(\mathsf{K}_{3}+d)\,,$ (A.1)

where the second inequality follows from (TI), and the last one follows from
the data-processing inequality for the $\operatorname{\mathsf{KL}}$
divergence.

For the second term, the Cauchy–Schwarz inequality leads to

$\displaystyle\operatorname{\mathbb{E}}\bigl{\langle}\nabla
W(X^{1}-X^{2})-\operatorname{\mathbb{E}}_{\pi}\nabla W(X^{1}-\cdot),\,\nabla
W(X^{1}-X^{3})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(X^{1}-\cdot)\bigr{\rangle}$
$\displaystyle\qquad=\operatorname{\mathbb{E}}\bigl{\langle}\nabla
W(X^{1}-X^{2})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(X^{1}-\cdot),\,\operatorname{\mathbb{E}}_{\mu^{3\mid 1:2}(\cdot\mid
X^{1:2})}\nabla W(X^{1}-\cdot)-\operatorname{\mathbb{E}}_{\pi}\nabla
W(X^{1}-\cdot)\bigr{\rangle}$
$\displaystyle\qquad\leq\beta_{W}^{2}\sqrt{\operatorname{\mathbb{E}}[\lVert
X-Y\rVert^{2}]}\sqrt{\operatorname{\mathbb{E}}\mathcal{W}_{2}^{2}\bigl{(}\mu^{3\mid
1:2}(\cdot\mid X^{1:2}),\,\pi\bigr{)}}$
$\displaystyle\qquad\underset{\text{(i)}}{\leq}\beta_{W}^{2}\sqrt{4C_{\mathsf{LSI}}(\pi)\,(\mathsf{K}_{3}+d)}\sqrt{2C_{\mathsf{LSI}}(\pi)\operatorname{\mathbb{E}}\operatorname{\mathsf{KL}}\bigl{(}\mu^{3\mid
1:2}(\cdot\mid X^{1:2})\bigm{\|}\pi\bigr{)}}$
$\displaystyle\qquad\underset{\text{(ii)}}{\leq}3\beta_{W}^{2}C_{\mathsf{LSI}}(\pi)\sqrt{\mathsf{K}_{3}+d}\sqrt{\mathsf{K}_{3}}$
(A.2) $\displaystyle\qquad\leq
3\beta_{W}^{2}C_{\mathsf{LSI}}(\pi)\,(\mathsf{K}_{3}+d)\,,$

where in (i) we applied the bound (A.1) as well as (TI), and in (ii) we used
the chain rule for the $\mathsf{KL}$ divergence.

We return to the analysis of the term $\mathsf{B}$. In a similar way, we
obtain

$\displaystyle\mathsf{B}$
$\displaystyle=\operatorname{\mathbb{E}}\Bigl{[}\Bigl{\lVert}\int\nabla
W(X^{1}-\cdot)\,\bigl{(}\mathrm{d}\mu^{k+1|1:k}(\cdot\mid
X^{1:k})-\mathrm{d}\pi\bigr{)}\Bigr{\rVert}^{2}\Bigr{]}\leq\beta_{W}^{2}\operatorname{\mathbb{E}}\mathcal{W}_{2}^{2}\bigl{(}\mu^{k+1|1:k}(\cdot\mid
X^{1:k}),\,\pi\bigr{)}$ $\displaystyle\leq
2\beta_{W}^{2}C_{\mathsf{LSI}}(\pi)\,(\mathsf{K}_{k+1}-\mathsf{K}_{k})\,.$

####  A.1.2 Induction

Putting our bounds on $\mathsf{A}$ and $\mathsf{B}$ together, we obtain for
$N\geq 30$,

$\mathsf{K}_{k}\leq\frac{30k^{3}\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)}{\sigma^{4}N^{2}}\,(\mathsf{K}_{3}+d)+\frac{8k\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)}{\sigma^{4}}\,(\mathsf{K_{k+1}-\mathsf{K}_{k})}\,.$
(A.3)

In particular, the case of $k=N$ involves our bounds only on $\mathsf{A}$,
leading to

$\mathsf{K}_{N}\leq\frac{30N\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)}{\sigma^{4}}\,(\mathsf{K}_{3}+d)\,.$

By grouping together the $\mathsf{K}_{k}$ terms in (A.3),

$\displaystyle\mathsf{K}_{k}\leq\underset{=:\mathcal{C}_{k}}{\underbrace{\frac{8k\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)/\sigma^{4}}{1+8k\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)/\sigma^{4}}}}\,\Bigl{(}\mathsf{K}_{k+1}+\bigl{(}\frac{2k}{N}\bigr{)}^{2}\,(\mathsf{K}_{3}+d)\Bigr{)}\,.$
(A.4)

Iterating this inequality down to $k=3$, for
$\rho\coloneqq\nicefrac{{\sigma^{4}}}{{8\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)}}$,

$\displaystyle\mathsf{K}_{3}$
$\displaystyle\leq\Bigl{(}\prod_{k=3}^{N-1}\mathcal{C}_{k}\Bigr{)}\,\frac{30N\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)}{\sigma^{4}}\,(\mathsf{K}_{3}+d)+\sum_{k=3}^{N-1}\Bigl{(}\prod_{\ell=3}^{k}\mathcal{C}_{\ell}\Bigr{)}\,\bigl{(}\frac{2k}{N}\bigr{)}^{2}\,(\mathsf{K}_{3}+d)$
$\displaystyle\leq\underbrace{\biggl{[}\Bigl{(}\prod_{k=3}^{N-1}\mathcal{C}_{k}\Bigr{)}\,\frac{4N}{\rho}+\sum_{k=3}^{N-1}\Bigl{(}\prod_{\ell=3}^{k}\mathcal{C}_{\ell}\Bigr{)}\,\bigl{(}\frac{2k}{N}\bigr{)}^{2}\biggr{]}}_{\eqqcolon
c_{N}}\,(\mathsf{K}_{3}+d)\,.$

Now we show $c_{N}<1/2$, which implies $\mathsf{K}_{3}\leq 2c_{N}d$. We
require the following lemma.

######  Lemma 15.

For $3\leq i\leq k\leq N$,

$\prod_{\ell=i}^{k}\mathcal{C}_{\ell}\leq\Bigl{(}\frac{i+\rho}{k+1+\rho}\Bigr{)}^{\rho}\,.$

Proof. For $\mathcal{C}_{\ell}=\frac{\ell\rho^{-1}}{1+\ell\rho^{-1}}$, we have

$C\coloneqq\log\prod_{\ell=i}^{k}\mathcal{C}_{\ell}=\sum_{\ell=i}^{k}\log\Bigl{(}1-\frac{1}{1+\ell\rho^{-1}}\Bigr{)}\leq-\sum_{\ell=i}^{k}\frac{1}{1+\ell\rho^{-1}}\,.$

As the summand is decreasing in $\ell$, it follows that

$C\leq-\sum_{\ell=i}^{k}\int_{\ell}^{\ell+1}\frac{1}{1+x\rho^{-1}}\,\mathrm{d}x=-\int_{i}^{k+1}\frac{1}{1+x\rho^{-1}}\,\mathrm{d}x=-\rho\log\frac{k+1+\rho}{i+\rho}\,.$

Therefore,

$\displaystyle\prod_{\ell=i}^{k}\mathcal{C}_{\ell}=\exp
C\leq\Bigl{(}\frac{k+1+\rho}{i+\rho}\Bigr{)}^{-\rho}\,,$

which proves the lemma. ∎

Using Lemma 15, we obtain

$c_{N}\leq
4\,{(3+\rho)}^{\rho}\,\Bigl{(}\frac{N^{1-\rho}}{\rho}+\frac{1}{N^{2}}\sum_{k=3}^{N-1}k^{2-\rho}\Bigr{)}\,.$

Under Assumption 3, i.e., $\rho\geq 3$, we may assume $\rho=3$ since we can
always take a worse bound on the constants $\beta_{W}$ so that $\rho=3$. As
seen shortly, the rate does not improve even if $\rho>3$.555Alternatively, one
can show the bound in Lemma 15 decreases in $\rho$, so we can just substitute
$\rho=3$ therein. For $\rho=3$ and $N\geq 100$, we therefore obtain

$c_{N}\leq
864\,\Bigl{(}\frac{1}{3N^{2}}+\frac{1}{N^{2}}\sum_{k=3}^{N-1}\frac{1}{k}\Bigr{)}\leq\frac{1}{2}\,,$

and thus

$\displaystyle\mathsf{K}_{3}\lesssim\frac{d\log N}{N^{2}}\,.$ (A.5)

####  A.1.3 Bootstrapping

Substituting the bound (A.5) for $\mathsf{K}_{3}$ into the recursive
inequality (A.4), we end up with a suboptimal rate of
$\widetilde{\mathcal{O}}(k^{3}/N^{2})$ for $\mathsf{K}_{k}$. To improve the
bound, we substitute our established bound (A.5) into (A.2), which results in
an improved recursive inequality. Indeed,

$\displaystyle\mathsf{A}$ $\displaystyle\lesssim
k\beta_{W}^{2}C^{2}_{\mathsf{LSI}}(\pi)\,({\operatorname{\mathsf{KL}}_{3}}+d)+k^{2}\beta_{W}^{2}C_{\mathsf{LSI}}(\pi)\sqrt{{\operatorname{\mathsf{KL}}_{3}}+d}\sqrt{\operatorname{\mathsf{KL}}_{3}}\lesssim
dk\beta_{W}^{2}C_{\mathsf{LSI}}(\pi)\sqrt{\log N}$

and therefore

$\displaystyle\mathsf{K}_{k}$
$\displaystyle\leq\widetilde{\mathcal{O}}\Bigl{(}\frac{dk^{2}\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)}{\sigma^{4}N^{2}}\Bigr{)}+\frac{8k\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)}{\sigma^{4}}\,(\mathsf{K}_{k+1}-\mathsf{K}_{k})\,.$

For $k=N$ this yields

$\displaystyle\mathsf{K}_{N}$
$\displaystyle\leq\widetilde{\mathcal{O}}\Bigl{(}\frac{d\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)}{\sigma^{4}}\Bigr{)}\,.$

Regrouping $\mathsf{K}_{k}$ as before, we obtain

$\mathsf{K}_{k}\leq\mathcal{C}_{k}\,\Bigl{(}\mathsf{K}_{k+1}+\widetilde{\mathcal{O}}\bigl{(}\frac{dk}{N^{2}}\bigr{)}\Bigr{)}\,.$

Iterating this down to $k=N$,

$\displaystyle\mathsf{K}_{k}$
$\displaystyle\leq\Bigl{(}\prod_{\ell=k}^{N-1}\mathcal{C}_{\ell}\Bigr{)}\,\mathsf{K}_{N}+\sum_{\ell=k}^{N-1}\Bigl{(}\prod_{j=k}^{\ell}\mathcal{C}_{j}\Bigr{)}\,\widetilde{\mathcal{O}}\bigl{(}\frac{d\ell}{N^{2}}\bigr{)}$
$\displaystyle\underset{\text{(i)}}{\leq}\widetilde{\mathcal{O}}\Bigl{(}\frac{k^{3}}{N^{3}}\,\frac{d\beta_{W}^{2}C_{\mathsf{LSI}}^{2}(\pi)}{\sigma^{4}}+\sum_{\ell=k}^{N-1}\frac{k^{3}}{\ell^{3}}\,\frac{d\ell}{N^{2}}\Bigr{)}\underset{\text{(ii)}}{\leq}\widetilde{\mathcal{O}}\bigl{(}\frac{dk^{2}}{N^{2}}\bigr{)}\,,$

where in (i) we used Lemma 15 with $\rho=3$, and (ii) follows from $\rho\geq
3$ and $\sum_{\ell\geq k}\ell^{-2}\leq k^{-1}$. Therefore, for some fixed $k$
it suffices to take $N=100\vee\widetilde{\Omega}(k\sqrt{d}/\varepsilon)$ to
achieve $\varepsilon^{2}$-bias in $\mathsf{KL}$, completing the proof of
Theorem 3.

###  A.2 Strongly Convex Case

The following propagation of chaos argument for the strongly log-concave case
is based on [Szn91]. Let $(X^{1:N}_{t})_{t\geq 0}$ denote the stochastic
process following the finite-particle stochastic differential equation
($\mathsf{pMV}_{N}$). Let the corresponding semigroup be denoted
${(\mathcal{T}_{t})}_{t\geq 0}$, defined as follows. For any test function
$f:\mathbb{R}^{d\times N}\to\mathbb{R}$,

$\displaystyle\mathcal{T}_{t}f(x^{1:N})=\operatorname{\mathbb{E}}[f(X^{1:N}_{t})\mid
X_{0}^{1:N}=x^{1:N}]\,.$

Then, the following simple lemma proves Wasserstein contraction for the
finite-particle system.

######  Lemma 16.

Under Assumption 4 and for
$N\geq\frac{\alpha_{V}-{\alpha_{W}^{-}}}{\alpha_{V}+{(\alpha_{W})}_{-}}$,
${(\mathcal{T}_{t})}_{t\geq 0}$ is a contraction in the $2$-Wasserstein
distance with exponential rate at least $\alpha/2$, where
$\alpha\coloneqq\alpha_{V}+{\alpha_{W}^{-}}$. In other words, for any measures
$\mu_{0}^{1:N}$, $\nu_{0}^{1:N}$ in $\mathcal{P}_{2}(\mathbb{R}^{d\times N})$,

$\displaystyle\mathcal{W}_{2}(\mu_{0}^{1:N}\mathcal{T}_{t},\,\nu_{0}^{1:N}\mathcal{T}_{t})$
$\displaystyle\leq\exp(-\alpha
t/2)\,\mathcal{W}_{2}(\mu_{0}^{1:N},\nu_{0}^{1:N})\,.$

Proof. Note that ${(\mathcal{T}_{t})}_{t\geq 0}$ corresponds to the time-
scaled (by factor $\sigma^{2}/2$) Langevin diffusion with stationary
distribution $\mu^{1:N}$, which is
$\frac{2}{\sigma^{2}}\,(\alpha_{V}+\frac{N}{N-1}\,{\alpha_{W}^{-}})$-strongly
log-concave by Lemma 7. The condition on $N$ ensures that this is at least
$\alpha/\sigma^{2}$. Consequently, it is well-known (e.g., via synchronous
coupling) that the diffusion is a contraction in the Wasserstein distance with
rate at least $\alpha/2$. ∎

We next bound the error incurred in one step from applying the finite-particle
semigroup to $\pi^{\otimes N}$.

######  Lemma 17.

Under Assumptions 1 and 4, for any $\lambda>0$, $\mathcal{T}_{h}$ induces the
following error in Wasserstein distance:

$\displaystyle\mathcal{W}_{2}^{2}(\pi^{\otimes
N}\mathcal{T}_{h},\,\pi^{\otimes N})$
$\displaystyle\leq\frac{(1+\lambda^{-1})\,\beta_{W}^{2}\sigma^{2}dh^{2}}{\alpha}\exp\Bigl{(}\frac{(1+\lambda)\,\beta_{W}^{2}h^{2}}{2}\Bigr{)}\,.$

Proof. We resort to a coupling argument, noting that $\pi$ is stationary under
($\mathsf{pMV}$). Starting with $\pi^{\otimes N}$, we evolve
${(X^{1:N}_{t})}_{t\geq 0}$ and ${(Y^{1:N}_{t})}_{t\geq 0}$ according to
($\mathsf{pMV}_{N}$) and ($\mathsf{pMV}$) respectively, i.e.,
$X^{1:N}_{t}\sim\pi^{\otimes N}\mathcal{T}_{t}$ and
$Y^{1:N}_{t}\sim\pi^{\otimes N}$. This argument is adapted from the original
propagation of chaos proof by [Szn91].

We can compute the evolution under a synchronous coupling as:

$\displaystyle\mathrm{d}(X^{i}_{t}-Y^{i}_{t})$ $\displaystyle=-\bigl{(}\nabla
V(X^{i}_{t})-\nabla
V(Y^{i}_{t})\bigr{)}\,\mathrm{d}t-\frac{1}{N-1}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\bigl{(}\nabla
W(X^{i}_{t}-X^{j}_{t})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(Y^{i}_{t}-\cdot)\bigr{)}\,\mathrm{d}t$ $\displaystyle=-\bigl{(}\nabla
V(X^{i}_{t})-\nabla
V(Y^{i}_{t})\bigr{)}\,\mathrm{d}t-\frac{1}{N-1}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\bigl{(}\nabla W(X^{i}_{t}-X^{j}_{t})-\nabla
W(Y^{i}_{t}-X^{j}_{t})\bigr{)}\,\mathrm{d}t$
$\displaystyle\qquad{}-\frac{1}{N-1}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}\bigl{(}\nabla W(Y^{i}_{t}-X^{j}_{t})-\nabla
W(Y^{i}_{t}-Y^{j}_{t})\bigr{)}\,\mathrm{d}t$
$\displaystyle\qquad{}-\frac{1}{N-1}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}\bigl{(}\nabla
W(Y^{i}_{t}-Y^{j}_{t})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(Y^{i}_{t}-\cdot)\bigr{)}\,\mathrm{d}t\,.$

Now let us denote by $\overline{\nabla W}(x,y)\coloneqq\nabla
W(x-y)-\operatorname{\mathbb{E}}_{\pi}\nabla W(x-\cdot)$ the centered gradient
(with respect to $\pi$). By Itô’s formula and Assumption 4,

$\displaystyle\mathrm{d}\lVert X_{t}^{i}-Y_{t}^{i}\rVert^{2}$
$\displaystyle=2\,\langle
X_{t}^{i}-Y_{t}^{i},\mathrm{d}(X_{t}^{i}-Y_{t}^{i})\rangle$
$\displaystyle\leq-2\,(\alpha_{V}+\alpha_{W})\,\lVert
X_{t}^{i}-Y_{t}^{i}\rVert^{2}\,\mathrm{d}t$
$\displaystyle\qquad{}-\frac{2}{N-1}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}\langle X_{t}^{i}-Y_{t}^{i},\nabla
W(Y^{i}_{t}-X^{j}_{t})-\nabla W(Y^{i}_{t}-Y^{j}_{t})\rangle\,\mathrm{d}t$
$\displaystyle\qquad{}-\frac{2}{N-1}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}\langle X_{t}^{i}-Y_{t}^{i},\nabla
W(Y^{i}_{t}-Y^{j}_{t})-\operatorname{\mathbb{E}}_{\pi}\nabla
W(Y^{i}_{t}-\cdot)\rangle\,\mathrm{d}t$
$\displaystyle\leq\frac{2\beta_{W}\,\lVert
X_{t}^{i}-Y_{t}^{i}\rVert}{N-1}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}{\lVert
X_{t}^{j}-Y_{t}^{j}\rVert}\,\mathrm{d}t+\frac{2\,\lVert
X_{t}^{i}-Y_{t}^{i}\rVert}{N-1}\,\Bigl{\lVert}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\overline{\nabla
W}(Y_{t}^{i},Y_{t}^{j})\Bigr{\rVert}\,\mathrm{d}t$

or

$\displaystyle\mathrm{d}\lVert X_{t}^{i}-Y_{t}^{i}\rVert$
$\displaystyle\leq\frac{\beta_{W}}{N-1}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}{\lVert
X_{t}^{j}-Y_{t}^{j}\rVert}\,\mathrm{d}t+\frac{1}{N-1}\,\Bigl{\lVert}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\overline{\nabla
W}(Y_{t}^{i},Y_{t}^{j})\Bigr{\rVert}\,\mathrm{d}t\,.$

Integrating and squaring,

$\displaystyle\lVert X_{t}^{i}-Y_{t}^{i}\rVert^{2}$
$\displaystyle\leq\Bigl{\lvert}\int_{0}^{t}\Bigl{(}\frac{\beta_{W}}{N-1}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}{\lVert
X_{s}^{j}-Y_{s}^{j}\rVert}+\frac{1}{N-1}\,\Bigl{\lVert}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\overline{\nabla
W}(Y_{s}^{i},Y_{s}^{j})\Bigr{\rVert}\Bigr{)}\,\mathrm{d}s\Bigr{\rvert}^{2}$
$\displaystyle\leq\frac{(1+\lambda)\,\beta_{W}^{2}t}{N-1}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\int_{0}^{t}\lVert
X_{s}^{j}-Y_{s}^{j}\rVert^{2}\,\mathrm{d}s+\frac{(1+\lambda^{-1})\,t}{{(N-1)}^{2}}\int_{0}^{t}\Bigl{\lVert}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}\overline{\nabla
W}(Y_{s}^{i},Y_{s}^{j})\Bigr{\rVert}^{2}\,\mathrm{d}s\,,$

where the last line follows from Young’s inequality.

Next, we take expectations. Note that $\overline{\nabla W}(\cdot,\cdot)$ is
centered in its second variable, so for any $j\neq k$,

$\displaystyle\operatorname{\mathbb{E}}\langle\overline{\nabla
W}(Y^{i}_{t},Y^{j}_{t}),\overline{\nabla W}(Y^{i}_{t},Y^{k}_{t})\rangle=0\,.$

Otherwise, we can bound the terms via

$\displaystyle\operatorname{\mathbb{E}}[\lVert\overline{\nabla
W}(Y^{i}_{t},Y^{j}_{t})\rVert^{2}]\leq\beta_{W}^{2}\operatorname{\mathbb{E}}_{\begin{subarray}{c}Y^{j}_{t}\sim\pi\\\
Z\sim\pi\end{subarray}}[\lVert
Y^{j}_{t}-Z\rVert^{2}]\leq\frac{\beta_{W}^{2}\sigma^{2}d}{\alpha}\,.$

Here, $Z$ is an independent draw from $\pi$ and so cannot be reduced via
coupling. The second inequality follows from a standard bound on the centered
second moment of a strongly log-concave measure, using the fact that $\pi$ is
$2\alpha/\sigma^{2}$-strongly log-concave [[, c.f.]]DalKarRio22NonStrongly.

Therefore, taking expectations and summing over the particles,

$\displaystyle\operatorname{\mathbb{E}}[\lVert
X_{t}^{1:N}-Y_{t}^{1:N}\rVert^{2}]$
$\displaystyle\leq(1+\lambda)\,\beta_{W}^{2}t\int_{0}^{t}\lVert
X_{s}^{1:N}-Y_{s}^{1:N}\rVert^{2}\,\mathrm{d}s+\frac{(1+\lambda^{-1})\,\beta_{W}^{2}\sigma^{2}dt^{2}}{\alpha}\,.$

By Grönwall’s inequality below,

$\displaystyle\operatorname{\mathbb{E}}[\lVert
X_{h}^{1:N}-Y_{h}^{1:N}\rVert^{2}]$
$\displaystyle\leq\frac{(1+\lambda^{-1})\,\beta_{W}^{2}\sigma^{2}dh^{2}}{\alpha}\exp\Bigl{(}\frac{(1+\lambda)\,\beta_{W}^{2}h^{2}}{2}\Bigr{)}\,.$

This concludes the proof. ∎

######  Lemma 18 (Grönwall’s Inequality).

For $T>0$, let $f:[0,T]\to\mathbb{R}_{\geq 0}$ be bounded. Suppose that the
following holds pointwise for some functions $a,b:[0,T]\to\mathbb{R}$, where
$a$ is increasing:

$f(t)\leq a(t)+\int_{0}^{t}b(s)f(s)\,\mathrm{d}s\,.$

Then,

$f(t)\leq a(t)\exp\Bigl{(}\int_{0}^{t}b(s)\,\mathrm{d}s\Bigr{)}\,.$

Composing Lemmas 16 and 17, we now prove our propagation of chaos results.

Proof of Theorem 4 Indeed, we have

$\displaystyle\mathcal{W}_{2}(\mu^{1:N},\pi^{\otimes
N})=\mathcal{W}_{2}(\mu^{1:N}\mathcal{T}_{h},\,\pi^{\otimes
N})\leq\mathcal{W}_{2}(\mu^{1:N}\mathcal{T}_{h},\,\pi^{\otimes
N}\mathcal{T}_{h})+\mathcal{W}_{2}(\pi^{\otimes
N}\mathcal{T}_{h},\,\pi^{\otimes N})$ $\displaystyle\qquad\leq\exp(-\alpha
h/2)\,\mathcal{W}_{2}(\mu^{1:N},\pi^{\otimes
N})+\sqrt{\frac{(1+\lambda^{-1})\,\beta_{W}^{2}\sigma^{2}dh^{2}}{\alpha}}\exp\Bigl{(}\frac{(1+\lambda)\,\beta_{W}^{2}h^{2}}{4}\Bigr{)}\,.$

Rearranging,

$\displaystyle\mathcal{W}_{2}(\mu^{1:N},\pi^{\otimes N})$
$\displaystyle\leq\frac{1}{1-\exp(-\alpha
h/2)}\sqrt{\frac{(1+\lambda^{-1})\,\beta_{W}^{2}\sigma^{2}dh^{2}}{\alpha}}\exp\Bigl{(}\frac{(1+\lambda)\,\beta_{W}^{2}h^{2}}{4}\Bigr{)}\,.$

Let $h\searrow 0$ first and then $\lambda\nearrow\infty$ to obtain

$\displaystyle\mathcal{W}_{2}^{2}(\mu^{1:N},\pi^{\otimes N})$
$\displaystyle\leq\frac{4\beta_{W}^{2}\sigma^{2}d}{\alpha^{3}}\,.$

Finally, when $k<N$, we use exchangeability (see Lemma 23 below) to conclude
the proof of (3.2).

For (3.3), by the Bakry–Émery condition we have
$C_{\mathsf{LSI}}(\pi)\leq\nicefrac{{\sigma^{2}}}{{2\alpha}}$, and
tensorization [[, c.f.]Proposition 5.2.7]bakry2014analysis leads to
$C_{\mathsf{LSI}}(\pi^{\otimes N})\leq\nicefrac{{\sigma^{2}}}{{2\alpha}}$.
Thus, (TI) leads to

$\mathsf{KL}(\mu^{1:N}\mathbin{\|}\pi^{\otimes
N})\leq\frac{\sigma^{2}}{4\alpha}\,\mathsf{FI}(\mu^{1:N}\mathbin{\|}\pi^{\otimes
N})\,.$

However, one notes that the density of $\mu^{1:N}$ is log-smooth with
parameter $\frac{2}{\sigma^{2}}\,(\beta_{V}+\frac{N}{N-1}\,\beta_{W})$ (Lemma
7). Likewise, $\pi^{\otimes N}$ is log-smooth with parameter
$\frac{2}{\sigma^{2}}\,(\beta_{V}+\beta_{W})$. Now consider a functional
$\mathcal{F}$ on the space of probability measures on $\mathcal{P}_{2,\rm
ac}(\mathbb{R}^{d\times N})$ given by
$\mathcal{F}:\nu\mapsto\operatorname{\mathbb{E}}_{\nu}[\lVert\nabla\log\frac{\mu^{1:N}}{\pi^{\otimes
N}}\rVert^{2}]$. Note that $\log(\mu^{1:N}/\pi^{\otimes N})$ is smooth with
parameter at most
$\frac{4}{\sigma^{2}}\,(\beta_{V}+\frac{N}{N-1}\,\beta_{W})\leq\frac{8}{\sigma^{2}}\,(\beta_{V}+\beta_{W})$,
for $N\geq 2$.

Next, note that for $Y^{1:N}\sim\pi^{\otimes N}$,

$\displaystyle\mathcal{F}(\pi^{\otimes N})$
$\displaystyle=\operatorname{\mathbb{E}}_{\pi^{\otimes
N}}[\lVert\nabla\log\mu^{1:N}-\nabla\log\pi^{\otimes N}\rVert^{2}]$
$\displaystyle=\frac{4N}{\sigma^{4}\,{(N-1)}^{2}}\operatorname{\mathbb{E}}\Bigl{[}\Bigl{\lVert}\sum_{j=2}^{N}\Bigl{(}\nabla
W(Y^{1}-Y^{j})-\int\nabla
W(Y^{1}-\cdot)\,\mathrm{d}\pi\Bigr{)}\Bigr{\rVert}^{2}\Bigr{]}$
$\displaystyle=\frac{4N}{\sigma^{4}\,{(N-1)}^{2}}\operatorname{\mathbb{E}}\Bigl{[}\Bigl{\lVert}\sum_{j=2}^{N}\overline{\nabla
W}(Y^{1},Y^{j})\Bigr{\rVert}^{2}\Bigr{]}\,,$

by using exchangeability and the definition of $\overline{\nabla W}$.

Subsequently, one derives the following inequality using the Wasserstein
distance bound:

$\displaystyle\mathcal{F}(\mu^{1:N})$
$\displaystyle\leq\frac{128}{\sigma^{4}}\,{(\beta_{V}+\beta_{W})}^{2}\,\mathcal{W}_{2}^{2}(\mu^{1:N},\pi^{\otimes
N})+2\mathcal{F}(\pi^{\otimes N})$
$\displaystyle\leq\frac{512\beta_{W}^{2}d}{\alpha^{3}\sigma^{2}}\,{(\beta_{V}+\beta_{W})}^{2}+\frac{8N}{\sigma^{4}\,{(N-1)}^{2}}\operatorname{\mathbb{E}}\Bigl{[}\Bigl{\lVert}\sum_{j=2}^{N}\overline{\nabla
W}(Y^{1},Y^{j})\Bigr{\rVert}^{2}\Bigr{]}$
$\displaystyle\leq\frac{512\beta_{W}^{2}d}{\alpha^{3}\sigma^{2}}\,{(\beta_{V}+\beta_{W})}^{2}+\frac{16\beta_{W}^{2}N}{\sigma^{4}\,(N-1)}\operatorname{\mathbb{E}}[\lVert
Y^{1}-\operatorname{\mathbb{E}}Y^{1}\rVert^{2}]$
$\displaystyle\leq\frac{512\beta_{W}^{2}d}{\alpha^{3}\sigma^{2}}\,{(\beta_{V}+\beta_{W})}^{2}+\frac{16\beta_{W}^{2}d}{\alpha\sigma^{2}}\,,$

by using (3.2) and the fact that $\overline{\nabla W}(\cdot,\cdot)$ is a
centered random variable in its second argument. This concludes the proof for
$k=N$, and as in the $\mathcal{W}_{2}^{2}$ bound, Lemma 24 will conclude the
proof for $k<N$. ∎

###  A.3 General Functional Case

For any measure $\mu$, define its entropy as
$\operatorname{\mathsf{ent}}(\mu)=\int\log\mu\,\mathrm{d}\mu$. We now provide
a self-contained propagation of chaos argument in the general McKean–Vlasov
setting, following [CRW22]. We begin with the following entropy toast
inequality, i.e., half of the entropy sandwich inequality from [CRW22].

######  Lemma 19 (Entropy Toast Inequality).

Define the empirical total energy for an $N$-finite particle system as
follows. Given a measure
$\nu^{1:N}\in\mathcal{P}_{2,\rm{ac}}(\mathbb{R}^{d\times N})$,

$\displaystyle\mathcal{E}^{N}(\nu^{1:N})=N\int\mathcal{F}(\rho_{x^{1:N}})\,\nu^{1:N}(\mathrm{d}x^{1:N})+\frac{\sigma^{2}}{2}\operatorname{\mathsf{ent}}(\mu)\,.$

Under Assumption 5, it holds for all measures $\nu^{1:N}\in\mathcal{P}_{2,\rm
ac}(\mathbb{R}^{d\times N})$

$\frac{\sigma^{2}}{2}\operatorname{\mathsf{KL}}(\nu^{1:N}\mathbin{\|}\pi^{\otimes
N})\leq\mathcal{E}^{N}(\nu^{1:N})-N\mathcal{E}(\pi)\,,$

where $\mathcal{E}$ is the total energy ($\mathsf{gE}$) and $\pi$ is the
stationary measure (2.2).

Proof. By Assumption 5, we have

$\displaystyle\mathcal{E}^{N}(\nu^{1:N})-N\mathcal{E}(\pi)=N\operatorname{\mathbb{E}}_{x^{1:N}\sim\nu^{1:N}}[\mathcal{F}(\rho_{x^{1:N}})-\mathcal{F}(\pi)]+\frac{\sigma^{2}}{2}\,\bigl{(}\operatorname{\mathsf{ent}}(\nu^{1:N})-N\operatorname{\mathsf{ent}}(\pi)\bigr{)}$
$\displaystyle\qquad\geq\operatorname{\mathbb{E}}_{x^{1:N}\sim\nu^{1:N}}\Bigl{[}N\int\delta\mathcal{F}(\pi,z)\,(\rho_{x^{1:N}}(\mathrm{d}z)-\pi(\mathrm{d}z))\Bigr{]}+\frac{\sigma^{2}}{2}\,\bigl{(}\operatorname{\mathsf{ent}}(\nu^{1:N})-N\operatorname{\mathsf{ent}}(\pi)\bigr{)}$
$\displaystyle\qquad\underset{\text{(i)}}{=}-\frac{\sigma^{2}}{2}\operatorname{\mathbb{E}}_{x^{1:N}\sim\nu^{1:N}}\Bigl{[}N\int\log\pi(z)\,(\rho_{x^{1:N}}(\mathrm{d}z)-\pi(\mathrm{d}z))\Bigr{]}+\frac{\sigma^{2}}{2}\,\bigl{(}\operatorname{\mathsf{ent}}(\nu^{1:N})-N\operatorname{\mathsf{ent}}(\pi)\bigr{)}$
$\displaystyle\qquad=-\frac{\sigma^{2}}{2}\operatorname{\mathbb{E}}_{x^{1:N}\sim\nu^{1:N}}\Bigl{[}N\int\log\pi(z)\,\rho_{x^{1:N}}(\mathrm{d}z)\Bigr{]}+\frac{\sigma^{2}}{2}\operatorname{\mathsf{ent}}(\nu^{1:N})$
$\displaystyle\qquad=-\frac{\sigma^{2}}{2}\int\sum_{i=1}^{N}\log\pi(x^{i})\,\nu^{1:N}(\mathrm{d}x^{1:N})+\frac{\sigma^{2}}{2}\operatorname{\mathsf{ent}}(\nu^{1:N})\,,$

where (i) follows from that $\pi$ is the minimizer of $\mathcal{E}$. The last
line is just
$\frac{\sigma^{2}}{2}\operatorname{\mathsf{KL}}(\nu^{1:N}\mathbin{\|}\pi^{\otimes
N})$, so we are done. ∎

Proof of Theorem 5 We bound $\mathcal{E}^{N}(\mu^{1:N})-N\mathcal{E}^{N}(\pi)$
via the following argument. First, define the finite-particle mean-field
functional as
$\mathcal{F}^{N}(\nu^{1:N})=N\int\mathcal{F}(\rho_{x^{1:N}})\,\nu^{1:N}(\mathrm{d}x^{1:N})$.
In the sequel, we also use the following notation for conditional measures: if
$x^{-i}\coloneqq(x^{1:i-1},x^{i+1:N})\in\mathbb{R}^{d\times(N-1)}$,

$\mu^{1:N}(x^{1:N})=\mu^{i|-i}(x^{i}\mid x^{-i})\times\mu^{-i}(x^{-i})\,.$

We know that

$\displaystyle\mathcal{E}^{N}(\mu^{1:N})-N\mathcal{E}(\pi)=\mathcal{F}^{N}(\mu^{1:N})-N\mathcal{F}(\pi)+\frac{\sigma^{2}}{2}\operatorname{\mathsf{ent}}(\mu^{1:N})-\frac{N\sigma^{2}}{2}\operatorname{\mathsf{ent}}(\pi)\,.$

Furthermore, by Assumption 5,

$\displaystyle\mathcal{F}^{N}(\mu^{1:N})-N\mathcal{F}(\pi)\leq
N\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\int\delta\mathcal{F}(\rho_{x^{1:N}},z)\,(\rho_{x^{1:N}}(\mathrm{d}z)-\pi(\mathrm{d}z))\,.$

Using the subadditivity of entropy, we can therefore write

$\displaystyle\mathcal{E}^{N}(\mu^{1:N})-N\mathcal{E}(\pi)\leq\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\Bigl{[}$
$\displaystyle\delta\mathcal{F}(\rho_{x^{1:N}},x^{i})-\int\delta\mathcal{F}(\rho_{x^{1:N}},\cdot)\,\mathrm{d}\pi$
$\displaystyle\qquad{}+\frac{\sigma^{2}}{2}\,\bigl{(}\operatorname{\mathsf{ent}}(\mu^{i\mid-i}(\cdot\mid
x^{-i}))-\operatorname{\mathsf{ent}}(\pi)\bigr{)}\Bigr{]}\,.$

To decouple the terms, we now replace each
$\delta\mathcal{F}(\rho_{x^{1:N}},\cdot)$ term with
$\delta\mathcal{F}(\rho_{x^{-i}},\cdot)$:

$\displaystyle\mathcal{E}^{N}(\mu^{1:N})-N\mathcal{E}(\pi)$
$\displaystyle\qquad\leq\underbrace{\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\Bigl{[}\delta\mathcal{F}(\rho_{x^{-i}},x^{i})-\int\delta\mathcal{F}(\rho_{x^{-i}},\cdot)\,\mathrm{d}\pi+\frac{\sigma^{2}}{2}\,\bigl{(}\operatorname{\mathsf{ent}}(\mu^{i\mid-i}(\cdot\mid
x^{-i}))-\operatorname{\mathsf{ent}}(\pi)\bigr{)}\Bigr{]}}_{\mathsf{A}}$
$\displaystyle\qquad\quad{}+{\underbrace{\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\Bigl{[}\delta\mathcal{F}(\rho_{x^{1:N}},x^{i})-\delta\mathcal{F}(\rho_{x^{-i}},x^{i})-\int\bigl{(}\delta\mathcal{F}(\rho_{x^{1:N}},\cdot)-\delta\mathcal{F}(\rho_{x^{-i}},\cdot)\bigr{)}\,\mathrm{d}\pi\Bigr{]}}_{\mathsf{B}}}\,.$

We consider the two terms in turn, beginning with the first.

Note that by Fubini’s theorem,

$\displaystyle\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\delta\mathcal{F}(\rho_{x^{-i}},x^{i})$
$\displaystyle=\operatorname{\mathbb{E}}_{x^{-i}\sim\mu^{-i}}\int\delta\mathcal{F}(\rho_{x^{-i}},\cdot)\,\mathrm{d}\mu^{i\mid-i}(\cdot\mid
x^{-i})\,.$

In order to relate the first term $\mathsf{A}$ to the
$\operatorname{\mathsf{KL}}$ divergence, for each
$x^{-i}\in\mathbb{R}^{d\times(N-1)}$ we introduce the probability measure
$\tau_{x^{-i}}\in\mathcal{P}_{2,\rm ac}(\mathbb{R}^{d})$ via

$\displaystyle\tau_{x^{-i}}\propto\exp\Bigl{(}-\frac{2}{\sigma^{2}}\,\delta\mathcal{F}(\rho_{x^{-i}},\cdot)\Bigr{)}\,.$

We can compute

$\displaystyle\operatorname{\mathsf{KL}}\bigl{(}\mu^{i\mid-i}(\cdot\mid
x^{-i})\bigm{\|}\tau_{x^{-i}}\bigr{)}$
$\displaystyle\qquad=\int\Bigl{(}\frac{2}{\sigma^{2}}\,\delta\mathcal{F}(\rho_{x^{-i}},\cdot)+\log\mu^{i\mid-i}(\cdot\mid
x^{-i})\Bigr{)}\,\mathrm{d}\mu^{i\mid-i}(\cdot\mid x^{-i})+\log
Z(\tau_{x^{-i}})\,,$

where $Z(\tau_{x^{-i}})$ is the normalization constant for $\tau_{x^{-i}}$,

$\displaystyle\log Z(\tau_{x^{-i}})$
$\displaystyle=\log\int\exp\Bigl{(}-\frac{2}{\sigma^{2}}\,\delta\mathcal{F}(\rho_{x^{-i}},z)\Bigr{)}\,\mathrm{d}z$
$\displaystyle=\log\int\exp\Bigl{(}-\log\pi(z)-\frac{2}{\sigma^{2}}\,\delta\mathcal{F}(\rho_{x^{-i}},z)\Bigr{)}\,\pi(\mathrm{d}z)$
$\displaystyle\geq-\frac{2}{\sigma^{2}}\int\delta\mathcal{F}(\rho_{x^{-i}},\cdot)\,\mathrm{d}\pi-\operatorname{\mathsf{ent}}(\pi)\,.$

Upon taking expectations, we obtain

$\displaystyle\mathsf{A}$
$\displaystyle\leq\frac{\sigma^{2}}{2}\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{-i}\sim\mu^{-i}}\operatorname{\mathsf{KL}}\bigl{(}\mu^{i\mid-i}(\cdot\mid
x^{-i})\bigm{\|}\tau_{x^{-i}}\bigr{)}\,.$

Moreover, we can recognize that $\tau_{x^{-i}}$ is a proximal Gibbs measure.
By Assumptions 6 and 7,

$\displaystyle\mathsf{A}$
$\displaystyle\leq\frac{\overline{C}_{\mathsf{LSI}}\,\sigma^{2}}{4}\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{-i}\sim\mu^{-i}}\operatorname{\mathsf{FI}}\bigl{(}\mu^{i\mid-i}(\cdot\mid
x^{-i})\bigm{\|}\tau_{x^{-i}}\bigr{)}$
$\displaystyle=\frac{\overline{C}_{\mathsf{LSI}}\,\sigma^{2}}{4}\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\Bigl{[}\Bigl{\lVert}\nabla_{x^{i}}\log\mu^{i\mid-i}(x^{i}\mid
x^{-i})+\frac{2}{\sigma^{2}}\,\nabla_{\mathcal{W}_{2}}\mathcal{F}(\rho_{x^{-i}},x^{i})\Bigr{\rVert}^{2}\Bigr{]}$
$\displaystyle=\frac{\overline{C}_{\mathsf{LSI}}\,\sigma^{2}}{4}\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\Bigl{[}\Bigl{\lVert}\nabla_{x^{i}}\log\mu^{1:N}(x^{1:N})+\frac{2}{\sigma^{2}}\,\nabla_{\mathcal{W}_{2}}\mathcal{F}(\rho_{x^{-i}},x^{i})\Bigr{\rVert}^{2}\Bigr{]}$
$\displaystyle=\frac{\overline{C}_{\mathsf{LSI}}}{\sigma^{2}}\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}[\lVert\nabla_{\mathcal{W}_{2}}\mathcal{F}(\rho_{x^{1:N}},x^{i})-\nabla_{\mathcal{W}_{2}}\mathcal{F}(\rho_{x^{-i}},x^{i})\rVert^{2}]$
$\displaystyle\leq\frac{\beta^{2}\overline{C}_{\mathsf{LSI}}}{\sigma^{2}}\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\mathcal{W}_{1}^{2}(\rho_{x^{1:N}},\rho_{x^{-i}})\,.$

To transport the mass from $\rho_{x^{1:N}}$ to $\rho_{x^{-i}}$, we take the
transport plan which moves $\frac{1}{N\,(N-1)}$ of the mass from $x^{i}$ to
each $x^{j}$, $j\neq i$. It yields

$\displaystyle\mathcal{W}_{1}(\rho_{x^{1:N}},\rho_{x^{-i}})$
$\displaystyle\leq\frac{1}{N\,(N-1)}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}{\lVert x^{i}-x^{j}\rVert}\,.$ (A.6)

Hence,

$\displaystyle\mathsf{A}$
$\displaystyle\leq\frac{\beta^{2}\overline{C}_{\mathsf{LSI}}}{\sigma^{2}N^{2}\,{(N-1)}^{2}}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\sum_{i=1}^{N}\Bigl{(}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}{\lVert x^{i}-x^{j}\rVert}\Bigr{)}^{2}$
$\displaystyle\leq\frac{\beta^{2}\overline{C}_{\mathsf{LSI}}}{\sigma^{2}N^{2}\,(N-1)}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\sum_{i\neq
j}{\lVert
x^{i}-x^{j}\rVert^{2}}=\frac{\beta^{2}\overline{C}_{\mathsf{LSI}}}{\sigma^{2}N}\operatorname{\mathbb{E}}_{x^{1:2}\sim\mu^{1:2}}[\lVert
x^{1}-x^{2}\rVert^{2}]\,.$

We then use the inequality

$\displaystyle\frac{1}{2}\operatorname{\mathbb{E}}_{x^{1:2}\sim\mu^{1:2}}[\lVert
x^{1}-x^{2}\rVert^{2}]$ $\displaystyle\leq
2\mathcal{W}_{2}^{2}(\mu^{1:2},\pi^{\otimes
2})+\operatorname{\mathbb{E}}_{x^{1:2}\sim\pi^{\otimes 2}}[\lVert
x^{1}-x^{2}\rVert^{2}]$
$\displaystyle\leq\frac{4}{N}\,\mathcal{W}_{2}^{2}(\mu^{1:N},\pi^{\otimes
N})+2\operatorname{\mathbb{E}}_{x\sim\pi}[\lVert
x-\operatorname{\mathbb{E}}x\rVert^{2}]$
$\displaystyle\leq\frac{8\overline{C}_{\mathsf{LSI}}}{N}\operatorname{\mathsf{KL}}(\mu^{1:N}\mathbin{\|}\pi^{\otimes
N})+2d\overline{C}_{\mathsf{LSI}}\,,$ (A.7)

where we used Lemma 23 and the Poincaré inequality for $\pi$. Hence,

$\displaystyle\mathsf{A}$
$\displaystyle\leq\frac{2\beta^{2}\overline{C}_{\mathsf{LSI}}}{\sigma^{2}N}\,\Bigl{(}\frac{8\overline{C}_{\mathsf{LSI}}}{N}\operatorname{\mathsf{KL}}(\mu^{1:N}\mathbin{\|}\pi^{\otimes
N})+2d\overline{C}_{\mathsf{LSI}}\Bigr{)}\,.$

Next, we turn toward term $\mathsf{B}$. First, define a function
$\zeta_{x^{1:N}}^{i}:\mathbb{R}^{d}\to\mathbb{R}$ by

$\displaystyle\zeta_{x^{1:N}}^{i}(y)\coloneqq\delta\mathcal{F}(\rho_{x^{-i}},y)-\delta\mathcal{F}(\rho_{x^{1:N}},y)\,.$

It is clear from Assumption 6 that this function is Lipschitz with constant
$2\beta\mathcal{W}_{1}(\rho_{x^{1:N}},\rho_{x^{-i}})$. Thus, we obtain using
this Lipschitzness, (A.6), and Young’s inequality,

$\displaystyle\mathsf{B}$
$\displaystyle=\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\int\bigl{(}\zeta_{x^{1:N}}^{i}(x^{i})-\zeta_{x^{1:N}}^{i}(z)\bigr{)}\,\pi(\mathrm{d}z)$
$\displaystyle\leq\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\int\frac{2\beta}{N\,(N-1)}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}{\lVert x^{j}-x^{i}\rVert\,\lVert
x^{i}-z\rVert}\,\pi(\mathrm{d}z)$
$\displaystyle\leq\frac{\beta}{N\,(N-1)}\operatorname{\mathbb{E}}_{x^{1:N}\sim\mu^{1:N}}\sum_{i\neq
j}{\lVert
x^{i}-x^{j}\rVert^{2}}+\frac{\beta}{N}\sum_{i=1}^{N}\operatorname{\mathbb{E}}_{(x^{i},z)\sim\mu^{1}\otimes\pi}[\lVert
x^{i}-z\rVert^{2}]$
$\displaystyle=\beta\operatorname{\mathbb{E}}_{x^{1:2}\sim\mu^{1:2}}[{\lVert
x^{1}-x^{2}\rVert^{2}}]+\beta\operatorname{\mathbb{E}}_{(x,z)\sim\mu^{1}\otimes\pi}[\lVert
x-z\rVert^{2}]\,.$

For the first term, we can apply (A.7), and for the second term, we can apply
(A.1). It yields

$\displaystyle\mathsf{B}$
$\displaystyle\leq\frac{20\beta\overline{C}_{\mathsf{LSI}}}{N}\operatorname{\mathsf{KL}}(\mu^{1:N}\mathbin{\|}\pi^{\otimes
N})+8\beta\overline{C}_{\mathsf{LSI}}d\,.$

Putting the bounds together with Lemma 19,

$\displaystyle\operatorname{\mathsf{KL}}(\mu^{1:N}\mathbin{\|}\pi^{\otimes
N})$ $\displaystyle\leq\frac{33\beta\overline{C}_{\mathsf{LSI}}d}{\sigma^{2}}$

for all $N\geq 160\beta\overline{C}_{\mathsf{LSI}}/\sigma^{2}$. The result for
$k\leq N$ follows from Lemma 24. ∎

##  Appendix B Isoperimetric Results for the Stationary Distributions

###  B.1 Convexity and Smoothness

Here, we verify the convexity and smoothness properties of $\mu^{1:N}$ in the
pairwise McKean–Vlasov setting.

Proof of Lemma 7 For $x^{1:N}=[x^{1},\dotsc,x^{N}]\in\mathbb{R}^{d\times N}$,
the Hessian of $\log(1/\mu^{1:N})$ can be explicitly computed as

$\displaystyle-\frac{\sigma^{2}}{2}\,\nabla^{2}\log\mu^{1:N}(x^{1:N})=\begin{bmatrix}\nabla^{2}V(x^{1})&0&\cdots&0\\\
0&\nabla^{2}V(x^{2})&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\nabla^{2}V(x^{N})\end{bmatrix}$
$\displaystyle\qquad+\frac{1}{N-1}\,\underbrace{\begin{bmatrix}\sum_{j=2}^{N}\nabla^{2}W(x^{1}-x^{j})&-\nabla^{2}W(x^{1}-x^{2})&\cdots&-\nabla^{2}W(x^{1}-x^{N})\\\
-\nabla^{2}W(x^{2}-x^{1})&\sum_{\begin{subarray}{c}j=1\\\ j\neq
2\end{subarray}}^{N}\nabla^{2}W(x^{2}-x^{j})&\cdots&-\nabla^{2}W(x^{2}-x^{N})\\\
\vdots&\vdots&\ddots&\vdots\\\
-\nabla^{2}W(x^{N}-x^{1})&-\nabla^{2}W(x^{N}-x^{2})&\cdots&\sum_{j=1}^{N-1}\nabla^{2}W(x^{N}-x^{j})\end{bmatrix}}_{=\mathsf{B}}\,.$

Clearly, the first block matrix has eigenvalues between $\alpha_{V}$ and
$\beta_{V}$. For the second block matrix $\mathsf{B}$, let us denote
$A_{i,j}\coloneqq\nabla^{2}W(x^{i}-x^{j})$ for $i,j\in[N]$. Note that
$A_{i,j}=A_{j,i}$ since $W$ is even, and each $A_{i,j}$ is clearly symmetric.

For $\mathsf{B}\in\mathbb{R}^{dN\times dN}$ the second matrix and
$y=[y^{1},\dotsc,y^{N}]\in\mathbb{R}^{dN}$, we have

$\displaystyle y^{\mathsf{T}}\mathsf{B}y$ $\displaystyle=\sum_{i\leq
N}y_{i}^{\mathsf{T}}\,\Bigl{(}\sum_{j\in[N]\backslash
i}A_{i,j}\Bigr{)}\,y_{i}-\sum_{\begin{subarray}{c}i,j\leq N\\\ i\neq
j\end{subarray}}y_{i}^{\mathsf{T}}A_{i,j}y_{j}$
$\displaystyle=\sum_{\begin{subarray}{c}i,j\leq N\\\
i<j\end{subarray}}\bigl{(}y_{i}^{\mathsf{T}}A_{i,j}y_{i}+y_{j}^{\mathsf{T}}A_{j,i}y_{j}-y_{i}^{\mathsf{T}}A_{i,j}y_{j}-y_{j}^{\mathsf{T}}A_{j,i}y_{i}\bigr{)}=\sum_{\begin{subarray}{c}i,j\leq
N\\\ i<j\end{subarray}}(y_{i}-y_{j})^{\mathsf{T}}A_{i,j}(y_{i}-y_{j})\,.$

Using $\alpha_{W}I_{d}\preceq A_{i,j}\preceq\beta_{W}I_{d}$ and

$\displaystyle\mathsf{M}\coloneqq\nabla^{2}_{y}\sum_{\begin{subarray}{c}i,j\leq
N\\\ i<j\end{subarray}}\lVert y_{i}-y_{j}\rVert^{2}$
$\displaystyle=2\,\begin{bmatrix}N-1&-1&\cdots&-1\\\ -1&N-1&\cdots&-1\\\
\vdots&\vdots&\ddots&\vdots\\\ -1&-1&\cdots&N-1\end{bmatrix}\otimes I_{d}\,,$

we have
$\frac{1}{2}\,\alpha_{W}\mathsf{M}\preceq\mathsf{B}\preceq\frac{1}{2}\,\beta_{W}\mathsf{M}$.
Since the circulant matrix in $\mathsf{M}$ is PSD due to diagonal dominance
and its largest eigenvalue is at most $N$, it follows that the eigenvalues of
$\mathsf{M}$ lie between $0$ and $2N$. Hence, the eigenvalues of $\mathsf{B}$
lie in the interval
$[\frac{N}{N-1}\,{\alpha_{W}^{-}},\frac{N}{N-1}\,\beta_{W}]$. ∎

###  B.2 Bounded Perturbations

In this section, we prove the isoperimetric results from §3.2.1. We again
introduce the conditional measure: if
$x^{-i}\coloneqq(x^{1:i-1},x^{i+1:N})\in\mathbb{R}^{d\times(N-1)}$ we define

$\mu^{1:N}(x^{1:N})=\mu^{i|-i}(x^{i}\mid x^{-i})\times\mu^{-i}(x^{-i})$

for the conditional distribution of the $i$-th particle and the distribution
of an $N$-particle system with the $i$-th particle marginalized out.

Proof of Lemma 8 We begin by proving the statement about $\pi$. The potential
of the invariant measure $\pi$ can also be written as

$\displaystyle\log\frac{1}{\pi(x)}$
$\displaystyle=\frac{2}{\sigma^{2}}\,\Bigl{(}V_{0}(x)+V_{1}(x)+\int\bigl{(}W_{0}(x-\cdot)+W_{1}(x-\cdot)\bigr{)}\,\mathrm{d}\pi\Bigr{)}$
$\displaystyle=\frac{2}{\sigma^{2}}\,\Bigl{(}V_{0}(x)+\int
W_{0}(x-\cdot)\,\mathrm{d}\pi\Bigr{)}+\frac{2}{\sigma^{2}}\,\Bigl{(}V_{1}(x)+\int
W_{1}(x-\cdot)\,\mathrm{d}\pi\Bigr{)}\,.$

This is the sum of a
$\frac{2}{\sigma^{2}}\,(\alpha_{V_{0}}+\alpha_{W_{0}})$-convex function with a
$\frac{2}{\sigma^{2}}\,(\operatorname{\mathsf{osc}}(V_{1})+\operatorname{\mathsf{osc}}(W_{1}))$-bounded
perturbation. Thus, $\pi$ satisfies (LSI) with the claimed parameter.

We now prove the statement about $\mu^{1:N}$. By the tensorization argument in
[BGL14, Proposition 5.2.7] or [BLM13, Theorem 4.10], we obtain the inequality
that for any smooth test function $g^{2}:\mathbb{R}^{d\times
N}\to\mathbb{R}_{+}$ and $X^{-i}\sim\mu^{-i}$,

$\operatorname{\mathsf{ent}}_{\mu^{1:N}}(g^{2})\leq\sum_{i=1}^{N}\operatorname{\mathbb{E}}\operatorname{\mathsf{ent}}_{\mu^{i|-i}(\cdot|X^{-i})}\bigl{(}g^{2}(X^{1:i-1},\cdot,X^{i+1:N})\bigr{)}\,,$
(B.1)

where each entropy term in RHS involves the distribution of the $i$-th
particle, conditioned on the remaining particles. This conditional measure has
a density of the form

$\displaystyle\mu^{i|-i}(\cdot\mid
x^{-i})\propto\exp\Bigl{(}-\frac{2}{\sigma^{2}}\,V(\cdot)-\frac{2}{\sigma^{2}\,(N-1)}\sum_{j\in[N]\backslash
i}W(\cdot-x^{j})\Bigr{)}\,,$

where both $V$ and $W$ are bounded perturbations of $\alpha_{V_{0}}$,
$\alpha_{W_{0}}$-strongly convex functions respectively, irrespective of the
conditional variables. Thus, by Holley–Stroock perturbation and the
Bakry–Émery condition, each $\mu^{i|-i}(\cdot\mid x^{-i})$ satisfies (LSI)
with parameter

$\displaystyle
C_{\mathsf{LSI}}\leq\frac{\sigma^{2}}{2}\,\bigl{(}\alpha_{V_{0}}+\frac{N}{N-1}\,{\alpha_{W_{0}}^{-}}\bigr{)}^{-1}\exp\Bigl{(}\frac{2}{\sigma^{2}}\,\bigl{(}\operatorname{\mathsf{osc}}(V_{1})+\operatorname{\mathsf{osc}}(V_{1})\bigr{)}\Bigr{)}\,.$

Therefore, we can further bound each entropy term in (B.1) by

$\displaystyle\operatorname{\mathsf{ent}}_{\mu^{1:N}}(g^{2})\leq\sum_{i=1}^{N}2C_{\mathsf{LSI}}\operatorname{\mathbb{E}}\operatorname{\mathbb{E}}_{Z\sim\mu^{i|-i}(\cdot\mid
X^{-i})}[\lvert\partial_{i}g(X^{1:i-1},Z,X^{i+1:N})\rvert^{2}]=2C_{\mathsf{LSI}}\operatorname{\mathbb{E}}_{\mu^{1:N}}[\lVert\nabla
g\rVert^{2}]\,.$

Therefore, $\mu^{1:N}$ satisfies (LSI) with parameter $C_{\mathsf{LSI}}$. ∎

###  B.3 Logarithmic Sobolev Inequalities via Perturbations

In this section, we state log-Sobolev inequalities for Lipschitz perturbations
of strongly log-concave measures, which is used for the general McKean–Vlasov
setting in §2.1.2.

######  Lemma 20 (LSI under Lipschitz Perturbations [CG22, Theorem 2.7]).

Let $\nu\propto\exp(-V)$ for an $\alpha_{V}$-strongly convex and
$\beta_{V}$-smooth function $V:\mathbb{R}^{d}\to\mathbb{R}$. Let
$H:\mathbb{R}^{d}\to\mathbb{R}$ be a $L$-Lipschitz continuous function. Then,
$\mu\propto\exp(-H)\,\nu$ satisfies (LSI) with constant
$C_{\mathsf{LSI}}(\mu)$ given by

$\displaystyle
C_{\mathsf{LSI}}(\mu)\leq\frac{4}{\alpha_{V}}+\Bigl{(}\frac{L}{\alpha_{V}}+\sqrt{\frac{2}{\alpha_{V}}}\Bigr{)}^{2}\,\Bigl{(}2+d+\frac{d}{2}\log\frac{\beta_{V}}{\alpha_{V}}+\frac{4L^{2}}{\alpha_{V}}\Bigr{)}\exp\Bigl{(}\frac{L^{2}}{2\alpha_{V}}\Bigr{)}\,.$

From this one derives the following log-Sobolev inequality for the proximal
Gibbs measure.

######  Lemma 21 (Uniform LSI for the Proximal Gibbs Measure [SWN23, Theorem
1]).

For the proximal Gibbs measure (2.3) in the setting of §2.1.2, under
Assumptions 6, and 8, we have that

$\displaystyle\sup_{\mu\in\mathcal{P}_{2,\rm
ac}(\mathbb{R}^{d})}C_{\mathsf{LSI}}(\pi_{\mu})\leq\overline{C}_{\mathsf{LSI}}\,,$

where $\alpha$ can be bounded by

$\displaystyle\overline{C}_{\mathsf{LSI}}\leq\frac{\sigma^{2}}{\beta}\exp\Bigl{(}\frac{8B^{2}}{\lambda\sigma^{2}}\sqrt{2d/\pi}\Bigr{)}\wedge\Bigl{\\{}\frac{2\sigma^{2}}{\beta}+\bigl{(}\frac{B}{\beta}+\frac{\sigma}{\sqrt{\beta}}\bigr{)}^{2}\,\bigl{(}2+d+\frac{8B^{2}}{\beta\sigma^{2}}\bigr{)}\exp\bigl{(}\frac{B^{2}}{\lambda\sigma^{2}}\bigr{)}\Bigr{\\}}\,.$

Finally, we provide an LSI for the finite-particle stationary distribution
$\mu^{1:N}$ with a constant independent of $N$, using an approximate
tensorization argument.

######  Lemma 22.

Consider the general McKean–Vlasov setting with
$\mathcal{F}(\mu)=\mathcal{F}_{0}(\mu)+\frac{\lambda}{2}\int\lVert\cdot\rVert^{2}\,\mathrm{d}\mu$
and suppose that Assumption 8 holds. Then, $\mu^{1:N}$ satisfies (LSI) with
parameter independent of $N$ given by

$\displaystyle
C_{\mathsf{LSI}}(\mu^{1:N})\leq\frac{4\sigma^{2}}{\lambda}+\Bigl{(}\frac{2B}{\lambda}+\frac{\sigma}{\sqrt{\lambda}}\Bigr{)}^{2}\,\Bigl{(}2+d+\frac{16B^{2}}{\lambda\sigma^{2}}\Bigr{)}\exp\Bigl{(}\frac{2B^{2}}{\lambda\sigma^{2}}\Bigr{)}\,.$

Proof. Since
$\mu^{1:N}(x^{1:N})\propto\exp(-\frac{2N}{\sigma^{2}}\,\mathcal{F}(\rho_{x^{1:N}}))$,
where $\rho_{x^{1:N}}=\frac{1}{N}\sum_{i=1}^{N}\delta_{x^{i}}$ is the
empirical measure, we can check that

$\mu^{i\mid-i}(x^{i}\mid
x^{-i})\propto\exp\Bigl{(}-\frac{2N}{\sigma^{2}}\,\mathcal{F}_{0}(\rho_{x^{1:N}})-\frac{\lambda}{\sigma^{2}}\,\lVert
x^{i}\rVert^{2}\Bigr{)}\,.$

Also, note that

$\displaystyle\nabla_{x^{i}}\mathcal{F}_{0}(\rho_{x^{1:N}})$
$\displaystyle=\frac{1}{N}\,\nabla_{\mathcal{W}_{2}}\mathcal{F}_{0}(\rho_{x^{1:N}},x^{i})\,.$

Using Assumption 8,
$x^{i}\mapsto\frac{2N}{\sigma^{2}}\,\mathcal{F}_{0}(\rho_{x^{1:N}})$ is
$\nicefrac{{2B}}{{\sigma^{2}}}$-Lipschitz.

Then, apply Lemma 20 with
$\alpha_{V}=\beta_{V}=\nicefrac{{\lambda}}{{\sigma^{2}}}$, to obtain

$\displaystyle C_{\mathsf{LSI}}\bigl{(}\mu^{i|-i}(\cdot\mid
x^{-i})\bigr{)}\leq\frac{4\sigma^{2}}{\lambda}+\Bigl{(}\frac{2B}{\lambda}+\frac{\sigma}{\sqrt{\lambda}}\Bigr{)}^{2}\,\Bigl{(}2+d+\frac{16B^{2}}{\lambda\sigma^{2}}\Bigr{)}\exp\Bigl{(}\frac{2B^{2}}{\lambda\sigma^{2}}\Bigr{)}\,.$

Applying the tensorization argument from the proof of Lemma 8 completes the
proof. ∎

##  Appendix C Explicit Calculations for the Gaussian Case

Here we provide complete details for Example 10: for any $k\leq N$,

$\frac{dk^{2}}{N^{2}}\lesssim\operatorname{\mathsf{KL}}(\mu^{1:k}\mathbin{\|}\pi^{\otimes
k})\lesssim\frac{dk^{2}}{N^{2}}\log N\,.$

Note that for $\mathsf{C}\in\mathbb{R}^{N\times N}$ with
$\mathsf{C}_{i,i}=N-1$ and $\mathsf{C}_{i,j}=-1$ if $i\neq j$,

$\mu^{1:N}=\mathcal{N}\Bigl{(}0,\underbrace{\frac{\sigma^{2}}{2}\,\bigl{(}I_{N}\otimes
A+\frac{\lambda}{N-1}\,\mathsf{C}\otimes
I_{d}\bigr{)}^{-1}}_{\eqqcolon\Sigma_{1}}\Bigr{)}\qquad\text{and}\qquad\pi=\mathcal{N}\Bigl{(}0,\underbrace{\frac{\sigma^{2}}{2}\,{(A+\lambda
I_{d})}^{-1}}_{\eqqcolon\Sigma_{2}}\Bigr{)}\,.$

The $k$-particle marginal $\mu^{1:k}$ is a Gaussian with zero mean and
covariance being the upper-left $(kN\times kN)$-block matrix of $\Sigma_{1}$,
which we denote by $\Sigma_{1,k}$. Clearly, $\pi^{\otimes k}$ is also a
Gaussian with zero mean and covariance $\Sigma_{2,k}\coloneqq
I_{k}\otimes\Sigma_{2}$. From a well-known formula for the
$\operatorname{\mathsf{KL}}$ divergence between two Gaussian distributions,

$\operatorname{\mathsf{KL}}(\mu^{1:k}\mathbin{\|}\pi^{\otimes
k})=\frac{1}{2}\,\bigl{(}-\log\det(\Sigma_{2,k}^{-1}\Sigma_{1,k})-dk+\mathop{\mathrm{missing}}{tr}(\Sigma_{2,k}^{-1}\Sigma_{1,k})\bigr{)}\,.$
(C.1)

Let $1_{p}\in\mathbb{R}^{p}$ be the $p$-dimensional vector with all entries
$1$. From $\mathsf{C}=NI_{N}-1_{N}1_{N}^{\mathsf{T}}$,

$\displaystyle\frac{2}{\sigma^{2}}\,\Sigma_{1}$
$\displaystyle=\Bigl{(}I_{N}\otimes\bigl{(}\underbrace{A+\frac{\lambda
N}{N-1}\,I_{d}}_{\eqqcolon
A_{\lambda}}\bigr{)}-\frac{\lambda}{N-1}\,(1_{N}1_{N}^{\mathsf{T}})\otimes
I_{d}\Bigr{)}^{-1}$ $\displaystyle\underset{\text{(i)}}{=}\Bigl{(}I_{N}\otimes
A_{\lambda}-\frac{\lambda}{N-1}\,(1_{N}\otimes
I_{d})(1_{N}^{\mathsf{T}}\otimes I_{d})\Bigr{)}^{-1}$
$\displaystyle\underset{\text{(ii)}}{=}\mathopen{}\mathclose{{}\left(I_{N}\otimes
A_{\lambda}}\right)^{-1}$ $\displaystyle\qquad{}-{(I_{N}\otimes
A_{\lambda})}^{-1}(1_{N}\otimes I_{d})\bigl{(}I_{d}+(1_{N}^{\mathsf{T}}\otimes
I_{d})\mathopen{}\mathclose{{}\left(I_{N}\otimes
A_{\lambda}}\right)^{-1}(1_{N}\otimes
I_{d})\bigr{)}^{-1}(1_{N}^{\mathsf{T}}\otimes
I_{d})\mathopen{}\mathclose{{}\left(I_{N}\otimes A_{\lambda}}\right)^{-1}$
$\displaystyle\underset{\text{(iii)}}{=}I_{N}\otimes
A_{\lambda}^{-1}-(1_{N}\otimes
A_{\lambda}^{-1})\bigl{(}I_{d}+(1_{N}^{\mathsf{T}}1_{N})\otimes
A_{\lambda}^{-1}\bigr{)}^{-1}(1_{N}^{\mathsf{T}}\otimes A_{\lambda}^{-1})$
$\displaystyle=I_{N}\otimes A_{\lambda}^{-1}-(1_{N}\otimes
A_{\lambda}^{-1}){(I_{d}+NA_{\lambda}^{-1})}^{-1}(1_{N}^{\mathsf{T}}\otimes
A_{\lambda}^{-1})$ $\displaystyle=I_{N}\otimes
A_{\lambda}^{-1}-(1_{N}1_{N}^{\mathsf{T}})\otimes{(A_{\lambda}^{2}+NA_{\lambda})}^{-1}\,,$

where in (i) we used $(A\otimes B)(C\otimes D)=(AC)\otimes(BD)$, (ii) follows
from the Woodbury matrix identity, and (iii) used $(A\otimes
B)^{-1}=A^{-1}\otimes B^{-1}$. Hence it follows that

$\frac{2}{\sigma^{2}}\,\Sigma_{1,k}=I_{k}\otimes
A_{\lambda}^{-1}-(1_{k}1_{k}^{\mathsf{T}})\otimes{(A_{\lambda}^{2}+NA_{\lambda})}^{-1}\,.$

By the spectral decomposition of $A$, we can write $A=UDU^{\mathsf{T}}$ for a
diagonal $D\in\mathbb{R}^{d\times d}$ and an orthogonal matrix
$U\in\mathbb{R}^{d\times d}$ such that $\\{\sigma_{i}\coloneqq
D_{i,i}\\}_{i\in[d]}$ correspond to the eigenvalues of $A$. Since
$\log\det(\cdot)$ and $\mathop{\mathrm{missing}}{tr}(\cdot)$ in (C.1) are
orthogonally invariant, let us look at the orthogonal conjugate of
$\Sigma_{2,k}^{-1}\Sigma_{1,k}$ by $I_{k}\otimes
U^{\mathsf{T}}\in\mathbb{R}^{dk\times dk}$. Using $(A\otimes
B)^{\mathsf{T}}=A^{\mathsf{T}}\otimes B^{\mathsf{T}}$ and denoting
$D_{\lambda}\coloneqq D+\frac{\lambda N}{N-1}\,I_{d}$,

$\displaystyle(I_{k}\otimes
U^{\mathsf{T}})\Sigma_{2,k}^{-1}\Sigma_{1,k}(I_{k}\otimes U)$
$\displaystyle\qquad=(I_{k}\otimes
U^{\mathsf{T}})\bigl{(}I_{k}\otimes(A+\lambda
I_{d})\bigr{)}\bigl{(}I_{k}\otimes
A_{\lambda}^{-1}-(1_{k}1_{k}^{\mathsf{T}})\otimes{(A_{\lambda}^{2}+NA_{\lambda})}^{-1}\bigr{)}(I_{k}\otimes
U)$ $\displaystyle\qquad=\bigl{(}I_{k}\otimes(D+\lambda
I_{d})\bigr{)}\bigl{(}I_{k}\otimes
D_{\lambda}^{-1}-(1_{k}1_{k}^{\mathsf{T}})\otimes{(D_{\lambda}^{2}+ND_{\lambda})}^{-1}\bigr{)}$
$\displaystyle\qquad=I_{k}\otimes\mathopen{}\mathclose{{}\left((D+\lambda
I_{d})D_{\lambda}^{-1}}\right)-\underbrace{(1_{k}1_{k}^{\mathsf{T}})}_{\eqqcolon
J_{k}}\otimes\underbrace{\bigl{(}(D+\lambda
I_{d}){(D_{\lambda}^{2}+ND_{\lambda})}^{-1}\bigr{)}}_{\eqqcolon S_{\lambda}}$
$\displaystyle\qquad=I_{dk}-\Bigl{(}\underbrace{\frac{\lambda}{N-1}\,I_{k}\otimes
D_{\lambda}^{-1}+J_{k}\otimes S_{\lambda}}_{\eqqcolon\mathsf{M}}\Bigr{)}\,.$

For $\sigma_{d}\coloneqq\min_{i\in[d]}\sigma_{i}$,
$\alpha\coloneqq\sigma_{d}+\lambda$, and $\varepsilon\coloneqq\lambda/(N-1)$,
we have $D_{\lambda}^{-1}\precsim\frac{1}{\alpha}\,I_{d}$ and
$S_{\lambda}\precsim\frac{1}{\alpha+N}\,I_{d}$ due to

$((D_{\lambda})^{-1})_{i,i}\leq\frac{1}{\alpha+\varepsilon}\qquad\text{and}\qquad(S_{\lambda})_{i,i}\leq\frac{\alpha}{(\alpha+\varepsilon)\,(\alpha+\varepsilon+N)}\,.$

Since the eigenvalues of $A\otimes B$ consist of all possible combinations
arising from the product of eigenvalues, one from $A$ and one from $B$, the
largest eigenvalue $\eta_{1}$ of $\mathsf{M}$ is less than $1$:

$\eta_{1}\leq\frac{\lambda}{N-1}\,\lVert D_{\lambda}^{-1}\rVert+k\,\lVert
S_{\lambda}\rVert\leq\frac{\varepsilon}{\alpha+\varepsilon}+\frac{\alpha
N}{(\alpha+\varepsilon)\,(\alpha+\varepsilon+N)}=\frac{\varepsilon+N}{\alpha+\varepsilon+N}\,.$

Denoting the eigenvalues of $\mathsf{M}$ by $\eta_{i}$, it follows from (C.1)
that

$\displaystyle 2\operatorname{\mathsf{KL}}(\mu^{1:k}\mathbin{\|}\pi^{\otimes
k})$
$\displaystyle=-\bigl{(}\log\det(I_{dk}-\mathsf{M})+dk-\mathop{\mathrm{missing}}{tr}(I_{dk}-\mathsf{M})\bigr{)}=-\sum_{i=1}^{dk}\bigl{(}\log(1-\eta_{i})+\eta_{i}\bigr{)}$
$\displaystyle=\sum_{i=1}^{dk}\sum_{n\geq 2}\frac{\eta_{i}^{n}}{n}\,.$

Then, we have a trivial lower bound of
$\frac{1}{2}\sum_{i=1}^{dk}\eta_{i}^{2}$, and as for the upper bound,

$\sum_{i=1}^{dk}\sum_{n\geq
2}\frac{\eta_{i}^{n}}{n}\leq\sum_{i=1}^{dk}\Bigl{(}\eta_{i}^{2}+\sum_{n\geq
1}\frac{\eta_{i}^{n+2}}{n}\Bigr{)}=\sum_{i=1}^{dk}\eta_{i}^{2}\log\bigl{(}\frac{e}{1-\eta_{i}}\bigr{)}\lesssim\bigl{(}1\vee\log\frac{N}{\alpha}\bigr{)}\mathop{\mathrm{missing}}{tr}(\mathsf{M}^{2})\,,$

where the last inequality follows from
$(1-\eta_{i})^{-1}\leq(1-\eta_{1})^{-1}$.

Using $\mathop{\mathrm{missing}}{tr}(A\otimes
B)=\mathop{\mathrm{missing}}{tr}(A)\cdot\mathop{\mathrm{missing}}{tr}(B)$, and
$D_{\lambda}^{-1}\precsim\frac{1}{\alpha}\,I_{d}$ and
$S_{\lambda}\precsim\frac{1}{\alpha+N}\,I_{d}$, we have

$\displaystyle\mathop{\mathrm{missing}}{tr}(\mathsf{M}^{2})$
$\displaystyle\lesssim\frac{\lambda^{2}}{{(N-1)}^{2}}\mathop{\mathrm{missing}}{tr}(I_{k})\mathop{\mathrm{missing}}{tr}(D_{\lambda}^{-2})+\mathop{\mathrm{missing}}{tr}(J_{k}^{2})\mathop{\mathrm{missing}}{tr}(S_{\lambda}^{2})$
$\displaystyle\lesssim\frac{\lambda^{2}}{\alpha^{2}}\,\frac{dk}{N^{2}}+\frac{dk^{2}}{{(\alpha+N)}^{2}}\lesssim\frac{dk^{2}}{N^{2}}\,.$

As for the lower bound, since $(S_{\lambda})_{i,i}\sim\frac{1}{N}$ for large
$N$, we have

$\mathop{\mathrm{missing}}{tr}(\mathsf{M}^{2})\gtrsim\frac{dk^{2}}{N^{2}}\,,$

which completes the proof.

##  Appendix D Additional Technical Lemmas

In our proofs, we used the following general lemmas on exchangeability.

######  Lemma 23.

Let $\mu^{1:N}$, $\nu^{1:N}$ be two exchangeable measures over
$\mathbb{R}^{d\times N}$. For any $k\leq N$,

$\displaystyle\mathcal{W}_{2}^{2}(\mu^{1:k},\nu^{1:k})$
$\displaystyle\leq\frac{k}{N}\,\mathcal{W}_{2}^{2}(\mu^{1:N},\nu^{1:N})\,.$

Proof. Let $(X^{1:N},Y^{1:N})$ be optimally coupled for $\mu^{1:N}$ and
$\nu^{1:N}$. For each subset $S\subseteq[N]$ of size $k$, it induces a
coupling $(X^{S},Y^{S})$ of $\mu^{1:k}$ and $\nu^{1:k}$ (by exchangeability).
In particular, the law of $(X^{\mathsf{S}},Y^{\mathsf{S}})$, where
$\mathsf{S}$ is an independent and uniformly random subset of size $k$, is
also a coupling of $\mu^{1:k}$ and $\nu^{1:k}$. Hence,

$\displaystyle\mathcal{W}_{2}^{2}(\mu^{1:k},\nu^{1:k})$
$\displaystyle\leq\operatorname{\mathbb{E}}[\lVert
X^{\mathsf{S}}-Y^{\mathsf{S}}\rVert^{2}]=\frac{1}{\binom{N}{k}}\sum_{\lvert
S\rvert=k}\operatorname{\mathbb{E}}[\lVert X^{S}-Y^{S}\rVert^{2}]$
$\displaystyle=\frac{1}{\binom{N}{k}}\sum_{i=1}^{N}\sum_{\lvert
S\rvert=k\,:\,i\in S}\operatorname{\mathbb{E}}[\lVert
X^{i}-Y^{i}\rVert^{2}]=\frac{\binom{N-1}{k-1}}{\binom{N}{k}}\sum_{i=1}^{N}\operatorname{\mathbb{E}}[\lVert
X^{i}-Y^{i}\rVert^{2}]$
$\displaystyle=\frac{k}{N}\,\mathcal{W}_{2}^{2}(\mu^{1:N},\nu^{1:N})\,,$

which completes the proof. ∎

######  Lemma 24 (Information Inequality [Csi84]).

If $\mathcal{X}^{1},\ldots,\mathcal{X}^{N}$ are Polish spaces and $\mu^{1:N}$,
$\nu^{1:N}$ are probability measures on
$\mathcal{X}^{1}\times\cdots\times\mathcal{X}^{N}$, where
$\nu^{1:N}=\nu^{1}\otimes\cdots\otimes\nu^{N}$ is a product measure, then for
the marginals $\mu^{i}$ of $\mu$, it holds that

$\sum_{i=1}^{N}\operatorname{\mathsf{KL}}(\mu^{i}\mathbin{\|}\nu^{i})\leq\operatorname{\mathsf{KL}}(\mu^{1:N}\mathbin{\|}\nu^{1:N})\,.$

In particular when $\mu^{1:N}$, $\nu^{1:N}$ are both exchangeable, this states
that
$\operatorname{\mathsf{KL}}(\mu^{1}\mathbin{\|}\nu^{1})\leq\frac{1}{N}\operatorname{\mathsf{KL}}(\mu^{1:N}\mathbin{\|}\nu^{1:N})$.

Note that Lemma 24 follows from the chain rule and convexity of the
$\operatorname{\mathsf{KL}}$ divergence.

##  Appendix E Sampling Guarantees

Here, we show how to obtain the claimed rates in §4. We begin with some
preliminary facts.

#####  $\operatorname{\mathsf{KL}}$ divergence guarantees.

To obtain our guarantees in $\operatorname{\mathsf{KL}}$ divergence, we use
the following lemma.

######  Lemma 25 ([Zha+23, Proof of Theorem 6]).

Let $\hat{\mu}$, $\mu$, and $\pi$ be three probability measures, and assume
that $\mu$ satisfies (LSI) with constant $C_{\mathsf{LSI}}(\mu)$. Then,

$\displaystyle\operatorname{\mathsf{KL}}(\hat{\mu}\mathbin{\|}\pi)$
$\displaystyle\leq
2\,\chi^{2}(\hat{\mu}\mathbin{\|}\mu)+\operatorname{\mathsf{KL}}(\mu\mathbin{\|}\pi)+\frac{C_{\mathsf{LSI}}(\mu)}{4}\operatorname{\mathsf{FI}}(\mu\mathbin{\|}\pi)\,.$

We instantiate the lemma with $\hat{\mu}^{1:N}$, $\mu^{1:N}$, and
$\pi^{\otimes N}$ respectively. In the setting of Theorem 4, it is seen that
$\operatorname{\mathsf{KL}}(\mu^{1:N}\mathbin{\|}\pi^{\otimes N})$ and
$C_{\mathsf{LSI}}(\mu^{1:N})\operatorname{\mathsf{FI}}(\mu^{1:N}\mathbin{\|}\pi^{\otimes
N})$ are of the same order and can be made at most $N\varepsilon^{2}$ if we
take $N\asymp\kappa^{4}d/\varepsilon^{2}$. Thus, if we have a sampler that
achieves $\chi^{2}(\hat{\mu}^{1:N}\mathbin{\|}\mu^{1:N})\leq
N\varepsilon^{2}$, it follows from exchangeability (Lemma 24) that
$\operatorname{\mathsf{KL}}(\hat{\mu}^{1}\mathbin{\|}\pi)\lesssim\varepsilon^{2}$.

##### Guarantees using the sharp propagation of chaos bound.

Here, we impose Assumptions 1, 2, 3, and 9. It follows from Theorem 3 that
$N=\widetilde{\Theta}(\sqrt{d}/\varepsilon)$ suffices in order to make
$\nicefrac{{\sqrt{\alpha}}}{{\sigma}}\,\mathcal{W}_{2}(\mu^{1},\pi)\leq\varepsilon$.
For the first term, we use exchangeability (Lemma 23) to argue that

$\displaystyle\mathcal{W}_{2}(\hat{\mu}^{1},\mu^{1})\leq
N^{-1/2}\,\mathcal{W}_{2}(\hat{\mu}^{1:N},\mu^{1:N})\,,$

and hence we invoke sampling guarantees to ensure that
$\nicefrac{{\sqrt{\alpha}}}{{\sigma}}\,\mathcal{W}_{2}(\hat{\mu}^{1:N},\mu^{1:N})\leq
N^{1/2}\varepsilon$ under (LSI).

  * •

LMC: We use the guarantee for Langevin Monte Carlo from [VW19].

  * •

MALA–PS: We use the guarantee for the Metropolis-adjusted Langevin algorithm
together with the proximal sampler from [AC23]. Note that the iteration
complexity is $\widetilde{\mathcal{O}}(\kappa d^{1/2}N^{1/2})$, and we
substitute in the chosen value for $N$.

  * •

ULMC–PS: Here, we use underdamped Langevin Monte Carlo to implement the
proximal sampler. To justify the sampling guarantee, note that since
$\log\mu^{1:N}$ is $\beta$-smooth, if we choose step size $h=\frac{1}{2\beta}$
for the proximal sampler, then the RGO is $\beta$-strongly log-concave and
$3\beta$-log-smooth. According to [AC23, Proof of Theorem 5.3], it suffices to
implement the RGO in each iteration to accuracy
$N^{1/2}\varepsilon/\kappa^{1/2}$ in $\sqrt{\operatorname{\mathsf{KL}}}$.
Then, from [Zha+23], this can be done via ULMC with complexity
$\widetilde{\mathcal{O}}(\kappa^{1/2}d^{1/2}/\varepsilon)$. Finally, since the
number of outer iterations of the proximal sampler is
$\widetilde{\mathcal{O}}(\kappa)$, we obtain the claim.

  * •

ULMC+: Here, we use either the randomized midpoint discretization [SL19] or
the shifted ODE discretization [FLO21] of the underdamped Langevin diffusion.
We also replace the LSI assumptions (Assumptions 2 and 9) with strong
convexity (Assumption 4).

##### Guarantees under strong displacement convexity.

Here, we impose Assumptions 1 and 4. As discussed above, to obtain
$\operatorname{\mathsf{KL}}$ guarantees, we require log-concave samplers that
can achieve $\chi^{2}(\hat{\mu}^{1:N}\mathbin{\|}\mu^{1:N})\leq
N\varepsilon^{2}$. For $\mathcal{W}_{2}$ guarantees, by Theorem 4 we take
$N\asymp\kappa^{2}d/\varepsilon^{2}$ and we require log-concave samplers that
can achieve
$\nicefrac{{\sqrt{\alpha}}}{{\sigma}}\,\mathcal{W}_{2}(\hat{\mu}^{1:N},\mu^{1:N})\leq
N^{1/2}\varepsilon$.

  * •

LMC: For Langevin Monte Carlo, we use the $\chi^{2}$ guarantee from [Che+22a]
and the $\mathcal{W}_{2}$ guarantee from [DMM19].

  * •

ULMC: For underdamped Langevin Monte Carlo, we use the $\chi^{2}$ guarantee
from [AC23].

  * •

ULMC+: Here, we use the $\mathcal{W}_{2}$ guarantees for either the randomized
midpoint discretization [SL19] or the shifted ODE discretization [FLO21] of
the underdamped Langevin diffusion.

##### Guarantees in the general McKean–Vlasov setting.

In the setting of Theorem 5, we take $N\asymp\kappa d/\varepsilon^{2}$. We use
the same sampling guarantees under (LSI) as in the prior settings.

We also note that in order to apply the log-concave sampling guarantees, we
must check that $\mu^{1:N}$ is log-smooth. This follows from Assumption 6.
Indeed,

$\displaystyle\lVert\nabla\log\mu^{1:N}(x^{1:N})-\nabla\log\mu^{1:N}(y^{1:N})\rVert$
$\displaystyle=\frac{2}{\sigma^{2}}\sqrt{\sum_{i=1}^{N}{\lVert\nabla_{\mathcal{W}_{2}}\mathcal{F}(\rho_{x^{1:N}},x^{i})-\nabla_{\mathcal{W}_{2}}\mathcal{F}(\rho_{y^{1:N}},y^{i})\rVert^{2}}}$
$\displaystyle\leq\frac{2\sqrt{2}\,\beta}{\sigma^{2}}\sqrt{\sum_{i=1}^{N}\bigl{(}\lVert
x^{i}-y^{i}\rVert^{2}+\mathcal{W}_{1}^{2}(\rho_{x^{1:N}},\rho_{y^{1:N}})\bigr{)}}$
$\displaystyle\leq\frac{4\beta}{\sigma^{2}}\,\lVert x^{1:N}-y^{1:N}\rVert\,.$