<figure><img src="image_1.png" /></figure>

HAL
open science

# Self-interacting approximation to McKean-Vlasov long-time limit: a Markov chain Monte Carlo method

Kai Du, Zhenjie Ren, Florin Suciu, Songbo Wang

► To cite this version:

Kai Du, Zhenjie Ren, Florin Suciu, Songbo Wang. Self-interacting approximation to McKean-Vlasov long-time limit: a Markov chain Monte Carlo method. Journal de Mathématiques Pures et Appliquées, In press, 205, pp.103782. 10.1016/j.matpur.2025.103782. hal-05309416

HAL Id: hal-05309416

https://hal.science/hal-05309416v1

Submitted on 11 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Self-interacting approximation
to McKean–Vlasov long-time limit:
a Markov chain Monte Carlo method

Kai Du¹, Zhenjie Ren², Florin Suciu³ and Songbo Wang⁴

¹Shanghai Center for Mathematical Sciences, Fudan University, Shanghai, China

²LaMME, Université Évry Paris-Saclay, Évry, France

³CEREMADE, Université Paris Dauphine-PSL, Paris, France

⁴CMAP, École polytechnique, IP Paris, Palaiseau, France

October 11, 2025

Abstract

For a certain class of McKean–Vlasov processes, we introduce proxy pro-
cesses that substitute the mean-field interaction with self-interaction, em-
ploying a weighted occupation measure. Our study encompasses two key
achievements. First, we demonstrate the ergodicity of the self-interacting dy-
namics, under broad conditions, by applying the reflection coupling method.
Second, in scenarios where the drifts are negative intrinsic gradients of convex
mean-field potential functionals, we use entropy and functional inequalities to
demonstrate that the stationary measures of the self-interacting processes ap-
proximate the invariant measures of the corresponding McKean–Vlasov pro-
cesses. As an application, we show how to learn the optimal weights of a
two-layer neural network by training a single neuron.

Résumé

Pour une certaine classe de processus de McKean-Vlasov, nous introdui-
sons des processus proxy qui remplacent l'interaction de champ moyen par
une auto-interaction, en utilisant une mesure d'occupation pondérée. Notre
étude comprend deux contributions principales. Premièrement, nous démon-
trons l'ergodicité de la dynamique auto-interactive, sous des conditions gé-
nérales, en appliquant la méthode de couplage par réflexion. Deuxièmement,
dans les cas où les dérives sont des gradients intrinsèques négatifs de fonc-
tionnelles de potentiel de champ moyen convexes, nous utilisons l'entropie et
des inégalités fonctionnelles pour montrer que les mesures stationnaires des
processus auto-interactifs approchent les mesures invariantes des processus
de McKean-Vlasov correspondants. Comme application, nous montrons com-
ment apprendre les poids optimaux d'un réseau de neurones à deux couches
en entraînant un seul neurone.

2020 Mathematics Subject Classification. Primary 60K35; Secondary 62L20, 60F05.

Keywords. McKean–Vlasov process, self-interaction, gradient flow, long-time behavior.
# 1 Introduction

In this paper we develop a novel method to approximate the invariant measure of the non-degenerate McKean-Vlasov dynamics

$$dX_t = b(\text{Law}(X_t), X_t) dt + dB_t, \quad (1)$$

where B is a standard Brownian motion in $\mathbb{R}^d$. The McKean-Vlasov dynamics characterize the mean field limit of interacting particles, and they have widespread applications, encompassing fields such as granular materials [7, 10, 11], mathematical biology [31], statistical mechanics [36], and synchronization of oscillators [33]. More recently, there has been a growing interest in the role of such dynamics in the context of mean field optimization for training neural networks [37, 16, 30, 39, 38, 15, 17].

In order to simulate the invariant measure of (1), we turn to the corresponding N-particle approximation, i.e., the dynamics

$$dX_t^i = b\left(\frac{1}{N}\sum_{i=1}^{N} \delta_{X_t^i}, X_t^i\right) dt + dB_t^i, \quad \text{for } i = 1, \dots, N, \quad (2)$$

where $(B^i)_{1\le i\le N}$ are independent Brownian motions. It is expected that the empirical measure $\frac{1}{N}\sum_{i=1}^N \delta_{X_t^i}$ of the N-particle system can approximate the McKean-Vlasov long-time limit when N and t are both large enough. For fixed N, to ensure control over the distance between the McKean-Vlasov dynamics in (1) and the N-particle system in (2), throughout the entire time horizon, the literature has proposed various uniform-in-time propagation of chaos results under different scenarios, see for example, [28, 23, 19, 34, 14].

When the drift b does not depend on the marginal distribution, $\text{Law}(X_t)$, the diffusion process is Markovian. Under general conditions, we can leverage Birkhoff's ergodic theorem [8] to approximate its invariant measure using the occupation measure

$$\overline{m}_t := \frac{1}{t} \int_0^t \delta_{X_s} \, ds,$$

as $t \to \infty$. In scenarios where the drift takes the form of gradient $b(x) = -\nabla U(x)$, this ergodic property of the Markov diffusion lays the groundwork for various Markov chain Monte Carlo methods, such as Metropolis adjusted Langevin algorithm [42, 41] and unadjusted Langevin algorithm [24]. Motivated by the Markovian ergodicity, the recent paper [21] studied the following self-interacting process:

$$dX_t = b(\overline{m}_t, X_t) dt + dB_t, \quad (3)$$

where the dependence on the marginal distribution in McKean-Vlasov diffusion (1) is replaced by the occupation measure, that is, the empirical mean of the trajectory $(X_s)_{s\in[0,t]}$. In [21], the authors successfully demonstrated that, in the regime of weak interaction, where the dependence on the marginal distribution is sufficiently small, the occupation measures $(\overline{m}_t)_{t\ge 0}$ of the self-interacting process (3) also converge towards the invariant measure of (1) as $t \to \infty$. Remarkably, from a practical point of view, simulating the occupation measure of the self-interacting process (3) only requires a single particle, which distinguishes it from the conventional N-particle approximation (2). It is worth noting that the investigation into
the long-time behavior of the self-interacting diffusions can be traced back to the pioneering works of Cranston and Le Jan [18] and Raimond [40].

Building upon the discovery in [21], this paper ventures into uncharted territory, where the mean field interaction need not to be inherently weak. We propose to study the self-interacting particle with exponentially decaying dependence on its trajectory:

$$ 
\begin{aligned}
dX_t &= b(m_t, X_t) dt + dB_t, \\
dm_t &= \lambda(\delta_{X_t} - m_t) dt,
\end{aligned} 
\quad (4) $$

where $\lambda$ is a positive constant. Integrating the second equation of (4), we find

$$ m_t = e^{-\lambda t} m_0 + \int_0^t \lambda e^{-\lambda(t-s)} \delta_{X_s} ds, $$

that is to say, the measure $m_t$ is an exponentially weighted occupation measure with emphasis on the recent past. The dynamics (4) is a time-homogeneous Markov process and we show its exponential ergodicity in the first part of the paper. Although the state space where the random variable $(X_t, m_t)$ lives is infinite-dimensional, we have a non-degenerate noise in the $X$ component and an almost sure contraction in the $m$ component, which render such ergodicity possible. Under suitable conditions for the drift $b$, we show in Theorem 1 by a reflection coupling approach that the Markov process is contractive for a Wasserstein distance. This implies that the stationary measure $\rho^\lambda$ for the Markov process exists, and is unique and globally attractive. Notably, the conditions that we impose on $b$ do not imply the uniqueness of invariant measure for the McKean–Vlasov (1), and cover the case of the ferromagnetic Curie–Weiss model in Section 4. Here, we also remark that the exact weight for the measure $m_t$ is not important and we work with a time-homogeneous Markov structure only for convenience. We could, for example, take the alternative weighting

$$ \frac{\lambda}{1 - e^{-\lambda t}} \int_0^t e^{-\lambda(t-s)} \delta_{X_s} ds $$

to remove the dependency on the initial value.

We then proceed to investigate properties of the stationary $\rho^\lambda$ in Theorem 3. A crucial feature of the self-interacting process (4) is that when $\lambda \to 0$, the dynamics of the measure $m_t$ becomes slow, while the rate of the $X_t$ dynamics remains roughly unchanged. Suppose additionally that for some mean field energy functional $F: \mathcal{P}(\mathbb{R}^d) \to \mathbb{R}$, the drift is its negative intrinsic gradient:

$$ b(m, x) = -D_m F(m, x) = -\nabla_x \frac{\delta F}{\delta m}(m, x). $$

The double time-scale structure allows us to speculate that $\text{Law}(X_t)$ rapidly relaxes to the local-in-time equilibrium

$$ \hat{m}_t := \frac{\exp(-2 \frac{\delta F}{\delta m}(m_t, x))}{\int \exp(-2 \frac{\delta F}{\delta m}(m_t, x')) dx'} dx $$

so that the measure evolution is effectively described by

$$ \partial_t m_t = \lambda(\hat{m}_t - m_t). \quad (5) $$
This dynamics is called entropic fictitious play in a previous work of F. Chen and two of the authors [13] and this point of view plays a key role in various literatures, notably the series of works of Benaïm, Ledoux and Raimond [3, 4, 5, 6] and the article of Kleptsyn and Kurtzmann [32]. The main novelty of our method is that we provide a quantitative justification of the intuition presented above, and we are no longer restricted to the case of two-body interaction potentials. More precisely, letting $(X, m)$ be a random variable distributed as the stationary measure $\rho^\lambda$, we provide an entropy bound in Proposition 7 that measures in a way the distance between the conditional distribution Law$(X|m)$ and $\hat{m}$, relying crucially on the log-Sobolev inequality for $\hat{m}$. This method requires unfortunately a finite-dimensional dependency of the mean field in the energy functional $F$, which we explain in Remark 5 in detail. The entropy bound, together with an inherent gradient structure of the dynamics, is then used in the rest of the proof of Theorem 3 to show that the random measure $m$ solves approximately the stationary condition for the entropic fictitious play (5):

$$ \hat{m} = m. $$

In the case that the energy $F$ is convex, the equation above has a unique solution $m_*$, which is also the invariant measure for the McKean–Vlasov dynamics (1), and we show that the stationary measure $\rho^\lambda$ is in fact close to $m_* \otimes \delta_{m_*}$ for small $\lambda$.

The self-interacting dynamics (4) can be also thought as an intermediate scheme between the entropic fictitious play (5), which corresponds to the limit $\lambda \to 0$, and the linear dynamics

$$ dX_t = b(\delta_{X_t}, X_t) dt + dB_t, $$

which corresponds to the limit $\lambda \to \infty$. From a computational point of view, the linear dynamics is easy to sample and relaxes rapidly, but in the long time does not yield the McKean–Vlasov’s long-time limit. The entropic fictitious play reaches high precision in the long time, but at each time step, it requires usually a costly Monte Carlo run to sample the Gibbs measure $\hat{m}_t$. The self-interacting dynamics lies exactly in between by mixing the two time scales.

As a final note on the terminology, although the words “stationary” and “invariant” have almost identical meanings in the context of stochastic process, we always say “invariant measure” when referring to the McKean–Vlasov process (1), and “stationary measure” when referring to the self-interacting process (4). We hope this artificial distinction would reduce possible confusions for the readers.

The rest of the paper is organized as follows. The main results are stated in Section 2. Before moving to their proofs, we apply our results to the training of a two-layer neural network in Section 3 and to a Curie–Weiss model for ferromagnets in Section 4. Ergodicity of the self-interacting dynamics, i.e., Theorem 1, is proved in Section 5. In Section 6, we prove Theorem 3, which characterizes the stationary measure of the self-interacting process. Finally, a technical result and its proof, and the numerical algorithm are included in the appendices.

## 2 Main results

We state and discuss our main results in this section. First, we study the contrac- tivity of the self-interacting process (4), and in particular, the exponential conver- gence to its unique stationary measure is obtained. Then, focusing on the gradient
case, we quantify the distance between the self-interacting stationary measure and the corresponding McKean–Vlasov invariant measure. Finally, applying both the results, we propose an annealing scheme so that the self-interacting dynamics converges to the McKean–Vlasov invariant measure.

To avoid extra assumptions on the drift functional b, we will always assume the existence and the uniqueness of strong solution to the self-interacting process (4) without explicitly mentioning it in the rest of the paper.

**Assumption.** Given any filtered probability space supporting a Brownian motion $(B_t)_{t \ge 0}$ and satisfying the usual conditions, for any initial conditions $(X_0, m_0)$ measurable to the initial $\sigma$-algebra and taking value in $\mathbb{R}^d \times \mathcal{P}(\mathbb{R}^d)$, there exists a unique adapted stochastic process $(X_t, m_t)_{t \ge 0}$ such that for all $t \ge 0$,

$$X_t = \int_0^t b(m_s, X_s) \, ds + B_t + X_0,$$

$$m_t = \lambda \int_0^t (\delta_{X_s} - m_s) \, ds + m_0.$$

One may easily find various sufficient conditions for the assumption above. For example, if $b: \mathcal{P}(\mathbb{R}^d) \times \mathbb{R}^d \to \mathbb{R}^d$ is $W_1$-Lipschitz continuous in measure and Lipschitz continuous in space, then by Cauchy–Lipschitz arguments, we know that the assumption is satisfied.

## 2.1 Contractivity of the self-interacting diffusion

We first present the results on the contractivity of the self-interacting dynamics (4), from which follows the convergence to its unique stationary measure. We start with two basic definitions. First, we define moduli of continuity.

**Definition 1** (Modulus of continuity). We say that $\omega: [0, \infty) \to [0, \infty)$ is a modulus of continuity if it satisfies the following properties:

* $\omega(0) = 0$;

* $\omega$ is continuous and non-decreasing;

* for all $h, h' \ge 0$, we have $\omega(h + h') \le \omega(h) + \omega(h')$.

Note that a modulus of continuity necessarily has at most linear growth according to our definition. We also introduce the notion of semi-monotonicity following Eberle [25].

**Definition 2** (Semi-monotonicity). We say that $\kappa: (0, \infty) \to \mathbb{R}$ is a semi-monotonicity function for a vector field $v: \mathbb{R}^d \to \mathbb{R}^d$ if

$$(v(x) - v(x')) \cdot (x - x') \le -\kappa(|x - x'|)|x - x'|^2$$

holds for every $x, x' \in \mathbb{R}^d$ with $x' \ne x$. We say $\kappa$ is a uniform semi-monotonicity function of a family of vector fields if it is a semi-monotonicity function of each member.

In this subsection, we impose the following assumption on the drift of the McKean–Vlasov dynamics (1).
**Assumption 1.** The drift b satisfies the following conditions:

1. For any fixed $m \in \mathcal{P}(\mathbb{R}^d)$, the vector field $x \mapsto b(m, x)$ is uniformly equicontinuous and has a uniform semi-monotonicity function $\kappa_b$, given by $\kappa_b(x) = \kappa_0 - M_b/x$ for some $\kappa_0 > 0$ and $M_b \ge 0$.

2. There exist a bounded modulus of continuity $\omega : [0, \infty) \to [0, M_\omega]$ and a constant $L \ge 0$ such that

$$|b(m, x) - b(m', x)| \le L W_\omega(m, m')$$

for every $m, m' \in \mathcal{P}(\mathbb{R}^d)$ and every $x \in \mathbb{R}^d$. Here $W_\omega$ is the Wasserstein distance

$$W_\omega(m, m') = \inf_{\pi \in \Pi(m, m')} \int \omega(|x - x'|)\pi(dx \, dx').$$

Using reflection coupling, we derive the following result.

**Theorem 1.** Suppose Assumption 1 hold. Let $(X_t, m_t)_{t \ge 0}$, $(X'_t, m'_t)_{t \ge 0}$ be two processes following the dynamics (4) for some $\lambda > 0$ such that the first marginals of their initial values $X_0, X'_0$ have finite first moments. Define the following metric on $\mathbb{R}^d \times \mathcal{P}(\mathbb{R}^d)$:

$$d_{\lambda}((x, m), (x', m')) = |x - x'| + \frac{2L}{\lambda} W_{\omega}(m, m')$$

and denote the corresponding Wasserstein distance on $\mathcal{P}_1(\mathbb{R}^d \times \mathcal{P}(\mathbb{R}^d))$ by $W_{d_\lambda}$. Then, we have

$$W_{d_\lambda}(\text{Law}(X_t, m_t), \text{Law}(X'_t, m'_t)) \le Ce^{-ct}W_{d_\lambda}(\text{Law}(X_0, m_0), \text{Law}(X'_0, m'_0)),$$

where the constants C, c are given by

$$C = 1 + \frac{2M}{\sqrt{K_0}} \exp\left(\frac{M^2}{4K_0}\right),$$

$$c = \left( \frac{1}{K_0} + \frac{2M}{K_0^{3/2}} \exp\left(\frac{M^2}{4K_0}\right) \right)^{-1}$$

for $M = M_b + 2LM_\omega$ and $K_0 = \min(\kappa_0, \frac{\lambda}{2})$.

The proof of Theorem 1 is postponed to Section 5.

**Remark 1** (On the assumption). The first condition on the semi-monotonicity of the vector field $x \mapsto b(m, x)$ is stronger than those used in [25], in that we require a more gentle singularity in $\kappa(x)$ for $x$ close to 0. This is because in this work, we are concerned with a good convergence rate when the parameter $\lambda \to 0$ (see the following remark for more discussions) and it will become clear in the proof that this stronger requirement on the semi-monotonicity is necessary for our purpose. Nevertheless, this condition is not too difficult to fulfill. Upon defining $b_0(x) = -\kappa_0 x$ and $b_1(m, x) = b(m, x) - b_0(x)$, the first condition of Assumption 1 is equivalent to

$$(b_1(m, x) - b_1(m, x'), x - x') \le M_b,$$

and this holds true when $\sup_{(m,x) \in \mathcal{P}(\mathbb{R}^d) \times \mathbb{R}^d} |b_1(m, x)| \le M_b/2$ in particular.
Remark 2 (Rate of convergence). We discuss the rate of convergence across three parameter ranges.

In the first scenario, when the drift $b$ is $\kappa_0$-strongly monotone for some $\kappa_0 > 0$, i.e., $M_b = 0$, and when there is no mean field interaction ($L = 0$), we have $M = 0$ and $K_0 = \min(\kappa_0, \frac{\lambda}{2})$. Consequently, $C = 1$ and $c = \min(\kappa_0, \frac{\lambda}{2})$. It is worth noting that under these conditions, the component $X$ exhibits exponential contraction with a rate of $\kappa_0$, while $m$ contracts at a rate no greater than $\lambda$. In this case, the best contraction rate for the joint process is $\min(\kappa_0, \lambda)$. Thus our method yields a contraction rate that remains at least half of the optimal one.

In the second scenario, when $\lambda$ is small but non-zero (with self-interaction), we obtain $c \sim 2M\lambda^{3/2} \exp(-M^2/2\lambda)$ and $C \sim 2M\lambda^{-1/2} \exp(-M^2/2\lambda)$. We note that such exponentially slow convergence also arises in the kinetic Langevin process as the damping parameter approaches zero; see Eberle, Guillin and Zimmer [26, Section 2.6] for further discussion.

Finally, for $\lambda > 2\kappa_0$, the contractivity constants $C$, $c$ become independent of $\lambda$, consistent with the intuition that the self-interacting diffusion becomes linear in the large $\lambda$ limit.

Now we discuss a few examples satisfying Assumption 1.

*Example 1 (Two-body interaction).* Consider $b(m, x) = b_0(x) + b_1(m, x)$, where

$$b_1(m, x) = \int K(x, x')m(dx').$$

Suppose furthermore that

$$\sup_{x \in \mathbb{R}^d} \|K(x, \cdot)\|_{W^{1,\infty}} = \sup_{x \in \mathbb{R}^d} \max(\|K(x, \cdot)\|_{L^\infty}, \|\nabla K(x, \cdot)\|_{L^\infty}) \le M,$$

that is to say, the mapping $y \mapsto K(x, y)$ is $M$-bounded and $M$-Lipschitz continuous for every $x$. Thus, we have

$$|b(m, x) - b(m', x)| \le \left| \int K(x, x')(m - m')(dx') \right| \le M W_{\omega}(m, m')$$

for the modulus of continuity $\omega(x) = \min(x, 2)$. Therefore Assumption 1 is satisfied once $b_0$ is uniformly Lipschitz and has a semi-monotonicity function $\kappa_0(x) = \kappa_0 - M_1/x$ for some $\kappa_0 > 0$ and $M_1 \ge 0$.

We can generalize the example above to drifts that depend on the measure in a non-linear way.

*Example 2 ($C^1$ functional).* Suppose $m \mapsto b(m, x)$ is $C^1$ differentiable in the sense that there exists a continuous and bounded mapping $\frac{\delta b}{\delta m} : \mathcal{P}(\mathbb{R}^d) \times \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ such that

$$b(m, x) - b(m', x) = \int_0^1 \int \frac{\delta b}{\delta m}((1-t)m + tm', x, x')(m-m')(dx') dt$$

for all $m, m' \in \mathcal{P}(\mathbb{R}^d)$ and $x \in \mathbb{R}^d$. If $\sup_{m,x} \|\frac{\delta b}{\delta m}(m, x, \cdot)\|_{W^{1,\infty}}$ is finite and the vector fields $x \mapsto b(m, x)$ are uniformly Lipschitz and share a uniform semi-monotonicity function $\kappa_0(x) = \kappa_0 - M_1/x$ for some $\kappa_0 > 0$ and $M_1 \ge 0$, then by the same argument as in Example 1, Assumption 1 is satisfied.

We now examine the stationary measure of the self-interacting process (4).
**Definition 3.** We call a probability measure $P \in \mathcal{P}(\mathbb{R}^d \times \mathcal{P}(\mathbb{R}^d))$ stationary to the self-interacting diffusion (4) if the stochastic process $(X_t, m_t)_{t \ge 0}$ with initial value $\text{Law}(X_0, m_0) = P$ satisfies $\text{Law}(X_t, m_t) = P$ for all $t \ge 0$.

The definition above makes sense since we have assumed the existence and uniqueness of strong solution.

By the Banach fixed point theorem in the metric space $\mathcal{P}_1(\mathbb{R}^d \times \mathcal{P}(\mathbb{R}^d))$ and standard arguments, Theorem 1 implies the existence and uniqueness of stationary measure of the self-interacting process (4).

**Corollary 2.** Under Assumption 1, for every $\lambda > 0$, there exists a unique stationary measure of the self-interacting diffusion (4) in $\mathcal{P}(\mathbb{R}^d \times \mathcal{P}(\mathbb{R}^d))$ whose first marginal distribution on $\mathbb{R}^d$ has finite first moment.

Finally, we note that although Theorem 1, along with Corollary 2, implies that the self-interacting process (4) converges to its stationary measure exponentially, its assumptions are not sufficient to establish the uniqueness of invariant measure of the McKean–Vlasov process (1), as illustrated by the Curie–Weiss example in Section 4. So in general, there is no hope that the self-interacting stationary measure approximates the McKean–Vlasov invariant measure.

## 2.2 Stationary measure in the gradient case

In this subsection, we study the stationary measure of the self-interacting dynamics (4), provided that the drift $b$ is the negative intrinsic gradient of a finite-dimensional mean field functional, whose precise meaning will be explained in the following. In particular, we aim at proving that, in this case, the stationary measure of (4) provides an approximation to the invariant measure of the McKean–Vlasov dynamics (1). We fix a positive $\lambda$ in this subsection.

The first assumption in the subsection is that the drift $b$ corresponds to a gradient descent whose dependency in the mean field is finite-dimensional.

**Assumption 2 (Finite-dimensional mean field).** For a closed convex set $\mathcal{K} \subset \mathbb{R}^D$, there exists a function

$$ \ell = (\ell^1, \dots, \ell^D) \in \mathcal{C}^1(\mathbb{R}^d; \mathcal{K}) $$

whose gradient is of at most linear growth, and a function $\Phi \in \mathcal{C}^2(\mathbb{R}^D; \mathbb{R})$ whose Hessian $\nabla^2\Phi$ is bounded, such that the drift term $b$ reads

$$ \begin{aligned} b(m, x) &= -\nabla\Phi(\langle\ell, m\rangle) \cdot \nabla\ell(x) = -\nabla\Phi\left(\int \ell(x)m(dx)\right) \cdot \nabla\ell(x) \\ &= -\sum_{\nu=1}^{D} \nabla_{\nu}\Phi\left(\int \ell(x)m(dx)\right) \nabla\ell^{\nu}(x). \end{aligned} $$

In other words, for the mean field functional $F: \mathcal{P}_2(\mathbb{R}^d) \to \mathbb{R}$ defined by

$$ F(m) = \Phi(\langle\ell, m\rangle) = \Phi\left(\int \ell(x')m(dx')\right), $$

the drift $b$ is the negative intrinsic derivative:

$$ b(m, x) = -D_m F(m, x) = -\nabla_x \frac{\delta F}{\delta m}(m, x). $$
We shall also impose the following conditions on a family of probability measures related to $\Phi$ and $\ell$.

**Assumption 3 (Uniform LSI).** The probability measures $(\hat{m}_y)_{y \in \mathcal{K}}$ on $\mathbb{R}^d$ determined by

$$ \hat{m}_y(dx) \propto \exp(-2\nabla\Phi(y) \cdot \ell(x)) dx $$

are well defined and satisfy a uniform $C_{LS}$-logarithmic Sobolev inequality for some $C_{LS} \ge 0$. That is to say, for all regular enough probability measure $m \in \mathcal{P}(\mathbb{R}^d)$ and all $y \in \mathcal{K}$, we have

$$ \int \log \frac{dm}{d\hat{m}_y} dm := H(m|\hat{m}_y) \le \frac{C_{LS}}{4} I(m|\hat{m}_y) := \frac{C_{LS}}{4} \int \left| \nabla \log \frac{dm}{d\hat{m}_y} \right|^2 \hat{m}_y, $$

where $dm/d\hat{m}_y$ is the Radon–Nikodým derivative between measures.

**Remark 3.** As mentioned in the introduction, the uniform log-Sobolev inequality is crucial to our method as it is used to obtain the entropy estimate in Proposition 7. This condition is often perceived as a strong one, but still it can be verified if for example

$$ 2\nabla\Phi(y) \cdot \ell(x) = U(x) + G(y, x) $$

for some strongly convex $U: \mathbb{R}^d \to \mathbb{R}$ and some bounded $G: \mathcal{K} \times \mathbb{R} \to \mathbb{R}$. Indeed, in this case $\hat{m}_y$ is a uniformly bounded perturbation of the strongly log-concave measure $e^{-U(x)} dx / \int e^{-U(x')} dx'$, so it satisfies a uniform log-Sobolev inequality by the Bakry–Émery condition [2] and the Holley–Stroock perturbation [29]. We note that this condition has also been exploited recently to obtain long-time behaviors of mean field Langevin [15, 38] and its particle systems [14, 12].

Finally, we assume that the following quantitative bound on $\Phi$ and $\ell$.

**Assumption 4 (Bound).** The following quantity is finite:

$$ \begin{aligned} M_2 &:= \sup_{x \in \mathbb{R}^d, y \in \mathbb{R}^D} |\nabla^2 \Phi(y)^{1/2} \nabla \ell(x)|^2 \\ &= \sup_{x \in \mathbb{R}^d, y \in \mathbb{R}^D} \sup_{\substack{a \in \mathbb{R}^d \\ |a|=1}} a^\top \nabla \ell(x)^\top \nabla^2 \Phi(y) \nabla \ell(x) a. \end{aligned} $$

**Remark 4.** Under the three assumptions above, if $\Phi$ is additionally convex, then there exists a unique invariant measure $m_*$ of the McKean–Vlasov dynamics (1) by [14, Proposition 4.4 and Corollary 4.8], and the convergence to the invariant measure is exponentially fast [14, Theorem 2.1]. In fact, the convexity of $\Phi$ implies precisely the functional convexity of the mean field energy $F$ considered in [14].

The main discovery of this paper characterizes the distance between the stationary measure $P$ of the self-interacting dynamics (4) and $m_* \otimes \delta_{m_*}$.

**Theorem 3.** Let Assumptions 2, 3 and 4 hold true. Suppose that $P = P^\lambda \in \mathcal{P}_4(\mathbb{R}^d \times \mathcal{P}_4(\mathbb{R}^d))$ is a stationary measure of the self-interacting process (4) in the sense of Definition 3 that has finite fourth moment:

$$ \iint (|x|^4 + |x'|^4 m(dx')) P(dx \, dm) < \infty. $$

Let $(X, m)$ be a random variable distributed as $P$. Denote by $\rho^1$ and $\rho^2$ the probability distributions of the random variables $X$ and $\langle \ell, m \rangle = \int \ell(x)m(dx)$ respectively.
1. Suppose $\Phi$ is convex. Denote in the case by $m_*$ the unique invariant measure of the McKean-Vlasov dynamics (1). Define

$$
\begin{align*}
y_* &:= \langle \ell, m_* \rangle = \int \ell(x) m_*(dx), \\
\lambda_0 &:= \frac{1}{48 M_2 C_{LS}^2 (1 + 2M_2 C_{LS} (M_2^2 C_{LS}^2 + 1)^{1/2})}, \\
H &:= \frac{C_{LS}(D + 24M_2 C_{LS} d)\lambda}{2 - 96 M_2 C_{LS}^2 (1 + 2M_2 C_{LS} (M_2^2 C_{LS}^2 + 1)^{1/2})\lambda}.
\end{align*}
$$

Then, for $\lambda \in (0, \lambda_0)$, we have

$$
\begin{align*}
v(\rho^2) &:= \mathbb{E} \left[ \iint \int_0^1 \frac{\delta^2 F}{\delta m^2} \left( (1-t)m + tm_*, x', x'' \right) \mathrm{d}t (m-m_*)^{\otimes 2} (dx' \, dx'') \right] \\
&= \iint_0^1 (y-y_*)^{\top} \nabla^2 \Phi \left( (1-t)y + ty_* \right) (y-y_*) \mathrm{d}t \rho^2(dy) \\
&\le 4M_2 C_{LS} (M_2^2 C_{LS}^2 + 1) H,
\end{align*}
$$

and

$$
\begin{align*}
W_2^2(\rho^1, m_*) &\le (2C_{LS} + 4M_2 C_{LS}^2 (M_2^2 C_{LS}^2 + 1)^{1/2})H, \\
\|\rho^1 - m_*\|_{\mathrm{TV}}^2 &\le (4 + 8M_2 C_{LS} (M_2^2 C_{LS}^2 + 1)^{1/2})H.
\end{align*}
$$

2. If in addition to the convexity of $\Phi$,

$$
M_1 := \sup_{x \in \mathbb{R}^d, y \in \mathbb{R}^D} \ell(x)^{\top} \nabla^2 \Phi(y) \ell(x) < \infty,
$$

then the three inequalities above hold also for all $\lambda \in (0, \infty)$, with $H$ replaced
by

$$
H' := \frac{C_{LS}}{2} (D + 2M_1)\lambda.
$$

3. If $\Phi$ is concave, then for $\hat{y} := \langle \ell, \hat{m}_y \rangle$, we have

$$
- \int (\hat{y} - y)^{\top} \nabla^2 \Phi(y) (\hat{y} - y) \rho^2(dy) \le \frac{M_2 C_{LS}^2 D}{2} \lambda.
$$

The proof of Theorem 3 is postponed to Section 6.

Remark 5 (Dependence on the dimension *D*). Readers may have observed that, in our framework, the value of the functional *F*(*m*) may only depend on the *D*-dimensional vector *∫* ℓ(*x*)*m*(*dx*), and this corresponds to “cylindrical functions” considered in [1, Definition 5.1.11]. Given a continuous functional *F* on *P*(**R*<sup>*d*</sup>), for every compact subset *S* ⊂ **R*<sup>*d*</sup>, we can construct, according to the Stone-Weierstrass theorem, a sequence of functions ℓ<sub>*n*</sub>: **R*<sup>*d*</sup> → **R**<sup>*D*</sup> and Φ<sub>*n*</sub>: **R**<sup>*D*</sup> → **R** such that the cylindrical functionals *F*<sub>*n*</sub>(*m*) = Φ<sub>*n*</sub>(∫ ℓ<sub>*n*</sub> *dm*) approach *F* in the uniform topology of *C*(*P*(*S*)) (see [27, Lemma 2] for example). However, the dimension *D*<sub>*n*</sub> may tend to infinity when *n* → ∞. Since all the upper bounds in Theorem 3 depend on the dimension *D* linearly, our analysis and findings cannot be directly applied to more general functionals on *P*(**R*<sup>*d*</sup>).
Remark 6 (Meaning of $M_1$ and $M_2$). The second-order functional derivative of $F$ reads

$$ \frac{\delta^2 F}{\delta m^2}(m, x', x'') = \ell(x')^\top \nabla^2 \Phi(y) \ell(x''), $$

and in the case of convex $\Phi$, satisfies the Cauchy–Schwarz inequality:

$$ \begin{aligned} \left| \frac{\delta^2 F}{\delta m^2}(m, x', x'') \right| &= |\ell(x')^\top \nabla^2 \Phi(y) \ell(x'')| \\ &\le |\ell(x')^\top \nabla^2 \Phi(y) \ell(x')|^{1/2} |\ell(x'')^\top \nabla^2 \Phi(y) \ell(x'')|^{1/2} \\ &= \left| \frac{\delta^2 F}{\delta m^2}(m, x', x') \right|^{1/2} \left| \frac{\delta^2 F}{\delta m^2}(m, x'', x'') \right|^{1/2} \le M_1. \end{aligned} $$

Similarly, the second-order intrinsic derivative satisfies

$$ \begin{aligned} |D_m^2 F(m, x', x'')| &= |\nabla \ell(x')^\top \nabla^2 \Phi(y) \nabla \ell(x'')| \\ &\le |\nabla \ell(x')^\top \nabla^2 \Phi(y) \nabla \ell(x')|^{1/2} |\nabla \ell(x'')^\top \nabla^2 \Phi(y) \nabla \ell(x'')|^{1/2} \\ &= |D_m^2 F(m, x', x')|^{1/2} |D_m^2 F(m, x'', x'')|^{1/2} \le M_2. \end{aligned} $$

Moreover, by taking $x' = x''$ in the inequalities above, we find that $M_1$ and $M_2$ are the respective suprema of the second-order flat and intrinsic derivatives of the functional $F$.

We illustrate in the following example that the order of $\lambda$ when $\lambda \to 0$ for the variance of $\langle \ell, m \rangle$ under $P$ in Theorem 3 (i.e., the first claim) is optimal.

Example 3 (Optimality of the order of $\lambda$). Consider the mean field functional $F: P_2(\mathbb{R}^d) \to \mathbb{R}$ given by

$$ F(m) = \frac{1}{2} \int |x|^2 m(dx) + \frac{1}{2} \left| \int x m(dx) \right|^2. $$

By taking $D = d+1$, $\ell(x) = (x, |x|^2/2)^\top$, $\Phi(y_0, y_1) = y_0^2/2 + y_1$, the mean field functional $F$ is covered by the cylindrical setting (namely Assumption 2) of Theorem 3. Moreover, the function $\Phi$ is convex.

The corresponding gradient dynamics (1) is then characterized by the drift

$$ b(m, x) = -x - \int x' m(dx'), $$

and its unique invariant measure $m_*$ is $N(0, 1/2)$. We explicitly compute quantities related to the stationary measure $P^\lambda$ of the self-interacting dynamics (4) in the following. The dynamics reads

$$ \begin{aligned} dX_t &= -X_t dt - \int x' m_t(dx') dt + dB_t, \\ dm_t &= \lambda(\delta_{X_t} - m_t) dt, \end{aligned} $$

and has a unique strong solution by Cauchy–Lipschitz arguments. Upon defining $Y_{0,t} = \int x' m_t(dx')$, the process has the finite-dimensional projection:

$$ \begin{aligned} dX_t &= (-X_t - Y_{0,t}) dt + dB_t, \\ dY_{0,t} &= \lambda(-Y_{0,t} + X_t) dt. \end{aligned} $$
The finite-dimensional dynamics has a unique invariant measure that is a centered Gaussian with the following covariance structure:

$$ 
\begin{aligned}
\mathrm{E}[X \otimes X] &= \frac{\lambda + 2}{4(\lambda + 1)} \mathbb{1}_{d \times d}, \\
\mathrm{E}[Y_0 \otimes Y_0] &= \mathrm{E}[X \otimes Y_0] = \frac{\lambda}{4(\lambda + 1)} \mathbb{1}_{d \times d}.
\end{aligned}
$$

Hence, the exact distances, or bounds thereof, read

$$ 
\begin{aligned}
\mathrm{E}[|Y_0|^2] &= \frac{d\lambda}{4(1 + \lambda)}, \\
W_2^2(\mathrm{Law}(X), m_*) &= \frac{d}{2} \left( 1 - \left( 1 - \frac{\lambda}{2(1 + \lambda)} \right)^{1/2} \right)^2, \\
\|\mathrm{Law}(X) - m_*\|_{\mathrm{TV}}^2 &\in \left[ \frac{1}{10000} \frac{d\lambda^2}{4(1 + \lambda)^2}, \frac{9}{4} \frac{d\lambda^2}{4(1 + \lambda)^2} \right],
\end{aligned}
$$

where the last mutual bound is due to [20, Theorem 1.1].

Now we try to verify the assumptions of Theorem 3. By the Kolmogorov extension theorem, we can construct a stationary Markov process $(X_t, Y_{0,t})_{t \in \mathbb{R}}$, defined on the whole real line, such that $\mathrm{Law}(X_t, Y_{0,t}) = \mathrm{Law}(X, Y_0)$ for all $t \in \mathbb{R}$. Then, by defining

$$ m_t = \lambda \int_{-\infty}^{t} e^{-\lambda(t-s)} \delta_{X_s} \, ds, $$

we recover the solution $(X_t, m_t)_{t \in \mathbb{R}}$ to the original infinite-dimensional dynamics. By construction, $\mathrm{Law}(X_t, m_t)$ is stationary and has finite fourth moments. The rest of the assumptions of Theorem 3 can be satisfied with the constants $C_{LS} = 1/2$ and $M_2 = 1$. Since $\Phi$ is convex, by the first claim of the theorem, we get

$$ 
\begin{aligned}
\mathrm{E}[|Y_0|^2] &\le \frac{5(13d + 1)\lambda}{8 - 48(2 + \sqrt{5})\lambda}, \\
W_2^2(\mathrm{Law}(X), m_*) &\le \frac{(2 + \sqrt{5})(13d + 1)\lambda}{8 - 48(2 + \sqrt{5})\lambda}, \\
\|\mathrm{Law}(X) - m_*\|_{\mathrm{TV}}^2 &\le \frac{(2 + \sqrt{5})(13d + 1)\lambda}{2 - 12(2 + \sqrt{5})\lambda},
\end{aligned}
$$

for $\lambda < 1/6(2 + \sqrt{5})$. So the results of Theorem 3 give the optimal order of $\lambda$ when $\lambda \to 0$ for the variance of the $Y$ variable, but possibly sub-optimal ones for the Wasserstein and total-variation distances in the $X$ direction.

## 2.3 A class of dynamics

In this subsection, we present a class of dynamics to which both Theorems 1 and 3 are applicable. This class encompasses in particular the neural network example that will be discussed in the following Section 3.

**Assumption 5.** The drift functional writes

$$ b(m, x) = -\nabla U(x) - \nabla \Phi_0 \left( \int \ell_0(x')m(dx') \right) \cdot \nabla \ell_0(x). \quad (6) $$
for some functions $U : \mathbb{R}^d \to \mathbb{R}$, $\Phi_0 : \mathbb{R}^D \to \mathbb{R}$, $\ell_0 : \mathbb{R}^d \to \mathbb{R}^D$ satisfying the following conditions:

* the function $U$ is $C^2$ continuous with bounded Hessian, i.e., $\|\nabla^2 U\|_{L^\infty} < \infty$, and its gradient admits a semi-monotonicity function $\kappa(x) = \kappa_0 - M/x$ for some $\kappa_0 > 0$ and $M \ge 0$.

* the probability measure $Z^{-1} \exp(-2U(x)) dx$, with $Z = \int \exp(-2U(x)) dx$, is well defined in $\mathcal{P}(\mathbb{R}^d)$, and satisfies a $C_{LS,0}$-logarithmic Sobolev inequality.

* the function $\Phi_0$ is $\kappa_{\Phi_0}$-strongly convex for some $\kappa_{\Phi_0} > 0$ and belongs to $C^2(\mathbb{R}^D) \cap W^{2,\infty}(\mathbb{R}^D)$.

* the function $\ell_0$ belongs to $C^1(\mathbb{R}^d; \mathbb{R}^D) \cap W^{1,\infty}(\mathbb{R}^d; \mathbb{R}^D)$.

**Proposition 4.** Under Assumption 5, there exists a unique stationary measure $P \in \mathcal{P}_1(\mathbb{R}^d \times \mathcal{P}(\mathbb{R}^d))$ to the self-interacting dynamics (4). Moreover, there exists $C > 0$, independent of $\lambda$, such that for $(X, m)$ distributed as $P$,

$$E[|\langle \ell_0, m - m_* \rangle|^2] + W_2^2(\text{Law}(X), m_*) + \|\text{Law}(X) - m_*\|_{TV}^2 \le C\lambda,$$

where $m_*$ is the unique invariant measure to the McKean-Vlasov process.

*Proof of Proposition 4.* We first verify the conditions of Theorem 1 to establish the existence and uniqueness of the stationary measure. As the drift $b$ has derivative

$$\frac{\delta b}{\delta m}(m, x, x') = -\nabla^2 \Phi_0 \left( \int \ell_0(x'')m(dx'') \right) \cdot \nabla \ell_0(x) \cdot \ell_0(x'),$$

we have

$$\left\| \frac{\delta b}{\delta m}(m, x, \cdot) \right\|_{W^{1,\infty}} \le \|\nabla^2 \Phi_0\|_{L^\infty} \|\nabla \ell_0\|_{L^\infty} \|\ell_0\|_{W^{1,\infty}}$$

for every $m \in \mathcal{P}(\mathbb{R}^d)$ and every $x \in \mathbb{R}^d$. Then the dynamics falls into the class considered in Example 2. The conditions of Theorem 1 are satisfied. Applying Corollary 2 gives the existence and the uniqueness of the stationary measure $P \in \mathcal{P}_1(\mathbb{R}^d \times \mathcal{P}(\mathbb{R}^d))$.

We proceed to verify the conditions of Theorem 3. Introduce the functions $\Phi : \mathbb{R}^{D+1} \to \mathbb{R}$, $\ell : \mathbb{R}^d \to \mathbb{R}^{D+1}$ defined respectively by

$$\begin{aligned} \Phi(y_0, y_1) &= \Phi_0(y_0) + y_1, && \text{for } y_0 \in \mathbb{R}^D \text{ and } y_1 \in \mathbb{R}, \\ \ell(x) &= (\ell_0(x), U(x)), && \text{for } x \in \mathbb{R}^d. \end{aligned}$$

Here the range set $\mathcal{K}$ of the mapping $\ell$ is taken as the whole space $\mathbb{R}^D$. In this way, the drift $b$ reads

$$b(m, x) = -\nabla \Phi \left( \int \ell(x')m(dx') \right) \cdot \nabla \ell(x)$$

so Assumption 2 is satisfied. Now we show the stationary measure $P$ of the dynamics (4) has finite fourth moment. Consider the functional

$$E(x, m) = \kappa_0^{-1} |x|^4 + \int |x'|^4 m(dx').$$
Along the self-interacting dynamics (4), we have, by Itô's formula,

$$ \frac{d}{dt} \mathrm{E}[E(X_t, m_t)] = 4\kappa_0^{-1} \mathrm{E}[b(m_t, X_t) \cdot X_t |X_t|^2] + (2d + 4)\kappa_0^{-1} \mathrm{E}[|X_t|^2] \\ - \lambda \mathrm{E} \left[ \int |x'|^4 m'(dx') \right] + \lambda \mathrm{E}[|X_t|^4]. $$

As the vector field $x \mapsto b(m, x)$ has weak monotonicity function $\kappa_b(x) = \kappa_0 - M_b/x$, we have

$$ b(m, x) \cdot x \le -\frac{\kappa_0}{2} |x|^2 + C_2 $$

for every $(m, x) \in \mathcal{P}(\mathbb{R}^d) \times \mathbb{R}^d$, for some $C_2 \ge 0$. The functional $E$ verifies the Lyapunov condition as well:

$$ \frac{d}{dt} \mathrm{E}[E(X_t, m_t)] \le -c_3 \mathrm{E}[E(X_t, m_t)] + C_3 $$

for some $c_3 > 0$, $C_3 \ge 0$. In consequence, the stationary measure $P$ of finite first moment must have finite fourth moment. Since

$$ \hat{m}_y(dx) \propto \exp(-2\nabla\Phi(y) \cdot \ell(x)) dx = \exp(-2\nabla\Phi_0(y) \cdot \ell_0(x)) \exp(-2U(x)) dx $$

with $\sup_{x,y} |\nabla\Phi_0(y) \cdot \ell_0(x)| \le \|\nabla\Phi_0\|_{L^\infty} \|\ell_0\|_{L^\infty}$, by the Holley-Stroock perturbation argument [29], we know that the measure $\hat{m}_y$ verifies a uniform $C_{LS}$-LSI with

$$ C_{LS} = C_{LS,0} \exp(4 \|\nabla\Phi_0\|_{L^\infty} \|\ell_0\|_{L^\infty}), $$

so Assumption 3 is satisfied. The constants $M_1, M_2$ in Theorem 3 and Assumption 4 are finite as

$$ \begin{aligned} |\ell(x)^\top \nabla^2 \Phi(y) \ell(x)| &= |\ell_0(x)^\top \nabla^2 \Phi_0(y) \ell_0(x)| \le \|\nabla^2 \Phi_0\|_{L^\infty} \|\ell_0\|_{L^\infty}^2, \\ |\nabla \ell(x) \nabla^2 \Phi(y) \nabla \ell(x)| &= |\nabla \ell_0(x)^\top \nabla^2 \Phi_0(y) \nabla \ell_0(x)| \le \|\nabla^2 \Phi_0\|_{L^\infty}^2 \|\nabla \ell_0\|_{L^\infty}. \end{aligned} $$

All the conditions of Theorem 3 are satisfied. Since for all $y \in \mathbb{R}^D$,

$$ \kappa_{\Phi_0} |\langle \ell_0, m - m_* \rangle|^2 \le \langle \ell_0, m - m_* \rangle^\top \nabla^2 \Phi_0(y) \langle \ell_0, m - m_* \rangle = \langle \ell, m - m_* \rangle^\top \nabla^2 \Phi(y) \langle \ell, m - m_* \rangle, $$

the claim of the proposition follows from the second case stated in the theorem. $\square$

Remark 7 (Open question). Applying further the convergence result of Theorems 1, we can obtain a bound on the difference between the marginal distribution of the non-stationary self-interacting diffusion (4) and the invariant measure of the McKean-Vlasov dynamics (1). Specifically, let $(X_t, m_t)$ denote the self-interacting process (4). Theorem 1 yields

$$ W_1(\mathrm{Law}(X_t, m_t), \mathrm{Law}(X, m)) \le C e^{-ct}, $$

where $C, c$ are the contractivity constants in the theorem. Let $\varphi : \mathbb{R}^K \to \mathbb{R}$ be a 1-Lipschitz test function. Combining this with Proposition 4, we can bound the following $L^1$ simulation error:

$$ \mathrm{E}[|\langle \varphi \circ \ell, m_t - m_* \rangle|] = \mathrm{E}[|\langle \varphi \circ \ell, m_t - m + m - m_* \rangle|] = O(e^{-ct} + \sqrt{\lambda}). $$
As noted in Remark 2, the contraction rate c depends on $1/\lambda$ exponentially, ren-
dering the above error bound impractical for applications.

This naturally raises the question of whether the method and results of Theo-
rem 3, which address the static case (i.e., the comparison between Law($X$, $m$) and
$m_* \otimes \delta_{m_*}$), can be extended to the dynamical setting (comparing Law($X_t$, $m_t$) and
$m_* \otimes \delta_{m_*}$). Unfortunately, we are currently unable to establish the entropy estimate
in Section 6.2 for the parabolic problem, so our approach does not yet apply in this
context. We leave this as an open problem for future research.

3 Numerical application to training two-layer neu-
ral networks

Let us recall the structure of two-layer neural networks and introduce the mean
field model for it. Consider an *activation function* $\varphi : \mathbb{R} \to \mathbb{R}$ satisfying

$\varphi$ is continuous and non-decreasing,

$$
\lim_{\theta \to -\infty} \varphi(\theta) = 0, \quad \lim_{\theta \to +\infty} \varphi(\theta) = 1.
$$

Define $S = \mathbb{R} \times \mathbb{R}^{d_{\text{in}}} \times \mathbb{R}$, where the neurons take values. For each neuron $x = (c, a, b) \in S$ we define the feature map:

$$
\mathbb{R}^{d_{\text{in}}} \ni z \mapsto f(x; z) := \tau(c)\varphi(a \cdot z + b) \in \mathbb{R},
$$

where $\tau : \mathbb{R} \to [-L, L]$ is a truncation function with the truncation threshold $L \in (0, +\infty)$. The two-layer neural network is nothing but the averaged feature map parameterized by $N$ neurons $x^1, \dots, x^N \in S$:

$$
\mathbb{R}^{d_{\text{in}}} \ni z \mapsto f^N(x^1, \dots, x^N; z) = \frac{1}{N} \sum_{i=1}^{N} f(x^i; z) \in \mathbb{R}. \quad (7)
$$

The training of neural network aims to minimize the distance between the averaged
output (7) for K data points ($z_1, \dots, z_K$) and their labels ($y_1, \dots, y_K$), that is, the
objective function of the minimization reads

$$
F_{\text{NNet}}^{N}(x^{1}, \dots, x^{N}) = \frac{N}{2K} \sum_{i=1}^{K} |y_{i} - f^{N}(x^{1}, \dots, x^{N}; z_{i})|^{2}. \quad (8)
$$

This objective function is of high dimension and non-convex, and this difficulty
motivates the recent studies, see for example [37, 16, 30] among others, to lift
the finite-dimensional optimization (8) to the space of probability measures and to
consider the following mean field optimization:

$$
\inf_{m \in \mathcal{P}_2(S)} F_{\text{NNet}}(m), \quad \text{where} \quad F_{\text{NNet}}(m) := \frac{1}{2K} \sum_{i=1}^{K} |y_i - \mathbb{E}^{X \sim m}[f(X; z_i)]|^2.
$$

The mean field loss functional $F_{\text{NNet}}$ is apparently convex. Given that the activation and truncation functions $\varphi$, $\tau$ have bounded derivatives of up to fourth order, it has been proved in [14, Proposition 4.4] that the minimum of the entropy-regularized mean field optimization problem

$$
\inf_{m \in \mathcal{P}_2(S)} F_{\text{NNet}}(m) + \frac{\sigma^2}{2} H(m | \mathcal{N}(0, \gamma^{-1}))
$$
can be attained by the invariant measure of the mean field Langevin dynamics:

$$dX_t = -D_m F(\text{Law}(X_t), X_t) dt + \sigma dW_t, \quad (9)$$

where the mean field potential functional reads

$$F(m) := F_{\text{NNet}}(m) + \frac{\sigma^2 \gamma}{2} \int |x|^2 m(dx).$$

By defining

$$\ell^0(x) := |x|^2, \quad \ell^i(x) := f(x, z_i) - y_i, \quad \text{for } i = 1, \dots, K, \\ R^{K+1} \ni \theta = (\theta_0, \theta_1, \dots, \theta_K) \mapsto \Phi(\theta) := \frac{\sigma^2 \gamma}{2} \theta_0 + \frac{1}{2K} \sum_{i=1}^{K} |\theta_i|^2,$$

we note that the mean field potential functional is of the form:

$$F(m) = \Phi(\langle \ell, m \rangle), \quad \text{where } \ell := (\ell^0, \ell^1, \dots, \ell^K),$$

as in the gradient case investigated in Sections 2.2.

In order to simulate the invariant measure of the mean field Langevin dynamics (9), one usually turns to the corresponding N-particle system:

$$dX_t^j = (-\nabla_j F_{\text{NNet}}^N(X_t^1, \dots, X_t^N) - \sigma^2 \gamma X_t^j) dt + \sigma dW_t^j, \quad \text{for } j = 1, \dots, N \quad (10)$$

The uniform-in-time propagation of chaos results obtained in [14, 43] ensure that the marginal distributions $\text{Law}(X_t^1, \dots, X_t^N)$ of the N-particle system are close to those of the mean field Langevin dynamics uniformly on the whole time horizon, and can efficiently approximate mean field invariant measure provided that $t$ and $N$ are both large enough.

Note that the N-particle system (10) is nothing but a regularized and noised gradient descent flow for $N$ neurons. In contrast, the self-interacting diffusion

$$dX_t = -\sum_{i=0}^{K} \nabla_i \Phi(Y_t^0, Y_t^1, \dots, Y_t^K) \nabla^{\ell^i}(X_t) dt + \sigma dW_t \\ = -\left(\frac{1}{K} \sum_{k=1}^{K} Y_t^k \nabla f(X_t, z_k) + \sigma^2 \gamma X_t\right) dt + \sigma dW_t, \quad (11)$$

$$dY_t^i = \lambda(\ell^i(X_t) - Y_t^i) dt, \quad \text{for } i = 1, \dots, K,$$

introduces an innovative alternative algorithm for training two-layer neural networks, in which the algorithm iterations impact only a single neuron.

**Setup.** We aim to train a neural network to approximate the non-linear elementary function $\mathbb{R}^2 \ni z = (z_1, z_2) \mapsto g(z) := \sin 2\pi z_1 + \cos 2\pi z_2$. Note that in this numerical example $d_{\text{in}} = 2$ and therefore $S = \mathbb{R}^{1+2+1}$. We draw $K$ points $\{z_i\}_{i=1}^K$ sampled according to the uniform distribution on $[0, 1]^2$ and compute the corresponding labels $y_k = g(z_k)$ to form our training data $\{z_k, y_k\}_{k=1}^K$. We fix the truncation function $\tau$ by $\tau(\theta) = (\theta \wedge 30) \vee -30$ and the sigmoid activation function $\varphi$ by $\varphi(\theta) = 1/(1 + \exp(-\theta))$. The Brownian noise has volatility $\sigma$ such that $\sigma^2/2 = 0.05$, and the regularization constant $\gamma$ is fixed at $\gamma = 0.0025$ in the
<figure><img src="image_3.png" /><figcaption>Figure 1: Averaged over 100 repetitions losses for fixed values of λ and for discrete annealing.</figcaption></figure>

experiment. The initial value $X_0 = (C_0, A_0, B_0)$, taking values in $S = \mathbb{R}^{1+2+1}$, is sampled from the normal distribution $m_0 = \mathcal{N}(0, 10^2 \times \text{Id}_{4 \times 4})$ in four dimensions. To observe the impact of the self-interacting coefficient $\lambda$, we run the simulation of (11) for different $\lambda$ equal to $4^{-p}$ for $p = 4, \dots, 8$. Furthermore, to compare with these results with fixed $\lambda$, we choose the non-increasing piecewise constant function $\lambda_a$ such that $\lambda_a(t) = 4^{-i}$ on successive intervals of length $\delta T_i = 4^i$, for $i = 2, \dots, 11$, and train the neural network along the annealing scheme:

$$dX_t = -\left(\frac{1}{K} \sum_{k=1}^{K} Y_t^k \nabla f(X_t, z_k) + \sigma^2 \gamma X_t\right) dt + \sigma dW_t, \quad (12)$$

$$dY_t^i = \lambda_a(t)(\ell^i(X_t) - Y_t^i) dt, \quad \text{for } i = 1, \dots, K.$$

We simulate both the constant and dynamic-$\lambda$ self-interacting diffusions (11), (12) by the Euler scheme, as described in Section B, on a long time horizon till the terminal time $T = 10^6$, with the discrete time step $\delta t = 0.5$.

**Results and discussions.** We repeat the simulations mentioned above for fixed $\lambda$'s and dynamic $\lambda_a$ all for 100 times and compute the averaged loss $\Phi(Y_t)$ at each time $t$ and plot its evolution in Figure 1. On the annealing scheme curve, we include triangles to visually indicate the points at which there are changes on the values of $\lambda_a$. We observe that the value of $\Phi(Y_t)$ first decays exponentially, and the speed of initial decay depending on the value of $\lambda$. More precisely, the bigger is the value of $\lambda$, the faster is the initial decay. In particular, this remains true for the annealed process as it starts from a bigger value $\lambda_a(0)$. We notice that such phenomenon is in line with our theoretical results in Theorem 1. In the long run, when fixing $\lambda$, the value of $\Phi(Y_t)$ stabilizes at a level sensible with respect to the value of $\lambda$. We notice that the smaller is the value of $\lambda$, the better is the final performance.
These facts are in line with our discovery in Theorems 1 and 3. Finally, the training exhibits the best performance when implementing the discrete annealing. The loss $\Phi(Y_t)$ continues to decrease as the value of $\lambda_a(t)$ diminishes.

# 4 A Curie–Weiss model

In this section, we present and study a simple Curie–Weiss model, i.e., a mean field model of ferromagnets, which has possibly multiple invariant measures. In particular, we show that in this case, the last claim of Theorem 3 provides informations on the concentration of the self-interacting stationary measure.

Let $\ell_0 : \mathbb{R} \to \mathbb{R}$ be a smooth, odd, increasing function in $C^1 \cap W^{1,\infty}$. For a probability measure $m \in \mathcal{P}(\mathbb{R})$, consider the mean field energy

$$F(m) = -\frac{1}{2} \left( \int \ell_0(x) m(dx) \right)^2 + \frac{1}{2} \int x^2 m(dx).$$

By setting

$$\ell(x) := (\ell_0(x), |x|^2/2)^\top, \\ \Phi(y_0, y_1) := -\frac{1}{2}|y_0|^2 + y_1,$$

we have $F(m) = \Phi(\langle \ell, m \rangle)$. So we are in the cylindrical setting of Assumption 2 with the range set of $\ell$ being defined by

$$\mathcal{K} := [-\|\ell_0\|_{L^\infty}, \|\ell_0\|_{L^\infty}] \times \mathbb{R}.$$

Moreover, as the corresponding intrinsic derivative, or the drift, writes

$$b(m, x) = -D_m F(m, x) = \langle \ell_0, m \rangle \ell'_0(x) - x,$$

we can verify Assumption 1 with the modulus of continuity

$$\omega(r) = \sup_{x,x' \in \mathbb{R}: |x-x'| \le r} |\ell_0(x) - \ell_0(x')|.$$

This implies, by Corollary 2, that the self-interacting process has a unique invariant measure, which we denote by $\rho^\lambda = \rho$. Arguing as in the proof of Proposition 4, we can verify the conditions of Theorem 3 for the choice of $\Phi$, $\ell$ and $\rho$ above. Especially, the probability measures

$$\hat{m}_{(y_0, y_1)}(dx) \propto \exp \left( -\frac{1}{2} |x|^2 + y_0 \ell_0(x) \right) dx$$

satisfy a uniform LSI thanks to the boundedness of $y_0$ and the Holley–Stroock perturbation lemma.

Before applying Theorem 3, we first give a characterization of the invariant measure for the McKean–Vlasov dynamics (1). We discuss a special case and then move to general discussions.

1. If $\|\ell'\|_{L^\infty} < 1$, then by the result of [44], we already know that the McKean–Vlasov dynamics (1) has a unique invariant measure, which is the centered Gaussian of variance $1/2$, i.e., $\mathcal{N}(0, 1/2)$. This case corresponds to the weak interaction regime studied in [21].
2. In the general case where $\|\ell'\|_{L^\infty}$ is not restricted, the invariant measures of (1) correspond to fixed points of the one-dimensional mapping

$$ \mathbb{R} \ni y_0 \mapsto \Pi_0(y_0) = \hat{y}_0 \in \mathbb{R} $$

given by

$$ \Pi_0(y_0) = \frac{\int \ell_0(x) \exp(2y_0\ell_0(x) - |x|^2) \, dx}{\int \exp(2y_0\ell_0(x) - |x|^2) \, dx} $$

That is to say, if $y_0$ satisfies $\Pi_0(y_0) = y_0$, then the probability measure proportional to $\exp(2y_0\ell_0(x) - |x|^2) \, dx$ is invariant to (1); and vice versa. Due to the oddness of $\ell_0$, the mapping $\Pi_0$ is odd. In particular, $\Pi_0(0) = 0$ and we know that $\mathcal{N}(0, 1/2)$ is always invariant.

If $\Pi'_0(0) > 1$, then by the fact that $\|\Pi_0\|_{L^\infty} \le \|\ell_0\|_{L^\infty} < \infty$ and the intermediate value theorem, there exists $y_0 > 0$ such that $\Pi_0(y_0) = y_0$ and $\Pi_0(-y_0) = -y_0$. That is to say there exists at least three invariant measures, two of which are, in physicists' language, "symmetry-breaking" phases, and the centered Gaussian measure is the "symmetry-preserving" phase. Moreover, the centered measure corresponding to $y_0 = 0$ should be unstable as it is a local maximum point for the free energy landscape $y_0 \mapsto \frac{1}{2}y_0^2 - \frac{1}{2}\log \int \exp(2y_0\ell_0(x) - |x|^2) \, dx$.

We now turn to the study of the stationary measure $\rho^\lambda$ of the self-interacting process (4). Since $\Phi(y) = -|y_0|^2/2 + y_1$ is linear in its second coordinate, the last claim of Theorem 3 implies that

$$ \int |\hat{y}_0 - y_0|^2 \rho^2(dy) = - \int (\hat{y} - y)^\top \nabla^2 \Phi(y) (\hat{y} - y) \rho^2(dy) \le C\lambda, $$

where $C$ is a constant depending only on $\ell_0$ and as we recall, $y_0$ is the first coordinate of $y$. In other words, the stationary measure $\rho^\lambda$ solves the self-consistency equation

$$ \hat{y}_0 = y_0 $$

up to an error of order $O(\sqrt{\lambda})$. Denote the set of fixed points of $y_0 \mapsto \Pi_0(y_0) = \hat{y}_0$ by $S$. If the set $S$ is finite, and if for each $s \in S$ we have

$$ \Pi'_0(s) \ne 1, $$

then there exists $c > 0$ such that for all $y_0 \in \mathbb{R}$,

$$ |\hat{y}_0 - y_0| \ge c \min(d(y_0, S), 1), $$

where $d(\cdot, S)$ indicates the distance to the set $S$. In this case, we have

$$ \int \min(d(y_0, S), 1)^2 \rho^2(dy) = O(\lambda). $$

That is to say, for small $\lambda$, the random variable $Y = (Y_0, Y_1)$, distributed as the second component of the stationary measure $\rho^\lambda$, has $Y_0$ concentrated near the solutions to the self-consistency equation, which correspond to invariant measures of the McKean-Vlasov dynamics. However, the last claim of Theorem 3 is not sufficient to show that $Y_0$ is only concentrated around, or in a way "selects", the physically stable phases that are minimizers of a free energy functional. We refer readers to [5] for qualitative results on such selection mechanism.
# 5 Proof of Theorem 1

We first note that as the metric space $(\mathcal{P}(\mathbb{R}^d), W_\omega)$ is separable, we do not have measurability issues. We refer readers to [35, Chapter 2] for details.

Recall that the self-interacting dynamics (4) writes

$$dX_t = b(m_t, X_t) dt + dB_t,$$

$$d m_t = \lambda(-m_t + \delta_{X_t}) dt$$

and similarly for the other dynamics $(X', m')$ driven by another Brownian motion $B'$. Fix $n \in \mathbb{N}$. Let $\text{rc} : \mathbb{R}^d \times \mathbb{R}^d \to [0, 1]$ be a Lipschitz continuous function such that $\text{sc} := \sqrt{1 - \text{rc}^2}$ is also Lipschitz continuous and

$$\text{rc}(x, x') = \begin{cases} 1 & \text{if } |x - x'| \ge 2n^{-1}, \\ 0 & \text{if } |x - x'| \le n^{-1}. \end{cases}$$

We couple the two dynamics $(X, m)$, $(X', m')$ by

$$d B_t = \text{rc}(X_t, X'_t) d B_t^{\text{rc}} + \text{sc}(X_t, X'_t) d B_t^{\text{sc}},$$

$$d B'_t = \text{rc}(X_t, X'_t)(1 - 2e_t e_t^\top) d B_t^{\text{rc}} + \text{sc}(X_t, X'_t) d B_t^{\text{sc}},$$

where $B^{\text{rc}}$, $B^{\text{sc}}$ are independent Brownian motions and $e_t$ is the $d$-dimensional vector defined by

$$e_t = \begin{cases} \frac{X_t - X'_t}{|X_t - X'_t|} & \text{if } X_t \neq X'_t, \\ (1, 0, \dots, 0)^\top & \text{otherwise.} \end{cases}$$

Subtracting the dynamical equations of $X$, $X'$ and denoting $\delta X = X - X'$, we obtain

$$d\delta X_t = (b(m_t, X_t) - b(m'_t, X'_t)) dt + 2 \text{rc}(X_t, X'_t) e_t dW_t,$$

where $W_t := \int_0^t e_t^\top d B_t^{\text{rc}}$ and is a one-dimensional Brownian motion by Lévy's characterization. The absolute value of the semimartingale $\delta X_t$ admits the decomposition $d|\delta X_t| = dA_t^{|\delta X|} + dM_t^{|\delta X|}$ with

$$dA_t^{|\delta X|} \le -|\delta X_t|\kappa(|\delta X_t|) dt + LW_\omega(m_t, m'_t) dt,$$

$$dM_t^{|\delta X|} = 2 \text{rc}(X_t, X'_t) dW_t.$$

For the $m$ variable, we have

$$d(m_t - m'_t) = -\lambda(m_t - m'_t) dt + \lambda(\delta_{X_t} - \delta_{X'_t}) dt.$$

Thus, by considering the (random) $W_\omega$-optimal coupling at each $t$, we get

$$dW_\omega(m_t, m'_t) \le -\lambda W_\omega(m_t, m'_t) dt + \lambda\omega(|\delta X_t|) dt.$$

Hence the process

$$r_t = |\delta X_t| + \frac{2L}{\lambda} W_\omega(m_t, m'_t)$$

admits the decomposition $dr_t = dA_t^r + dM_t^r$ with

$$dA_t^r \le (-|\delta X_t|\kappa(|\delta X_t|) + 2L\omega(|\delta X_t|) - LW_\omega(m_t, m'_t)) dt,$$

$$dM_t^r = 2 \text{rc}(X_t, X'_t) dW_t.$$
Let $f : [0, \infty) \to [0, \infty)$ be a $C^2$ function to be determined in the following. We define $\rho_t = f(r_t)$. The process $\rho_t$ admits the decomposition $d\rho_t = dA_t^\rho + dM_t^\rho$ with

$$ 
\begin{aligned}
dA_t^\rho &\le (-|\delta X_t|\kappa(|\delta X_t|) + 2L\omega(|\delta X_t|) - LW_\omega(m_t, m'_t))f'_-(r_t) dt \\
&\quad + 2 \operatorname{rc}(X_t, X'_t)^2 f''(r_t) dt \\
&\le -r_t \tilde{K}(r_t) f'_-(r_t) dt + 2 \operatorname{rc}(X_t, X'_t)^2 f''(r_t) dt
\end{aligned}
$$

for the function $\tilde{K} : (0, \infty) \to \mathbb{R}$ defined by

$$ 
\begin{aligned}
\tilde{K}(r) &:= \inf_{\substack{x+2L\lambda^{-1}y=r \\ x,y>0}} \frac{x\kappa(x) - 2L\omega(x) + Ly}{r} \\
&\ge \inf_{\substack{x+2L\lambda^{-1}y=r \\ x,y>0}} \frac{\kappa_0 x + Ly - M_b - 2LM_\omega}{r} \\
&\ge \min\left(\kappa_0, \frac{\lambda}{2}\right) - \frac{M_b + 2LM_\omega}{r} \\
&:= K_0 - \frac{M}{r} =: K(r).
\end{aligned}
$$

Thus, we have shown

$$ dA_t^\rho \le -r_t K(r_t) f'_-(r_t) dt + 2 \operatorname{rc}(X_t, X'_t)^2 f''(r_t) dt. $$

Following Du et al. [22], we choose the function $f : [0, \infty) \to [0, \infty)$ by requiring

$$ f'(r) = \frac{1}{2} \int_{r}^{\infty} s \exp \left( -\frac{1}{2} \int_{r}^{s} \tau K(\tau) \, d\tau \right) \, ds. $$

and $f(0) = 0$. The function $f$ is indeed $C^2$ continuous, and, according to [22, Lemma 5.1], it is also concave and satisfies

$$ 2f''(r) - rK(r)f'(r) + r = 0 $$

and

$$ \frac{1}{K_0} \le f'(r), \quad \frac{f(r)}{r} \le f'(0) $$

for all $r > 0$. Plugging in the expression for $K$, we obtain

$$ 
\begin{aligned}
f'(0) &= \frac{1}{2} \int_{0}^{\infty} s \exp\left(-\frac{K_0 s^2}{4} + \frac{Ms}{2}\right) ds \\
&= \frac{1}{2} \exp\left(\frac{M^2}{4K_0}\right) \int_{0}^{\infty} s \exp\left(-\frac{K_0\left(s - \frac{M}{K_0}\right)^2}{4}\right) ds \\
&= \frac{1}{2} \exp\left(\frac{M^2}{4K_0}\right) \int_{-M/\sqrt{2K_0}}^{\infty} \left(\frac{2t}{K_0} + M \frac{2^{1/2}}{K_0^{3/2}}\right) \exp\left(-\frac{t^2}{2}\right) dt \\
&\le \frac{1}{2} \exp\left(\frac{M^2}{4K_0}\right) \left(\frac{2}{K_0} \exp\left(-\frac{M^2}{4K_0}\right) + 2M \frac{\pi^{1/2}}{K_0^{3/2}}\right) \\
&\le \frac{1}{K_0} + \frac{2M}{K_0^{3/2}} \exp\left(\frac{M^2}{4K_0}\right).
\end{aligned}
$$
For $|\delta X_t| \ge 2n^{-1}$, we have $\mathrm{rc}_t(X_t, X'_t) = 1$, and by the previous construction,

$$ \mathrm{d}f(r_t) = \mathrm{d}A_t^\rho + \mathrm{d}M_t^\rho \le -r_t \mathrm{d}t + \mathrm{d}M_t^\rho \le -\frac{f(r_t)}{K_0^{-1} + 2MK_0^{-3/2} \exp(M^2/4K_0)} \mathrm{d}t + \mathrm{d}M_t^\rho. $$

For $|\delta X_t| < 2n^{-1}$, we proceed differently. Let $\omega_b$ denote the uniform modulus of equicontinuity of the vector fields $x \mapsto b(m, x)$. The absolute continuous part of $r_t$ satisfies

$$ \begin{aligned} \mathrm{d}A_t^r &\le (\omega_b(2n^{-1}) + 2L\omega(2n^{-1}) - LW_\omega(m_t, m'_t)) \mathrm{d}t \\ &= -LW_\omega(m_t, m'_t) \mathrm{d}t + o(1) \mathrm{d}t \\ &\le -\frac{\lambda}{2}r_t \mathrm{d}t + o(1) \mathrm{d}t, \end{aligned} $$

where $o(1)$ denotes a term that tends to 0 when $n \to \infty$.

Now we combine the two cases. Define a sequence of functions by

$$ f_n(r) = \begin{cases} f(r), & \text{if } r \ge 2n^{-1}, \\ \frac{f(2n^{-1})}{2n^{-1}}r, & \text{if } r < 2n^{-1}. \end{cases} $$

Each function $f_n$ is concave and satisfies, by the arguments above,

$$ \mathrm{d} \mathrm{E}[f_n(r_t)] \le -c \mathrm{E}[f_n(r_t)] \mathrm{d}t + o(1) \mathrm{d}t, $$

where

$$ \begin{aligned} c &:= \min\left(\frac{1}{K_0^{-1} + 2MK_0^{-3/2} \exp(M^2/4K_0)}, \frac{\lambda}{2}\right) \\ &= \frac{1}{K_0^{-1} + 2MK_0^{-3/2} \exp(M^2/4K_0)}. \end{aligned} $$

Applying Grönwall's lemma, we obtain

$$ \mathrm{E}[f_n(r_t)] \le e^{-ct} \mathrm{E}[f_n(r_0)] + o(1). $$

Since $\lim_{n \to \infty} \mathrm{E}[f_n(r_t)] = \mathrm{E}[f(r_t)]$ by dominated convergence, taking the limit $n \to \infty$ completes the proof. □

# 6 Proof of Theorem 3

This section consists of four subsections. We show a series of intermediate, and sometimes technical, lemmas and propositions in the first three subsections before proving the theorem in the last subsection.

## 6.1 Elliptic equation for the stationary measure

We first note that the stationary measure $P$ solves a partial differential equation in the following weak sense.
**Proposition 5.** Let $C_b^{1,2}$ denote the class of functionals $\phi: \mathcal{P}(\mathbb{R}^d) \times \mathbb{R}^d \to \mathbb{R}$ with continuous spatial derivatives $\nabla_x\phi, \nabla_x^2\phi$ and a bounded linear functional derivative $\frac{\delta\phi}{\delta m}: \mathcal{P}(\mathbb{R}^d) \times \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ satisfying

$$ \forall m \in \mathcal{P}(\mathbb{R}^d), \forall x \in \mathbb{R}^d, \quad |\nabla_x\phi(m, x)| + |\nabla_x^2\phi(m, x)| \le C, $$

$$ \forall m \in \mathcal{P}(\mathbb{R}^d), \forall x, x' \in \mathbb{R}^d, \quad \left| \frac{\delta\phi}{\delta m}(m, x, x') \right| \le C, $$

for some $C \ge 0$. Under the settings of the theorem, let $(X, m)$ be a random variable distributed as the stationary measure $P$. Then we have, for all $\phi \in C_b^{1,2}$,

$$ \mathbb{E} \left[ \frac{1}{2} \Delta_x \phi(m, X) + b(m, X) \cdot \nabla_x \phi(m, X) + \lambda \int \frac{\delta \phi}{\delta m}(m, X, x')(\delta_X - m)(dx') \right] = 0. $$

We omit the proof of the proposition, which follows directly from expanding the difference $\mathbb{E}[\phi(m_t, X_t)] - \mathbb{E}[\phi(m_0, X_0)]$ by Itô-type calculus.

The infinite-dimensional nature of the PDE above makes its analysis difficult, and in the following we approach the problem by studying a finite-dimensional projection of it. Under Assumption 2, define the functions

$$ \beta(y, x) = -\nabla\Phi(y) \cdot \nabla\ell(x) = -\sum_{\nu=1}^{D} \nabla_\nu\Phi(y)\nabla\ell^\nu(x), $$

$$ V(y, x) = \nabla\Phi(y) \cdot \ell(x) = \sum_{\nu=1}^{D} \nabla_\nu\Phi(y)\ell^\nu(x). $$

They verify $\beta(y, x) = -\nabla_x V(y, x)$. Note that, if $m_y$ is a measure satisfying $\int \ell(x)m_y(dx) = y$, then we have

$$ \beta(y, x) = b(m_y, x), $$

$$ V(y, x) = \frac{\delta F}{\delta m}(m_y, x). $$

Let $(X, m)$ be distributed as the stationary measure $P$ and consider the random variable $Y = \langle \ell, m \rangle = \int \ell(x)m(dx)$ valued in $\mathcal{K}$. Applying Proposition 5 to functionals of the following form:

$$ \phi(m, x) = \varphi \left( x, \int \ell(x')m(dx') \right), $$

where $\varphi \in \mathcal{D}(\mathbb{R}^d \times \mathbb{R}^D)$, we get that the joint distribution $\rho = \rho^\lambda := \text{Law}(X, Y) \in \mathcal{P}(\mathbb{R}^d \times \mathcal{K})$ satisfies the stationary degenerate Fokker-Planck equation

$$ \frac{1}{2} \Delta_x \rho - \nabla_x \cdot (\beta(y, x) \rho) - \lambda \nabla_y \cdot ((\ell(x) - y) \rho) = 0. \quad (13) $$

in the sense of distributions. The fact that $P$ has finite fourth moment implies that its projection $\rho$ satisfies the following moment condition:

$$ \int (|x|^4 + |y|^2) \rho(dx dy) < \infty. \quad (14) $$

From the Fokker-Planck equation (13), we get the following result.
**Lemma 6.** Under the setting of the theorem, for every function $\varphi \in C^1(\mathbb{R}^D; \mathbb{R})$ whose gradient $\nabla\varphi$ is of linear growth, we have

$$\iint \nabla\varphi(y) \cdot (\ell(x) - y)\rho(dx dy) = 0.$$

*Proof of Lemma 6.* Since its gradient $\nabla\varphi$ is of linear growth, the function $\varphi$ is of quadratic growth. Thanks to the fact that $\rho$ satisfies the moment condition (14), we can take the duality with $\varphi$ in the static Fokker-Planck equation (13), from which the desired result follows. $\square$

## 6.2 Entropy estimate

In this subsection, we show an entropy estimate for the stationary measure $\rho$ by studying the Fokker-Planck equation (13).

Denote the first and second marginal distributions of $\rho$ by $\rho^1, \rho^2$ respectively. For $\rho^2$-almost all $y \in \mathbb{R}^D$, denote also the conditional distribution of the first variable by $\rho^{1|2}(\cdot|y) : \mathcal{B}(\mathbb{R}^d) \to \mathbb{R}$. Define

$$\hat{\rho}(dx dy) := \hat{m}_y(dx)\rho^2(dy) = \frac{1}{Z_y} \exp(-2V(y, x)) dx \rho^2(dy), \quad (15)$$

for $Z_y = \int \exp(-2V(y, x)) dx$. Recall that $\hat{m}_y$ is the probability measure on $\mathbb{R}^d$ that satisfies a uniform LSI according to Assumption 3.

**Proposition 7.** Under the setting of the theorem, we have

$$H(\rho|\hat{\rho}) \le \frac{C_{\text{LS}}}{2} \left( D + 2 \iint (\ell(x) - y)^{\top} \nabla^2 \Phi(y) (\ell(x) - y) \rho(dx dy) \right) \lambda,$$

where $\hat{\rho}$ is the measure defined by (15).

For convenience, we set

$$I := \iint (\ell(x) - y)^{\top} \nabla^2 \Phi(y) (\ell(x) - y) \rho(dx dy), \quad (16)$$

so the claim of the proposition reads

$$H(\rho|\hat{\rho}) \le \frac{C_{\text{LS}}}{2} (D + 2I)\lambda.$$

*Proof of Proposition 7.* Let $g^\epsilon$ be the Gaussian kernel in $D$ dimensions:

$$g^\epsilon(y) = (2\pi\epsilon)^{-D/2} \exp\left(-\frac{|y|^2}{2\epsilon}\right).$$

We define the partial convolution $\rho^\epsilon = \rho \star_y g^\epsilon$, and according to (13), it solves the non-degenerate elliptic equation

$$\frac{1}{2} \Delta_x \rho^\epsilon - \nabla_x \cdot (\beta^\epsilon(y, x) \rho^\epsilon) - \ell(x) \cdot \nabla_y \rho^\epsilon + \lambda \epsilon \Delta_y \rho^\epsilon + \lambda \nabla_y \cdot (y \rho^\epsilon) = 0, \quad (17)$$

where $\beta^\epsilon$ is defined by

$$\beta^\epsilon = \frac{(\beta \rho) \star_y g^\epsilon}{\rho^\epsilon}.$$
Indeed, we have

$$
\begin{align*}
(y\rho) \star_y g^\epsilon &= \int y' \rho(x', y') g^\epsilon(y - y') \, dy' \\
&= \int (y' - y + y) \rho(x', y') g^\epsilon(y - y') \, dy' \\
&= \varepsilon \int \rho(x', y') \nabla_y g^\epsilon(y - y') \, dy' + y \rho^\epsilon \\
&= \varepsilon \nabla_y \rho^\epsilon + y \rho^\epsilon.
\end{align*}
$$

Thus, $(\nabla_y \cdot (y\rho)) \star_y g^\epsilon = \nabla_y \cdot ((y\rho) \star_y g^\epsilon) = \varepsilon \Delta_y \rho^\epsilon + \nabla_y \cdot (y\rho^\epsilon)$. By [9, Lemma 3.1.1], we have $\|\beta^\epsilon\|_{L^2(\rho^\epsilon)} \le \|\beta\|_{L^2(\rho)} < \infty$. Then we can apply [9, Theorem 3.1.2] to (17) and obtain that $\rho^\epsilon \in W^{1,1}(\mathbb{R}^{d+D})$ and satisfies

$$
\begin{align*}
\frac{1}{2} \iint \frac{|\nabla_x \rho^\epsilon|^2}{\rho^\epsilon} + \lambda \epsilon \iint \frac{|\nabla_y \rho^\epsilon|^2}{\rho^\epsilon} \\
&= \iint \nabla_x \rho^\epsilon \cdot \beta^\epsilon + \lambda \iint \nabla_y \rho^\epsilon \cdot \ell - \lambda \iint \nabla_y \rho^\epsilon \cdot y.
\end{align*}
$$

As the function $\ell$ depends only on the $x$ argument, for the second term we have

$$
\iint \nabla_y \rho^\epsilon \cdot \ell = \iint \left( \iint \nabla_y \rho^\epsilon(x, y) dy \right) \ell(x) dx = \int 0 \cdot \ell(x) dx = 0.
$$

where the first equality is due to Fubini and the second to the fact that $\nabla \rho^\epsilon \in L_x^1(L_y^1)$. For the last term, similarly, since the function $((x, y) \mapsto \rho^\epsilon(x, y)y) \in W^{1,1}$, we have $\iint \nabla_y \cdot (\rho^\epsilon y) = 0$ and therefore,

$$
- \iint \nabla_y \rho^\epsilon \cdot y = \iint (\nabla_y \cdot y) \rho^\epsilon = D.
$$

That is to say, we have established

$$
\frac{1}{2} \iint \frac{|\nabla_x \rho^\epsilon|^2}{\rho^\epsilon} + \lambda \epsilon \iint \frac{|\nabla_y \rho^\epsilon|^2}{\rho^\epsilon} = \iint \nabla_x \rho^\epsilon \cdot \beta^\epsilon + \lambda D. \quad (18)
$$

This equality implies a uniform-in-$\epsilon$ bound on $\iint |\nabla_x \rho^\epsilon|^2 / \rho^\epsilon$ by Cauchy-Schwarz.
Using a sequence of functions in $C_c^\infty$ approaching $V(y, x)$ in (13), we also get

$$
\begin{align}
\iint \frac{1}{2} \beta \cdot \nabla_x \rho^\epsilon - \iint \lambda \epsilon \nabla_y V \cdot \nabla_y \rho^\epsilon - \iint \beta \cdot \beta^\epsilon \rho^\epsilon \nonumber \\
\qquad + \iint \lambda \nabla_y V(y, x) \cdot (\ell(x) - y) \rho^\epsilon (dx dy) = 0. \tag{19}
\end{align}
$$

The second term of (19) is upper bounded by

$$
\begin{align*}
\lambda \epsilon \iint |\nabla_y V \cdot \nabla_y \rho^\epsilon| &\le \lambda \epsilon \|\nabla_y V\|_{L^2(\rho^\epsilon)} \left( \iint \frac{|\nabla_y \rho^\epsilon|^2}{\rho^\epsilon} \right)^{1/2} \\
&\le \sqrt{\lambda \epsilon} \|\nabla_y V\|_{L^2(\rho^\epsilon)} \sqrt{\frac{1}{2} \iint \frac{|\nabla_x \rho^\epsilon|^2}{\rho^\epsilon} + \frac{1}{2} \|\beta^\epsilon\|_{L^2(\rho^\epsilon)}^2 + \lambda D},
\end{align*}
$$
where the second inequality is due to (18). Using the uniform-in-$\epsilon$ bound of $\iint |\nabla_x \rho^\epsilon|^2/\rho^\epsilon$, we obtain that the second term of (19) converges to 0 when $\epsilon \to 0$. The third term of (19) satisfies

$$\iint \beta \cdot \beta^\epsilon \rho^\epsilon = \iint \beta \cdot ((\beta \rho) \star g^\epsilon) = \iint (\beta \star g^\epsilon) \cdot \beta \rho \to \iint |\beta|^2 \rho$$

when $\epsilon \to 0$. Hence, by (18) and (19), we have

$$\begin{aligned}
\frac{1}{2} \iint |\nabla_x \log \rho^\epsilon - 2\beta|^2 \rho^\epsilon (dx dy) & \\
& \le \lambda D + \iint \nabla_x \rho^\epsilon \cdot \beta^\epsilon + \iint \beta \cdot \nabla_x \rho^\epsilon - 2 \iint \nabla_x \rho^\epsilon \cdot \beta \\
& \quad + 2\lambda \iint \nabla_y V(y,x) \cdot (\ell(x) - y) \rho^\epsilon (dx dy) + o(1) \\
& = \lambda D + 2\lambda \iint \nabla_y V(y,x) \cdot (\ell(x) - y) \rho^\epsilon (dx dy) + o(1),
\end{aligned}$$

where the last equality is due to the fact that

$$\left|\iint \nabla_x \rho^\epsilon \cdot (\beta^\epsilon - \beta)\right| \le \left(\iint \frac{|\nabla_x \rho^\epsilon|^2}{\rho^\epsilon}\right)^{1/2} \left(\iint |\beta - \beta^\epsilon|^2 \rho^\epsilon\right)^{1/2}$$

and

$$\begin{aligned}
\iint |\beta - \beta^\epsilon|^2 \rho^\epsilon &= \iint |\beta|^2 \rho^\epsilon - 2 \iint \beta \cdot \beta^\epsilon \rho^\epsilon + \iint |\beta^\epsilon|^2 \rho^\epsilon \\
&= \iint (|\beta|^2 \star_y g^\epsilon) \rho - 2 \iint (\beta \star_y g^\epsilon) \cdot \beta \rho + \iint |\beta^\epsilon|^2 \rho^\epsilon \to 0,
\end{aligned}$$

when $\epsilon \to 0$ by previous arguments and [9, Lemma 3.1.1]. We also have

$$\begin{aligned}
\iint \nabla_y V(y,x) \cdot (\ell(x) - y) \rho^\epsilon (dx dy) &= \iint (\nabla_y V(y,x) \cdot (\ell(x) - y)) \star_y g^\epsilon \rho (dx dy) \\
&\to \iint \nabla_y V(y,x) \cdot (\ell(x) - y) \rho (dx dy)
\end{aligned}$$

when $\epsilon \to 0$. Finally, by Lemma 6 for the function $\varphi(y) = y \cdot \nabla \Phi(y) - \Phi(y)$, we have

$$\begin{aligned}
\iint \nabla_y V(y,x) \cdot (\ell(x) - y) \rho(dx dy) & \\
&= \iint \ell(x)^\top \nabla^2 \Phi(y) (\ell(x) - y) \rho(dx dy) \\
&= \iint (\ell(x) - y)^\top \nabla^2 \Phi(y) (\ell(x) - y) \rho(dx dy) = I,
\end{aligned}$$

where the last equality is exactly the definition (16) of $I$. Thus, we have shown

$$\frac{1}{2} \iint \left| \nabla_x \log \frac{\rho^\epsilon(x, y)}{\hat{m}_y(x)} \right|^2 \rho^\epsilon (dx dy) \le (D + 2I)\lambda + o(1).$$
Note that by the lower semicontinuity of (partial) Fisher information,

$$ \liminf_{\epsilon \to 0} \iint \left| \nabla_x \log \frac{\rho^{\epsilon}(x, y)}{\hat{m}_y(x)} \right|^2 \rho^{\epsilon}(dx dy) \ge \iint \left| \nabla_x \log \frac{\rho(x, y)}{\hat{m}_y(x)} \right|^2 \rho(dx dy). $$

We refer readers to [12, Proof of Lemma A.1] for details. Taking the limit $\epsilon \to 0$, we obtain

$$ \frac{1}{2} \iint \left| \nabla_x \log \frac{\rho(x, y)}{\hat{m}_y(x)} \right|^2 \rho(dx dy) \le (D + 2I)\lambda. $$

Since $\rho^2$ is supported on $\mathcal{K}$, for $\rho^2$-almost all $y \in \mathbb{R}^D$, we have the following by the uniform LSI for $(\hat{m}_y)_{y \in \mathcal{K}}$:

$$ H(\rho|\hat{\rho}) = \int H(\rho^{1/2}(\cdot|y)|\hat{m}_y)\rho^2(dy) \le \frac{C_{\text{LS}}}{4} \iint \left| \nabla_x \log \frac{\rho(x, y)}{\hat{m}_y(x)} \right|^2 \rho(dx dy), $$

which completes the proof. $\square$

<figure><img src="image_4.png" /></figure>

**Remark 8.** If we formally integrate the static Fokker-Planck equation (13) with $\log(\rho/\hat{\rho})$ and integrate by parts, we obtain

$$ \frac{1}{2} \iint \left| \nabla_x \log \frac{\rho(x, y)}{\hat{m}_y(x)} \right|^2 \rho(dx dy) = \lambda D + 2\lambda \iint \ell(x)^{\top} \nabla^2 \Phi(y) (\ell(x) - y) \rho(dx dy). \quad (20) $$

However, the equality must not hold in all circumstances. Indeed, if one artificially increases the dimension of $\ell$ and $\Phi$ by defining the new functions

$$ \begin{aligned} \tilde{\Phi}(y_0, y_1) &= \Phi(y_0), \\ \tilde{\ell}(x) &= (\ell(x), 0), \end{aligned} $$

the right hand side of (20) increases while the left hand side stays unchanged. This phenomenon is caused by the fact that the equation (13) is degenerate elliptic and lacks Laplacian in the y directions. To illustrate this effect, consider the first-order equation

$$ \nabla_y \cdot (y\rho) = 0 $$

in d dimensions. This equation has a probability solution $\rho = \delta_0$, the Dirac mass at the origin. Formally integrating the equation with $\log \rho$ and integrating by parts, we have

$$ 0 = \int \log \rho \nabla_y \cdot (y\rho) = - \int \frac{\nabla_y \rho}{\rho} \cdot y\rho = - \int \nabla_y \rho \cdot y = \int \rho \nabla_y \cdot y = \int \rho d = d, $$

which is absurd.

To complete the entropy estimate, we provide in the following upper bounds for the integral $I$.

**Proposition 8.** Under the setting of the theorem, the integral $I$ in (16) satisfies the upper bound:

$$ I \le 12M_2 \left( dC_{\text{LS}} + 2W_2^2 (\text{Law}(X), m_*) \right), \quad (21) $$
where $X$ is the first component of the random variable $(X, m)$ following the stationary distribution $P = P^\lambda$, and $m_*$ is the unique invariant measure of the McKean-Vlasov process (1). If additionally $\nabla^2\Phi$ is convex and the quantity

$$M_1 := \sup_{x \in \mathbb{R}^d, y \in \mathbb{R}^D} \ell(x)^\top \nabla^2 \Phi(y) \ell(x)$$

is finite, then we have the alternative upper bound:

$$I \le M_1. \qquad (22)$$

*Proof of Proposition 8.* First let us treat the simpler case where $M_1 < \infty$. By applying Lemma 6 to the function $\varphi(y) = y \cdot \nabla\Phi(y) - \Phi(y)$, we get

$$I = \iint (\ell(x)^\top \nabla^2 \Phi(y) \ell(x) - y^\top \nabla^2 \Phi(y) y) \rho(dx dy) \\ \le \iint \ell(x)^\top \nabla^2 \Phi(y) \ell(x) \rho(dx dy) \le M_1.$$

So the second claim of the proposition is proved.

Without the assumption $M_1 < \infty$, we note that, for the second-order functional derivative

$$\frac{\delta^2 F}{\delta m^2}(m, x', x'') = \ell(x'')^\top \nabla^2 \Phi(\langle \ell, m \rangle) \ell(x),$$

we have

$$I = \mathbb{E} \left[ \iint \frac{\delta^2 F}{\delta m^2}(m, x', x'')(\delta_X - m)(dx')(\delta_X - m)(dx'') \right],$$

where $(X, m)$ is a random variable following the stationary distribution $P$, and $(X, \langle \ell, m \rangle)$ has the distribution $\rho$. We observe

$$\begin{aligned} |D_m^2 F(m, x', x'')| &= |\nabla_{x'} \nabla_{x''} \frac{\delta^2 F}{\delta m^2}(m, x', x'')| \\ &= |\nabla \ell(x'')^\top \nabla^2 \Phi(\langle \ell, m \rangle) \nabla \ell(x')| \le M_2. \end{aligned}$$

Then, by applying Lemma 10 in Appendix A to a sequence of $C^2$ functions approaching $(x', x'') \mapsto \frac{\delta^2 F}{\delta m^2}(m, x', x'')$, we get

$$\begin{aligned} I &\le M_2 \mathbb{E}[W_2^2(\delta_X, m)] \\ &= M_2 \mathbb{E} \left[ \int |X - x'|^2 m(dx') \right] \\ &\le 2M_2 \mathbb{E} \left[ |X - \mathbb{E}[X]|^2 + \int |x' - \mathbb{E}[X]|^2 m(dx') \right]. \end{aligned}$$

Let $\phi : \mathcal{P}_2(\mathbb{R}^d) \to \mathbb{R}$ be the functional defined by

$$\phi(m) = \int |x' - \mathbb{E}[X]|^2 m(dx').$$

We consider the sequence of $C_b^1$ “soft cut-offs” that approach $\phi$:

$$\phi_n(m) = \sum_{i=1}^{d} \int n^2 \tanh^2 \left( \frac{x'^i - \mathbb{E}[X^i]}{n} \right) m(dx'), \quad \text{for } n \in \mathbb{N}.$$
Then, by applying the sequence $\phi_n$ to Proposition 5 of stationary measure and taking the limit $n \to \infty$, we get

$$0 = \lambda \mathbb{E} \left[ \int \frac{\delta \phi}{\delta m} (m, x')(\delta_X - m)(dx') \right] \\ = \lambda \mathbb{E} \left[ |X - \mathbb{E}[X]|^2 \right] - \lambda \mathbb{E} \left[ \int |x' - \mathbb{E}[X]|^2 m(dx') \right].$$

Thus, we have derived

$$I \le 4M_2 \mathbb{E}[|X - \mathbb{E}[X]|^2] =: 4M_2 \mathrm{Var} X,$$

where $\mathrm{Var} X$ denotes the sum of the variances of each component of the random vector $X$. It remains only to find an upper bound for $\mathrm{Var} X$. Note that, using the definition of Wasserstein distance and the triangle inequality, and letting $X_*$ be distributed as $m_*$, we get

$$\mathrm{Var} X = W_2^2(\mathrm{Law}(X), \delta_{\mathrm{E}[X]}) \\ \le 3(W_2^2(\mathrm{Law}(X), m_*) + W_2^2(m_*, \delta_{\mathrm{E}[X_*]}) + W_2^2(\delta_{\mathrm{E}[X_*]}, \delta_{\mathrm{E}[X]})) \\ \le 3(\mathrm{Var} X_* + 2W_2^2(\mathrm{Law}(X), m_*)),$$

while the variance of $X_*$ is upper bounded by the Poincaré inequality:

$$\mathrm{Var} X_* = \sum_{i=1}^d \left( \int |x^i|^2 m_*(dx) - \left( \int x^i m_*(dx) \right)^2 \right) \\ \le C_{\mathrm{LS}} \sum_{i=1}^d \int |\nabla x^i|^2 m_*(dx) = C_{\mathrm{LS}} d.$$

We then conclude by combining the three inequalities above. □

## 6.3 Construction of another measure

In this subsection, we construct another measure in order to exploit the convexity of $\Phi$, used for the proof of the first and second claims of the theorem. Readers only interested in the last claim of the theorem can now directly go to the next subsection.

Let $\mu = \mu^\lambda$ be the probability measure on $\mathbb{R}^d \times \mathbb{R}^D$ characterized by the following formula:

$$\langle f, \mu \rangle = \int f(x, y)\mu(dx dy) \\ = \mathbb{E} \left[ \int f(x, \langle l, m \rangle) m(dx) \right] = \mathbb{E} \left[ \int f(x, Y) m(dx) \right] \quad (23)$$

for all bounded and measurable $f: \mathbb{R}^d \times \mathbb{R}^D \to \mathbb{R}$. By taking $f$ depending only on the $y$ variable, we first realize that the second marginals of $\rho$ and $\mu$ agree, that is,

$$\rho^2 = \mu^2.$$

In addition, we show the following important properties of $\mu$.
**Proposition 9.** Under the setting of the theorem, for every $C^2$ differentiable $\Psi : \mathbb{R}^D \to \mathbb{R}$ with bounded Hessian, we have

$$ \iint \nabla \Psi(y) \cdot \ell(x)(\mu - \rho)(dx \, dy) = 0. $$

In particular, the respective first marginals $\mu^1$, $\rho^1$ of $\mu$, $\rho$ satisfy

$$ \int \ell(x)(\mu^1 - \rho^1)(dx) = 0. $$

Moreover, denoting by $\mu^{1|2}(\cdot|\cdot) : \mathcal{B}(\mathbb{R}^d) \times \mathbb{R}^D \to \mathbb{R}$ the conditional measure of $\mu$ given its second variable, we have for $\rho^2$-almost all $y \in \mathbb{R}^D$,

$$ \langle \ell, \mu^{1|2}(\cdot|y) \rangle = y. $$

*Proof of Proposition 9.* Consider the functional

$$ \phi(m) = \int f(x', \langle k, m \rangle)m(dx'), $$

where $f \in C_b^1(\mathbb{R}^d \times \mathbb{R}^D; \mathbb{R})$ and $k \in C_b(\mathbb{R}^d; \mathbb{R}^D)$. Then its linear functional derivative reads

$$ \frac{\delta \phi}{\delta m}(m, x') = f(x', \langle k, m \rangle) + \int \nabla_y f(x'', \langle k, m \rangle) \cdot k(x')m(dx''), $$

so $\phi$ belongs to the $C_b^1$ class. Then, applying Proposition 5 to the functional $\phi$, we get

$$ \begin{aligned} 0 &= \mathrm{E}\left[\int \frac{\delta\phi}{\delta m}(m, x')(\delta_{X}-m)(dx')\right] \\ &= \mathrm{E}[f(X, m)] - \mathrm{E}\left[\int f(x', \langle k, m \rangle)m(dx')\right] \\ &\quad + \mathrm{E}\left[\int \nabla_y f(x', \langle k, m \rangle) \cdot (k(X) - \langle k, m \rangle)m(dx')\right] \end{aligned} $$

By approximation, the equality above holds for $k = \ell$ and for all $C^1$-continuous $f$ with the following growth bounds:

$$ \begin{aligned} |f(x, y)| &\le M(1 + |x|^4 + |y|^2), \\ |\nabla_y f(x, y)| &\le M(1 + |x|^2 + |y|), \end{aligned} $$

that is to say, we have

$$ \langle f, \rho - \mu \rangle + \mathrm{E}\left[\int \nabla_y f(x', Y)(\ell(X) - Y)m(dx')\right] = 0, $$

where, as before, $Y = \langle \ell, m \rangle$. Specializing to $f(x, y) = \nabla\Psi(y) \cdot \ell(x)$, we obtain

$$ \begin{aligned} \langle f, \mu - \rho \rangle &= \mathrm{E}\left[\int \ell(x)^{\top} \nabla^2 \Psi(Y)(\ell(X) - Y)m(dx)\right] \\ &= \mathrm{E}[Y^{\top} \nabla^2 \Psi(Y)(\ell(X) - Y)] = 0, \end{aligned} $$
where the last equality is due to Lemma 6, as for $\varphi(y) := \nabla\Psi(y) \cdot y - \Psi(y)$, we have $\nabla\varphi(y) = \nabla^2\Psi(y)y$. So the first claim is proved. Taking $\Psi(y) = y'$, for $\nu = 1, \dots, D$, yields the second claim.

For the last claim, we take $f(x, y) = \ell(x)g(y)$ for $g: \mathbb{R}^D \to \mathbb{R}$ of linear growth in the defining equation (23) of $\mu$. Then, we get

$$ \int g(y) \left( \int \ell(x) \mu^{1/2}(dx|y) \right) \mu^2(dy) = \iint f(x, y) \mu(dx dy) = \mathbb{E}[Yg(Y)] \\ = \int g(y) y \mu^2(dy). $$

The desired property follows from the arbitrariness of $g$. $\square$

## 6.4 Proving the theorem

Having established the entropy estimate and constructed the auxiliary measure, we finally move to the central part of the proof, which consists of six steps. The aim of the first five steps is to show the first and the second claims of the theorem, and in the last step we prove the last claim.

**Step 1: Control of the symmetrized entropy.** We aim at controlling the symmetrized entropy

$$ \int (H(\hat{m}_y|m_*) + H(m_*|\hat{m}_y)) \rho^2(dy) $$

in this step. First observe

$$ \int (H(\hat{m}_y|m_*) + H(m_*|\hat{m}_y)) \rho^2(dy) \\ = 2 \iint (V(y_*, x) - V(y, x))(\hat{m}_y(dx) - m_*(dx)) \rho^2(dy). \quad (24) $$

In order to control the right hand side above, we turn to the probability measure $\mu$ introduced in (23). Recall that $m_*$ is the invariant measure of the McKean-Vlasov (1), and $y_* := \langle \ell, m_* \rangle$. The convexity of $\Phi$ implies the convexity of $F$ as a functional, and as a result, for $\rho^2$-almost all $y \in \mathbb{R}^D$, we have the tangent inequalities

$$ \int V(\langle \ell, \mu^{1/2}(\cdot|y) \rangle, x) (\mu^{1/2}(dx|y) - m_*(dx)) \\ \geq F(\mu^{1/2}(\cdot|y)) - F(m_*) \\ \geq \int V(y_*, x) (\mu^{1/2}(dx|y) - m_*(dx)). \quad (25) $$

Thanks to the last claim of Proposition 9, the leftmost term satisfies, for $\rho^2$-almost all $y \in \mathbb{R}^D$,

$$ \int V(\langle \ell, \mu^{1/2}(\cdot|y) \rangle, x) (\mu^{1/2}(dx|y) - m_*(dx)) = \int V(y, x) (\mu^{1/2}(dx|y) - m_*(dx)). $$

Hence, integrating the tangent inequalities (25) above by $\rho^2$, we get

$$ \iint V(y, x) (\mu^{1/2}(dx|y) - m_*(dx)) \rho^2(dy) \\ \geq \iint V(y_*, x) (\mu^{1/2}(dx|y) - m_*(dx)) \rho^2(dy). \quad (26) $$
Using $\mu^2 = \rho^2$ and applying Proposition 9 to $V(y, x) = \nabla\Phi(y) \cdot \ell(x)$, we know that the left hand side of (26) satisfies

$$
\begin{align*}
\iint V(y, x)(\mu^{1/2}(dx|y) - m_*(dx))\rho^2(dy) \\
&= \iint V(y, x)\mu(dx\,dy) - \iint V(y, x)m_*(dx)\rho^2(dy) \\
&= \iint V(y, x)\rho(dx\,dy) - \iint V(y, x)m_*(dx)\rho^2(dy) \\
&= \iint V(y, x)(\rho^{1/2}(dx|y) - m_*(dx))\rho^2(dy).
\end{align*}
$$

The right hand side of (26) satisfies

$$
\begin{align*}
\iint V(y_*, x)(\mu^{1/2}(dx|y) - m_*(dx))\rho^2(dy) \\
&= \iint V(y_*, x)\mu(dx\,dy) - \iint V(y_*, x)m_*(dx) \\
&= \nabla\Phi(y_*) \cdot \int \ell(x)(\mu^1(dx) - m_*(dx)) \\
&= \nabla\Phi(y_*) \cdot \int \ell(x)(\rho^1(dx) - m_*(dx)) \\
&= \iint V(y_*, x)(\rho^{1/2}(dx|y) - m_*(dx))\rho^2(dy),
\end{align*}
$$

where the third equality is due to the last claim of Proposition 9. Thus, we have
derived

$$
\iint (V(y, x) - V(y_*, x))(\rho^{1/2}(dx|y) - m_*(dx))\rho^2(dy) \ge 0. \quad (27)
$$

Therefore, to dominate the right hand side of (24), it remains to control the fol-
lowing term. Using the Kantorovich duality, we get

$$
\begin{align*}
& \left| \iint (V(y, x) - V(y_*, x))(\rho^{1/2}(dx|y) - \hat{m}_y(dx))\rho^2(dy) \right| \\
&= \left| \int (\nabla\Phi(y) - \nabla\Phi(y_*)) \cdot \langle \ell, \rho^{1/2}(\cdot|y) - \hat{m}_y \rangle \rho^2(dy) \right| \\
&= \left| \iint_0^1 (y - y_*)^\top \nabla^2\Phi((1-t)y + ty_*) \langle \ell, \rho^{1/2}(\cdot|y) - \hat{m}_y \rangle dt \rho^2(dy) \right| \\
&\le \left( \iint_0^1 (y - y_*)^\top \nabla^2\Phi((1-t)y + ty_*) (y - y_*) dt \rho^2(dy) \right)^{1/2} \\
&\qquad \times \left( \iint_0^1 \left| \langle \nabla^2\Phi((1-t)y + ty_*)^{1/2} \ell, \rho^{1/2}(\cdot|y) - \hat{m}_y \rangle \right|^2 dt \rho^2(dy) \right)^{1/2} \\
&\le \left( \iint_0^1 (y - y_*)^\top \nabla^2\Phi((1-t)y + ty_*) (y - y_*) dt \rho^2(dy) \right)^{1/2} \\
&\qquad \times \sqrt{M_2} \left( \int W_1^2(\rho^{1/2}(\cdot|y), \hat{m}_y) \rho^2(dy) \right)^{1/2} \\
&:= \sqrt{M_2 v(\rho^2)} \left( \int W_1^2(\rho^{1/2}(\cdot|y), \hat{m}_y) \rho^2(dy) \right)^{1/2},
\end{align*}
$$
where, by Assumption 4, $|\nabla^2\Phi(y')^{1/2}\nabla\ell(x)| \le \sqrt{M_2}$ for all $x \in \mathbb{R}^d$, $y' \in \mathbb{R}^D$, and $v(\rho^2)$ is the quantity to be controlled in the first claim of the theorem. The uniform LSI for $(\hat{m}_y)_{y \in \mathcal{K}}$ implies a uniform Talagrand's transport inequality, from which we obtain

$$ \begin{aligned} \int W_1^2(\rho^{1/2}(\cdot|y), \hat{m}_y)\rho^2(dy) &\le \int W_2^2(\rho^{1/2}(\cdot|y), \hat{m}_y)\rho^2(dy) \\ &\le C_{LS} \int H(\rho^{1/2}(\cdot|y)|\hat{m}_y)\rho^2(dy) = C_{LS}H(\rho|\hat{\rho}), \end{aligned} $$

as $\rho^2$ is supported on $\mathcal{K}$. Combining the two inequalities above with (24) and (27), we obtain

$$ \int (H(\hat{m}_y|m_*) + H(m_*|\hat{m}_y))\rho^2(dy) \le 2\sqrt{M_2C_{LS}v(\rho^2)H(\rho|\hat{\rho})}. \quad (28) $$

**Step 2: Control of the conditional Wasserstein distance.** Now, using again the Talagrand's transport inequality for $\hat{m}_y$ and $m_*$ (note that $m_* = \hat{m}_{y_*}$ for $y_* = \langle \ell, m_* \rangle$), we get for $\rho^2$-almost all $y \in \mathbb{R}^D$,

$$ W_2^2(\hat{m}_y, m_*) \le \frac{C_{LS}}{2} (H(\hat{m}_y|m_*) + H(m_*|\hat{m}_y)), $$

while the triangle inequality and the transport inequality imply

$$ \begin{aligned} W_2^2(\rho^{1/2}(\cdot|y), m_*) &\le 2(W_2^2(\rho^{1/2}(\cdot|y), \hat{m}_y) + W_2^2(\hat{m}_y, m_*)) \\ &\le 2C_{LS}H(\rho^{1/2}(\cdot|y), \hat{m}_y) + 2W_2^2(\hat{m}_y, m_*). \end{aligned} $$

So, combining the three inequalities above and integrating with $\rho^2$, we find

$$ \int W_2^2(\rho^{1/2}(\cdot|y), m_*)\rho^2(dy) \le 2C_{LS}H(\rho|\hat{\rho}) + 2\sqrt{M_2C_{LS}^3v(\rho^2)H(\rho|\hat{\rho})}. \quad (29) $$

**Step 3: Control of $v(\rho^2)$ by $H(\rho|\hat{\rho})$.** Applying Proposition 9 to the function $\Psi(y) = \Phi(y) - \nabla\Phi(y_*) \cdot y$, where, as we recall, $y_* = \langle \ell, m_* \rangle$, we get

$$ \begin{aligned} 0 &= \iint \nabla\Psi(y) \cdot \ell(x)(\mu - \rho)(dx dy) \\ &= \iint (\nabla\Phi(y) - \nabla\Phi(y_*)) \cdot \ell(x)(\mu - \rho)(dx dy) \\ &= \iint_0^1 (y - y_*)^\top \nabla^2\Phi((1-t)y + ty_*) \langle \ell, \mu^{1/2}(\cdot|y) - \rho^{1/2}(\cdot|y) \rangle dt \rho^2(dy) \\ &= \iint_0^1 (y - y_*)^\top \nabla^2\Phi((1-t)y + ty_*) (y - \langle \ell, \rho^{1/2}(\cdot|y) \rangle) dt \rho^2(dy), \end{aligned} $$

where the last equality is due to the last claim of Proposition 9. In other words, we have

$$ \begin{aligned} &\iint_0^1 (y - y_*)^\top \nabla^2\Phi((1-t)y + ty_*) (y - y_*) dt \rho^2(dy) \\ &= \iint_0^1 (y - y_*)^\top \nabla^2\Phi((1-t)y + ty_*) \langle \ell, \rho^{1/2}(\cdot|y) - m_* \rangle dt \rho^2(dy), \end{aligned} $$
and this implies, by Cauchy-Schwarz,

$$
\begin{align*}
\iint_0^1 (y - y_*)^\top \nabla^2 \Phi((1-t)y + ty_*) (y - y_*) \, dt \, \rho^2(dy) & \\
&\le \iint_0^1 \langle \ell, \rho^{1/2}(\cdot|y) - m_* \rangle^\top \nabla^2 \Phi((1-t)y + ty_*) \langle \ell, \rho^{1/2}(\cdot|y) - m_* \rangle \, dt \, \rho^2(dy) \\
&= \iint_0^1 \left| \langle \nabla^2 \Phi((1-t)y + ty_*)^{1/2} \ell, \rho^{1/2}(\cdot|y) - m_* \rangle \right|^2 \, dt \, \rho^2(dy).
\end{align*}
$$

As $|\nabla^2 \Phi(y')^{1/2} \nabla \ell(x)| \le \sqrt{M_2}$ for all $y' \in \mathbb{R}^D$ and $x \in \mathbb{R}^d$, we have, by the Kantorovich duality,

$$
v(\rho^2) = \iint_0^1 (y - y_*)^\top \nabla^2 \Phi((1-t)y + ty_*) (y - y_*) \, dt \, \rho^2(dy) \\
\le M_2 \int W_1^2(\rho^{1/2}(\cdot|y), m_*) \rho^2(dy). \quad (30)
$$

Thus, using the fact that $W_1 \le W_2$ and the inequality (29), we obtain

$$
v(\rho^2) \le 2M_2 C_{LS} H(\rho|\hat{\rho}) + 2\sqrt{M_2^3 C_{LS}^3 v(\rho^2) H(\rho|\hat{\rho})}.
$$

Introduce the “adimensionalized” variable

$$
v = \frac{v(\rho^2)}{4M_2^3 C_{LS}^3 H(\rho|\hat{\rho})}.
$$

Then the inequality above reads

$$
v \le \frac{1}{2M_2^2 C_{LS}^2} + \sqrt{v} \le \frac{1}{2M_2^2 C_{LS}^2} + \frac{1}{2} + \frac{1}{2}v.
$$

Hence, we get

$$
v(\rho^2) \le 4M_2 C_{LS} (M_2^2 C_{LS}^2 + 1) H(\rho|\hat{\rho}). \quad (31)
$$

Step 4: Control of Wasserstein and TV distances by H(ρ|ρ̂). By inserting (31) into (29), and noting, by the definition of the Wasserstein distance,

$$
W_2^2(\rho^1, m_*) \le \int W_2^2(\rho^{1/2}(\cdot|y), m_*) \rho^2(dy),
$$

we get

$$
W_2^2(\rho^1, m_*) \le (2C_{LS} + 4M_2 C_{LS}^2 (M_2^2 C_{LS}^2 + 1)^{1/2}) H(\rho|\hat{\rho}). \quad (32)
$$

For the total variation distance, we observe that the Csiszár-Kullback-Pinsker inequality implies

$$
\begin{align*}
\int \|\rho^{1/2}(\cdot|y) - \hat{m}_y\|_{\text{TV}}^2 \rho^2(dy) &\le 2H(\rho|\hat{\rho}), \\
\int \|\hat{m}_y - m_*\|_{\text{TV}}^2 \rho^2(dy) &\le \int (H(\hat{m}_y|m_*) + H(m_*|\hat{m}_y)) \rho^2(dy).
\end{align*}
$$
By the triangle and Jensen's inequalities, we get

$$
\begin{align*}
\|\rho^1 - m_*\|_{\text{TV}}^2 &\le \int \|\rho^{1/2}(\cdot|y) - m_*\|_{\text{TV}}^2 \rho^2(dy) \\
&\le 2 \int \left( \|\rho^{1/2}(\cdot|y) - \hat{m}_y\|_{\text{TV}}^2 + \|\hat{m}_y - m_*\|_{\text{TV}}^2 \right) \rho^2(dy) \\
&\le 4H(\rho|\hat{\rho}) + 2 \int (H(\hat{m}_y|m_*) + H(m_*|\hat{m}_y))\rho^2(dy).
\end{align*}
$$

Then, by inserting (31) into (28), we get

$$
\|\rho^1 - m_*\|_{\text{TV}}^2 \le (4 + 8M_2C_{\text{LS}}(M_2^2C_{\text{LS}}^2 + 1)^{1/2})H(\rho|\hat{\rho}). \quad (33)
$$

Step 5: Control of $H(\rho|\hat{\rho})$ and conclusion for the convex case. In the case where

$$
M_1 := \sup_{x \in \mathbb{R}^d, y \in \mathbb{R}^D} \ell(x)^\top \nabla^2 \Phi(y) \ell(x) < \infty,
$$

by Proposition 7 and (22) in Proposition 8, we immediately get

$$
H(\rho|\hat{\rho}) \le \frac{C_{\text{LS}}}{2} (D + 2M_1)\lambda = H'.
$$

In the general case where $M_1$ is not necessarily finite, Proposition 7 and (21) in
Proposition 8 yield

$$
H(\rho|\hat{\rho}) \le \frac{C_{LS}}{2} (D + 24M_2 C_{LS} d) \lambda + 24M_2 C_{LS} W_2^2(\rho^1, m_*) \lambda.
$$

Together with the upper bound (32) of $W_2^2(\rho^1, m_*)$, we get

$$
H(\rho|\hat{\rho}) \le \frac{C_{\text{LS}}(D + 24M_2 C_{\text{LS}}d)\lambda}{2 - 96M_2 C_{\text{LS}}^2 (1 + 2M_2 C_{\text{LS}}(M_2^2 C_{\text{LS}}^2 + 1)^{1/2})\lambda} = H,
$$

for

$$
\lambda < \frac{1}{48 M_2 C_{\text{LS}}^2 \left(1 + 2 M_2 C_{\text{LS}} \left(M_2^2 C_{\text{LS}}^2 + 1\right)^{1/2}\right)} = \lambda_0.
$$

We obtain the desired estimates on $v(\rho^2)$, Wasserstein and TV distances, by insert-
ing the upper bounds of $H(\rho|\hat{\rho})$ for the respective cases into (31), (32), (33).

Now we work with a concave $\Phi$ and prove the last claim of the theorem.

Step 6: Case of concave $\Phi$. Observe first that the mapping $y \mapsto \nabla^2\Phi(y) \cdot \hat{y}$ is a gradient:

$$
\begin{align*}
\nabla^2 \Phi(y) \cdot \hat{y} &= \nabla^2 \Phi(y) \cdot \frac{\int \ell(x) \exp(-2\nabla\Phi(y) \cdot \ell(x)) \,dx}{\int \exp(-2\nabla\Phi(y) \cdot \ell(x)) \,dx} \\
&= -\frac{1}{2} \nabla_y \left( \log \int \exp(-2\nabla\Phi(y) \cdot \ell(x)) \,dx \right).
\end{align*}
$$

This identity is analogous to the fact in thermodynamics that when we derive the
free energy with respect to a variable, we get the statistical average of its response
variable. Moreover,

$$
\nabla^2 \Phi(y) \cdot y = \nabla_y (\nabla \Phi(y) \cdot y - \Phi(y)).
$$
Thus, taking the test function

$$ \varphi(y) = -\frac{1}{2} \log \int \exp(-2\nabla\Phi(y) \cdot \ell(x)) \, dx - \nabla\Phi(y) \cdot y + \Phi(y) $$

in Lemma 6, we get

$$ \iint (\ell(x) - y)^{\top} \nabla^2 \Phi(y) (\hat{y} - y) \rho(dx dy) = 0. $$

Consequently,

$$ \begin{aligned} & -\iint (\hat{y} - y)^{\top} \nabla^2 \Phi(y) (\hat{y} - y) \rho(dx dy) \\ &= -\iint (\hat{y} - \ell(x))^{\top} \nabla^2 \Phi(y) (\hat{y} - y) \rho(dx dy) \\ &= \int \langle (-\nabla^2 \Phi(y))^{1/2} \ell, \hat{m}_y - \rho^{1/2}(\cdot|y) \rangle \cdot (-\nabla^2 \Phi(y))^{1/2} (\hat{y} - y) \rho^2(dy) \\ &\leq \left( -\int (\hat{y} - y)^{\top} \nabla^2 \Phi(y) (\hat{y} - y) \rho^2(dy) \right)^{1/2} \\ &\qquad \times \left( \int \left| \langle (-\nabla^2 \Phi(y))^{1/2} \ell, \hat{m}_y - \rho^{1/2}(\cdot|y) \rangle \right|^2 \rho^2(dy) \right)^{1/2}, \end{aligned} $$

which implies that

$$ -\int (\hat{y} - y)^{\top} \nabla^2 \Phi(y) (\hat{y} - y) \rho^2(dy) \leq \int \left| \langle (-\nabla^2 \Phi(y))^{1/2} \ell, \hat{m}_y - \rho^{1/2}(\cdot|y) \rangle \right|^2 \rho^2(dy). $$

The term on the right satisfies

$$ \begin{aligned} \int \left| \langle (-\nabla^2 \Phi(y))^{1/2} \ell, \hat{m}_y - \rho^{1/2}(\cdot|y) \rangle \right|^2 \rho^2(dy) &\leq M_2 \int W_1^2(\hat{m}_y, \rho^{1/2}(\cdot|y)) \rho^2(dy) \\ &\leq M_2 C_{LS} \int H(\rho^{1/2}(\cdot|y) | \hat{m}_y) \rho^2(dy) \end{aligned} $$

by the uniform LSI for $(\hat{m}_y)_{y \in \mathcal{K}}$. The entropy estimate in Proposition 7 gives

$$ \int H(\rho^{1/2}(\cdot|y) | \hat{m}_y) \rho^2(dy) \leq \frac{C_{LS}}{2}(D + 2I)\lambda \leq \frac{C_{LS}D\lambda}{2}, $$

as in the case of concave $\Phi$, the term $I \leq 0$ by its definition (16). $\square$

<figure><img src="image_5.png" /></figure>

## A A lemma on Wasserstein duality

**Lemma 10.** Let $f: \mathbb{R}^{d_1} \times \mathbb{R}^{d_2} \to \mathbb{R}$ be a $C^2$ function such that the Euclidean operator norm of $\nabla_x \nabla_y f(x, y)$ satisfies

$$ |\nabla_x \nabla_y f(x, y)| \leq M $$

for all $x \in \mathbb{R}^{d_1}$ and $y \in \mathbb{R}^{d_2}$. Then, for all $\mu, \mu' \in \mathcal{P}_2(\mathbb{R}^{d_1})$ and $\nu, \nu' \in \mathcal{P}_2(\mathbb{R}^{d_2})$, we have

$$ \left| \iint f(x, y)(\mu - \mu')(dx)(\nu - \nu')(dy) \right| \leq MW_2(\mu, \mu')W_2(\nu, \nu'). $$
*Proof of Lemma 10.* Let $\pi \in \Pi(\mu, \mu')$ be the $W_2$-optimal transport plan between $\mu$ and $\mu'$, and $\pi' \in \Pi(\nu, \nu')$ be that between $\nu$ and $\nu'$. Construct the random variable $((X, X'), (Y, Y'))$ distributed as $\pi \otimes \pi'$. Then, we have

$$ \iint f(x, y)(\mu - \mu')(dx)(\nu - \nu')(dy) \\
= \mathbb{E}[f(X, Y) - f(X', Y) - f(X, Y') + f(X', Y')] \\
= \mathbb{E} \left[ \iint_{[0,1]^2} (X - X')^\top \nabla_x \nabla_y f((1-t)X + tX', (1-s)Y + sY')(Y - Y') dt ds \right]. $$

Therefore, taking absolute values, we get

$$ \left| \iint f(x, y)(\mu - \mu')(dx)(\nu - \nu')(dy) \right| \\
\leq M \mathbb{E}[|X - X'| |Y - Y'|] \\
\leq M \mathbb{E}[|X - X'|^2]^{1/2} \mathbb{E}[|Y - Y'|^2]^{1/2} \\
= MW_2(\mu, \mu')W_2(\nu, \nu'), $$

which concludes the proof.

<figure><img src="image_6.png" /></figure>

# B Algorithm

**Algorithm 1:** Training two-layer neural network by self-interacting diffusion

**Input** : Activation func. $\varphi$, truncation func. $\tau$, dataset $\{z_k, y_k\}_{k=1}^K$, volatility $\sigma$, regularization constant $\gamma$, initial distribution $m_0$, time step $\delta t$, time horizon $T$, non-increasing piecewise constant func. $\lambda : [0, T] \to (0, \infty)$

**Output:** $Y_T$

generate $X_0 = (C_0, A_0, B_0) \sim m_0$; $Y_0 = l(X_0)$;

**for** $t = 0, \delta t, 2\delta t, \dots, T - \delta t$ **do**

generate i.i.d. $N_t \sim N(0, 1)$;
$X_{t+\delta t} \leftarrow X_t - \left(\frac{1}{K} \sum_{k=1}^K Y_t^k \nabla f(X_t, z_k) + \sigma^2 \gamma X_t\right) \delta t + \sigma \sqrt{\delta t} N_t$;
**if** $\lambda(t) \equiv \lambda_*$ **is constant function then**
$Y_{t+\delta t} \leftarrow (1 - \lambda_* \delta t) Y_t + \lambda_* \delta t \ell(X_{t+\delta t})$;
// Update corresponding to fixed value of $\lambda$
**else**
$Y_{t+\delta t} \leftarrow (1 - \lambda(t) \delta t) Y_t + \lambda(t) \delta t \ell(X_{t+\delta t})$;
// Update corresponding to annealing $\lambda(t) = \lambda_a(t)$

**Acknowledgements.** The authors thank Fan Chen and Xiaozhen Wang for their assistance with the numerical implementation, and Yupeng Bai for identifying an error in Theorem 1 in a previous version of the paper.
**Funding.** The first author's research is supported by the National Natural Science Foundation of China (No. 12222103), and by the National Key R&D Program of China (No. 2022ZD0116401). The second author's research is supported by Finance For Energy Market Research Centre and France 2030 funding (ANR-21-EXES-0003). The third author's research is supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 945332.

References

[1] Luigi Ambrosio, Nicola Gigli and Giuseppe Savaré. *Gradient flows in metric spaces and in the space of probability measures*. Birkhäuser, Basel, 2nd edition, 2008. ISBN: 978-3-7643-8721-1.

[2] Dominique Bakry and Michel Émery. Diffusions hypercontractives. French. In *Sémin. de Probab. XIX*. Volume 1123, Lect. Notes Math. Springer-Verlag, Berlin etc., 1985.

[3] Michel Benaïm, Michel Ledoux and Olivier Raimond. Self-interacting diffusions. *Probab. Theory Relat. Fields*, 122(1):1–41, 2002. ISSN: 0178-8051. DOI: 10.1007/s004400100161.

[4] Michel Benaïm and Olivier Raimond. Self-interacting diffusions. II: Convergence in law. *Ann. Inst. Henri Poincaré, Probab. Stat.*, 39(6):1043–1055, 2003. ISSN: 0246-0203.

[5] Michel Benaïm and Olivier Raimond. Self-interacting diffusions. III: Symmetric interactions. *Ann. Probab.*, 33(5):1716–1759, 2005. ISSN: 0091-1798. DOI: 10.1214/009117905000000251.

[6] Michel Benaïm and Olivier Raimond. Self-interacting diffusions. IV: Rate of convergence. *Electron. J. Probab.*, 16, 66, 2011. ISSN: 1083-6489. DOI: 10.1214/EJP.v16-948.

[7] Dario Benedetto, Emanuele Caglioti, José A. Carrillo and Mario Pulvirenti. A non-Maxwellian steady distribution for one-dimensional granular media. *J. Stat. Phys.*, 91(5-6):979–990, 1998. ISSN: 0022-4715. DOI: 10.1023/A:1023032000560.

[8] George D. Birkhoff. Proof of the ergodic theorem. *Proc. Natl. Acad. Sci. USA*, 17:656–660, 1931. ISSN: 0027-8424. DOI: 10.1073/pnas.17.12.656.

[9] Vladimir I. Bogachev, Nicolai V. Krylov, Michael Röckner and Stanislav V. Shaposhnikov. *Fokker-Planck-Kolmogorov equations*, volume 207 of *Math. Surv. Monogr*. American Mathematical Society, Providence, RI, 2015. ISBN: 978-1-4704-2558-6. DOI: 10.1090/surv/207.

[10] François Bolley, Ivan Gentil and Arnaud Guillin. Uniform convergence to equilibrium for granular media. *Arch. Ration. Mech. Anal.*, 208(2):429–445, 2013. ISSN: 0003-9527. DOI: 10.1007/s00205-012-0599-z.

[11] José A. Carrillo, Robert J. McCann and Cédric Villani. Kinetic equilibration rates for granular media and related equations: entropy dissipation and mass transportation estimates. *Rev. Mat. Iberoam.*, 19(3):971–1018, 2003. ISSN: 0213-2230. DOI: 10.4171/RMI/376.
[12] Fan Chen, Yiqing Lin, Zhenjie Ren and Songbo Wang. Uniform-in-time propagation of chaos for kinetic mean field Langevin dynamics. *Electron. J. Probab.*, 29, 17, 2024. doi: 10.1214/24-EJP1079.

[13] Fan Chen, Zhenjie Ren and Songbo Wang. Entropic fictitious play for mean field optimization problem. *J. Mach. Learn. Res.*, 24, 211, 2023. URL: http://jmlr.org/papers/v24/22-0411.html.

[14] Fan Chen, Zhenjie Ren and Songbo Wang. Uniform-in-time propagation of chaos for mean field Langevin dynamics. *Ann. Inst. Henri Poincaré, Probab. Stat.*, 2024. arXiv: 2212.03050 [math.PR]. Forthcoming.

[15] Lénaïc Chizat. Mean-field Langevin dynamics: Exponential convergence and annealing. *Trans. Mach. Learn. Res.*, 2022. ISSN: 2835-8856. URL: https://openreview.net/forum?id=BDqzLH1gEm.

[16] Lénaïc Chizat and Francis Bach. On the global convergence of gradient descent for over-parameterized models using optimal transport. In *Adv. Neural Inf. Process. Syst*. Volume 31, pages 3036–3046. Curran Associates, Inc., 2018. URL: https://proceedings.neurips.cc/paper_files/paper/2018/file/a1afc58c6ca9540d057299ec3016d726-Paper.pdf.

[17] Giovanni Conforti, Anna Kazeykina and Zhenjie Ren. Game on random environment, mean-field Langevin system, and neural networks. *Math. Oper. Res.*, 48(1):78–99, 2023. doi: 10.1287/moor.2022.1252.

[18] Michael Cranston and Yves Le Jan. Self attracting diffusions: Two case studies. *Math. Ann.*, 303(1):87–93, 1995. ISSN: 0025-5831. DOI: 10.1007/BF01460980.

[19] François Delarue and Alvin Tse. Uniform in time weak propagation of chaos on the torus. *Ann. Inst. Henri Poincaré Probab. Stat.*, 61(2):1021–1074, 2025. ISSN: 0246-0203,1778-7017. DOI: 10.1214/23-aihp1451.

[20] Luc Devroye, Abbas Mehrabian and Tommy Reddad. The total variation distance between high-dimensional gaussians with the same mean. 2023. arXiv: 1810.08693 [math.ST].

[21] Kai Du, Yifan Jiang and Jinfeng Li. Empirical approximation to invariant measures for McKean-Vlasov processes: mean-field interaction vs self-interaction. *Bernoulli*, 29(3):2492–2518, 2023. ISSN: 1350-7265. DOI: 10.3150/22-BEJ1550.

[22] Kai Du, Yifan Jiang and Xiaochen Li. Sequential propagation of chaos. 2023. arXiv: 2301.09913 [math.PR].

[23] Alain Durmus, Andreas Eberle, Arnaud Guillin and Raphael Zimmer. An elementary approach to uniform in time propagation of chaos. *Proc. Am. Math. Soc.*, 148(12):5387–5398, 2020. ISSN: 0002-9939. DOI: 10.1090/proc/14612.

[24] Alain Durmus and Éric Moulines. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. *Ann. Appl. Probab.*, 27(3):1551–1587, 2017. ISSN: 1050-5164. DOI: 10.1214/16-AAP1238.

[25] Andreas Eberle. Reflection couplings and contraction rates for diffusions. *Probab. Theory Relat. Fields*, 166(3-4):851–886, 2016. ISSN: 0178-8051. DOI: 10.1007/s00440-015-0673-1.
[26] Andreas Eberle, Arnaud Guillin and Raphael Zimmer. Couplings and quan-
titative contraction rates for Langevin dynamics. *Ann. Probab.*, 47(4):1982-
2010, 2019. ISSN: 0091-1798. DOI: 10.1214/18-AOP1299.

[27] Wendell H. Fleming and Michel Viot. Some measure-valued Markov processes in population genetics theory. *Indiana Univ. Math. J.*, 28:817-843, 1979. ISSN: 0022-2518. DOI: 10.1512/iumj.1979.28.28058.

[28] Arnaud Guillin and Pierre Monmarché. Uniform long-time and propagation of chaos estimates for mean field kinetic particles in non-convex landscapes. *J. Stat. Phys.*, 185(2), 15, 2021. ISSN: 0022-4715. DOI: 10.1007/s10955-021-02839-6.

[29] Richard Holley and Daniel Stroock. Logarithmic Sobolev inequalities and stochastic Ising models. *J. Stat. Phys.*, 46(5-6):1159-1194, 1987. ISSN: 0022-4715. DOI: 10.1007/BF01011161.

[30] Kaitong Hu, Zhenjie Ren, David Šiška and Łukasz Szpruch. Mean-field Langevin dynamics and energy landscape of neural networks. *Ann. Inst. Henri Poincaré, Probab. Stat.*, 57(4):2043-2065, 2021. ISSN: 0246-0203. DOI: 10.1214/20-AIHP1140.

[31] Evelyn F. Keller and Lee A. Segel. Model for chemotaxis. *J. Theor. Biol.*, 30(2):225-234, 1971. ISSN: 0022-5193. DOI: 10.1016/0022-5193(71)90050-6.

[32] Victor Kleptsyn and Aline Kurtzmann. Ergodicity of self-attracting motion. *Electron. J. Probab.*, 17, 50, 2012. ISSN: 1083-6489. DOI: 10.1214/EJP.v17-2121.

[33] Yoshiki Kuramoto. Rhythms and turbulence in populations of chemical oscillators. *Phys. A, Stat. Mech. Appl.*, 106(1):128-143, 1981. ISSN: 0378-4371. DOI: 10.1016/0378-4371(81)90214-4.

[34] Daniel Lacker and Luc Le Flem. Sharp uniform-in-time propagation of chaos. *Probab. Theory Relat. Fields*, 187(1-2):443-480, 2023. ISSN: 0178-8051. DOI: 10.1007/s00440-023-01192-x.

[35] Michel Ledoux and Michel Talagrand. *Probability in Banach spaces. Isoperimetry and processes*, volume 23 of *Ergeb. Math. Grenzgeb., 3. Folge*. Springer-Verlag, Berlin etc., 1991. ISBN: 3-540-52013-9.

[36] Nicolas Martzel and Claude Aslangul. Mean-field treatment of the many-body Fokker-Planck equation. *J. Phys. A, Math. Gen.*, 34(50):11225-11240, 2001. ISSN: 0305-4470. DOI: 10.1088/0305-4470/34/50/305.

[37] Song Mei, Andrea Montanari and Phan-Minh Nguyen. A mean field view of the landscape of two-layer neural networks. *Proc. Natl. Acad. Sci. USA*, 115(33):e7665-e7671, 2018. ISSN: 0027-8424. DOI: 10.1073/pnas.1806579115.

[38] Atsushi Nitanda, Denny Wu and Taiji Suzuki. Convex analysis of the mean field Langevin dynamics. In *Proceedings of the 25th International Conference on Artificial Intelligence and Statistics*, volume 151 of *Proceedings of Machine Learning Research*, pages 9741-9757. PMLR, 2022. URL: https://proceedings.mlr.press/v151/nitanda22a.html.
[39] Atsushi Nitanda, Denny Wu and Taiji Suzuki. Particle dual averaging: optimization of mean field neural network with global convergence rate analysis. *J. Stat. Mech. Theory Exp.*, 2022(11), 114010, 2022. ISSN: 1742-5468. DOI: 10.1088/1742-5468/ac98a8.

[40] Olivier Raimond. Self-attracting diffusions: Case of the constant interaction. *Probab. Theory Relat. Fields*, 107(2):177-196, 1997. ISSN: 0178-8051. DOI: 10.1007/s004400050082.

[41] Gareth O. Roberts and Jeffrey S. Rosenthal. Optimal scaling of discrete approximations to Langevin diffusions. *J. R. Stat. Soc., Ser. B, Stat. Methodol.*, 60(1):255-268, 1998. ISSN: 1369-7412. DOI: 10.1111/1467-9868.00123.

[42] Gareth O. Roberts and Richard L. Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. *Bernoulli*, 2(4):341-363, 1996. ISSN: 1350-7265. DOI: 10.2307/3318418.

[43] Taiji Suzuki, Atsushi Nitanda and Denny Wu. Uniform-in-time propagation of chaos for the mean-field gradient Langevin dynamics. In *11th Int. Conf. Learn. Represent*. 2023. URL: https://openreview.net/forum?id=_JScUk9TBUn.

[44] Feng-Yu Wang. Distribution dependent SDEs for Landau type equations. *Stochastic Process. Appl.*, 128(2):595-621, 2018. ISSN: 0304-4149. DOI: 10.1016/j.spa.2017.05.006.