# Phase Diagram of Dropout for Two-Layer Neural Networks in the Mean-Field Regime

Lénaïc Chizat, Pierre Marion, Yerkin Yesbay

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

To cite this version:

Lénaïc Chizat, Pierre Marion, Yerkin Yesbay. Phase Diagram of Dropout for Two-Layer Neural Networks in the Mean-Field Regime. 2025. hal-05307271

HAL Id: hal-05307271

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

Preprint submitted on 10 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
Phase Diagram of Dropout for Two-Layer Neural Networks in
the Mean-Field Regime

Lénaïc Chizat¹*, Pierre Marion²†, Yerkin Yesbay³†

¹ École Polytechnique Fédérale de Lausanne (EPFL), Institute of Mathematics, Lausanne, Switzerland

²Inria, École Normale Supérieure, PSL Research University, Paris, France

³Télécom Paris, Institut Polytechnique de Paris, Palaiseau, France

Abstract

Dropout is a standard training technique for neural networks that consists of randomly deactivating units at each step of their gradient-based training. It is known to improve performance in many settings, including in the large-scale training of language or vision models. As a first step towards understanding the role of dropout in large neural networks, we study the large-width asymptotics of gradient descent with dropout on two-layer neural networks with the mean-field initialization scale. We obtain a rich asymptotic phase diagram that exhibits five distinct nondegenerate phases depending on the relative magnitudes of the dropout rate, the learning rate, and the width. Notably, we find that the well-studied “penalty” effect of dropout only persists in the limit with impractically small learning rates of order $O(1/width)$. For larger learning rates, this effect disappears and in the limit, dropout is equivalent to a “random geometry” technique, where the gradients are thinned randomly *after* the forward and backward pass have been computed. In this asymptotic regime, the limit is described by a mean-field jump process where the neurons’ update times follow independent Poisson or Bernoulli clocks (depending on whether the learning rate vanishes or not). For some of the phases, we obtain a description of the limit dynamics both in path-space and in distribution-space. The convergence proofs involve a mix of tools from mean-field particle systems and stochastic processes. Together, our results lay the groundwork for a renewed theoretical understanding of dropout in large-scale neural networks.

# 1 Introduction

Stochastic gradient descent (SGD) is the workhorse of deep learning optimization, but its raw form is rarely used in practice. Over the years, researchers have proposed many training techniques—adaptive learning rates, gradient clipping, weight decay, and more—on top of SGD to accelerate convergence or improve test performance. Among these heuristics, *dropout* [Srivastava et al., 2014] has stood the test of time because it is easy to implement, architecture-agnostic, and the associated performance improvement is empirically robust across tasks. It consists in randomly deactivating a subset of units/neurons during each forward-backward pass. This random *thinning* of the network induces a structured, non-centered, perturbation of the gradients whose theoretical analysis remains challenging.

The original intuition for dropout's performance is that “it prevents co-adaptation between units” and that it approximates the “combined opinions of exponentially many dropout

<footnote>*Correspondence: lenaic.chizat@epfl.ch</footnote>

<footnote>†Work mostly done while at the Institute of Mathematics, EPFL.</footnote>
nets” [Srivastava et al., 2014]. However, these intuitions remain hypothetical and the recent theoretical literature on dropout follows a quite different route, which we now outline.

Let $\theta$ be the parameters of a neural network (NN), $\mathcal{L}(\theta)$ the training loss and $\mathcal{L}^{\eta}(\theta)$ the training loss perturbed with a dropout mask encoded by the random variable $\eta$. Following e.g. Wei et al. [2020], one can decompose the noise induced by dropout on gradients as

$$ \nabla \mathcal{L}^{\eta}(\theta) - \nabla \mathcal{L}(\theta) = \underbrace{\mathbf{E}_{\eta}[\nabla \mathcal{L}^{\eta}(\theta)] - \nabla \mathcal{L}(\theta)}_{\text{bias}} + \underbrace{\nabla \mathcal{L}^{\eta}(\theta) - \mathbf{E}_{\eta}[\nabla \mathcal{L}^{\eta}(\theta)]}_{\text{centered noise}}. \quad (1) $$

Exchanging gradient and expectation, we notice that the bias is of the form $\nabla P(\theta)$ with

$$ P(\theta) = \mathbf{E}_{\eta}[\mathcal{L}^{\eta}(\theta)] - \mathcal{L}(\theta). \quad (2) $$

This function, known as *dropout penalty*, can be made explicit in simple settings, such as linearly-parameterized models or two-layer NNs with squared loss [Arora et al., 2021, Proposition 4, and references therein]. More generally, an informal but empirically accurate approximation of the penalty can be given in terms of the layer-to-layer Jacobians [Wei et al., 2020]. Much of the theoretical literature on dropout studies the optimization and generalization properties of this penalty, in particular by investigating its connection with more classical regularization terms such as the weight $\ell_2$ norm or path norm [notably, Wager et al., 2013, Helmbold and Long, 2015, Mou et al., 2018, Arora et al., 2021]. As for the stochastic term in (1), its effect has been studied in Mianjy and Arora [2020], Wei et al. [2020], Zhang and Xu [2024], Zhang et al. [2024].

Despite this extensive literature, the question of the scaling of dropout to large-width networks has received little attention. This is a particularly important issue given that dropout is still used by practitioners for large-scale models [Vaswani et al., 2017, Radford et al., 2018, Devlin et al., 2019, Dosovitskiy et al., 2021, Hu et al., 2021, Ouyang et al., 2022, Ramesh et al., 2022, Taylor et al., 2022]. This begs the question of revisiting the theory of dropout in the context of large-scale networks, which is our focus in this paper. To the best of our knowledge, we are the first work to study the scaling of dropout for large NNs from a theoretical angle. A couple of experimental papers investigated specifically the question of the role of dropout in large NNs for language and vision tasks [Liu et al., 2023, Xue et al., 2023]. These studies confirm the effectiveness of dropout in these settings.

This question falls into the broader effort of understanding the scaling properties of large-scale NNs. In the absence of dropout, prior work has shown that the training dynamics of large NNs already exhibit a rich asymptotic phase diagram [Chizat et al., 2019], with in particular the so-called neural tangent kernel (NTK) regime [Jacot et al., 2018], where (overparameterized) models may reach zero training loss albeit without building data-dependent representations, and the richer feature learning [Chizat and Bach, 2018, Mei et al., 2018, Sirignano and Spiliopoulos, 2020, Rotskoff and Vanden-Eijnden, 2022] associated to the so-called mean-field scaling or $\mu$-parameterization ($\mu$P) [Yang and Hu, 2021]. These scalings were investigated in many works for various architectures [e.g. Geiger et al., 2020, Yang and Hu, 2021, Vyas et al., 2023]. A current active topic in this line of research concerns the large-depth scaling of ResNets [see Dey et al., 2025, Chizat, 2025, and references therein].

**Contributions.** This paper is the first step in a broader effort to understand the theoretical properties of dropout in large-scale models. To establish a rigorous foundation, we adopt a mathematical and descriptive approach, focusing on a restricted yet illustrative setting: two-layer networks trained with gradient descent and dropout, using the squared loss and the mean-field scaling. Despite its simplicity, this setup already reveals surprising insights and challenges several widely-held beliefs about dropout.
Our main result (Theorem 1 in Section 3) characterizes all possible scaling limits of the dynamics, depending on the relative scalings of the width, learning rate, and dropout rate. The proof, given in Section 5, relies on different constructions depending on the scaling of parameters.

In Section 4, we analyze the qualitative properties of these limiting dynamics. In particular:

* All regimes involving the penalization term (2) require a learning rate of order $O(1/\text{width})$ which is computationally impractical. A more efficient alternative that does not require small learning rates is to explicitly add a deterministic penalty $P(\theta)$ to the risk. We show in Proposition 3 the equivalence of this procedure with dropout in the small learning rate regime.

* For larger learning rates, we show that the asymptotic effect of dropout is equivalent to gradient descent in a random metric—specifically, a random block-coordinate descent scheme (Proposition 3). We investigate the potential advantages of this regime and argue that it allows certain features to take larger steps than the local sharpness of the loss would otherwise permit (see Section 4.2.2). We also show that a time-averaging of the NTK (Section 4.2.3) occurs in this regime.

Let us now briefly comment on the mathematical techniques used in our proofs. The two main technical contributions correspond to cases (II) and (III) of Theorem 1.

In case (II), we establish convergence of the dynamics to a penalized Wasserstein gradient flow. The analysis relies on coupling techniques for mean-field particle systems [see e.g. McKean, 1967, Sznitman, 1991, Lacker, 2018], as well as martingale methods to control the continuous-time limit of stochastic terms. While similar tools have been combined previously [e.g. Mei et al., 2019], our setting is different and requires careful error decomposition, as some of the relevant limits do not commute. We derive this result as a special case of a more general mean-field stochastic approximation theorem (Section 5.3.2), which may be of independent interest. This approach allows to derive convergence both in path space and in distribution space.

In case (III), we prove convergence to a mean-field jump process. This is handled by constructing an asynchronous coupling between the discrete dynamics—driven by binomial clocks—and the limiting process—driven by a Poisson clock. We also address the well-posedness of the limiting dynamics, both at the SDE and PDE levels. As we were not able to obtain a suitable uniform-in-time particle-level error control, our convergence result is formulated only in terms of time-marginal distributions, and we leave the question of pathwise convergence open for this case.

## 2 Setup

**Two-layer neural networks.** Consider a two-layer neural network (2NN) of width $n \in \mathbb{N}^*$ of the form $f(\theta, z) = \frac{1}{n} \sum_{i=1}^n \phi(x^i, z)$ where $\theta = (x^1, \dots, x^n) \in (\mathbb{R}^p)^n$ are the weights, $z \in \mathbb{R}^d$ is the input and $\phi: \mathbb{R}^p \times \mathbb{R}^d \to \mathbb{R}$. The typical structure of a 2NN is recovered when $\phi(x, z) = a\sigma(b^\top z + c)$ where $\sigma: \mathbb{R} \to \mathbb{R}$ is the activation function and $x = (a, b, c) \in \mathbb{R} \times \mathbb{R}^d \times \mathbb{R}$.

Given a training set $(z_i, y_i)_{i=1}^m \in (\mathbb{R}^d \times \mathbb{R})^m$, it will be convenient to introduce the feature map $\phi: \mathbb{R}^p \to \mathbb{R}^m$ defined by $\phi(x) = (\phi(x, z_1), \dots, \phi(x, z_m))$. Analogously, we define $f(\theta) = (f(\theta, z_1), \dots, f(\theta, z_m))$, which we call the *predictor* parameterized by $\theta$. By construction, it holds

$$f(\theta) = \frac{1}{n} \sum_{i=1}^{n} \phi(x^i), \qquad \theta = (x^1, \dots, x^n).$$
All our theoretical results are proved under the following assumption:

**Assumption A (Regularity of $\phi$).** The function $\phi: \mathbb{R}^p \to \mathbb{R}^m$ is bounded and differentiable with a bounded and Lipschitz differential. We write $D\phi(x) \in \mathbb{R}^{m \times p}$ for the Jacobian matrix of $\phi$ at $x$.

Our boundedness and regularity assumptions on $\phi$ are quite strong (in particular they require bounded weights in the case of standard 2NN presented above) but they allow focusing on the core ideas in the proofs. We believe that it could be possible to relax these assumptions at the cost of an increased level of technicality (see e.g. Wojtowytsch [2020], Chizat and Bach [2020] for the mean-field limit of ReLU 2NN, where $\phi$ satisfies none of these assumptions).

**Dropout.** In our context, dropout consists in modifying the predictor at each step of the training process by applying a random mask to the sum defining $f$. The mask is encoded by a vector $\eta = (\eta^1, \dots, \eta^n) \in \mathbb{R}^n$ with i.i.d. entries and the *dropout predictor* is

$$f^{\eta}(\theta) = \frac{1}{n} \sum_{i=1}^{n} (1 + \eta^i) \phi(x^i), \quad \text{where} \quad \eta^i = \begin{cases} \frac{1-q}{q} & \text{with probability } q \\ -1 & \text{with probability } 1-q. \end{cases}$$

The *keep rate* $q \in [0, 1]$ is a hyperparameter that represents the probability of *not* masking a neuron. It is equal to one minus the dropout rate. Observe that $(1 + \eta^i)$ equals $1/q$ with probability $q$ and equals 0 with probability $1 - q$. We will repeatedly use the following identities:

$$\mathbf{E}[\eta^i] = 0, \quad \mathbf{E}[(\eta^i)^2] = \frac{1-q}{q}, \quad \mathbf{E}[(1+\eta^i)^2] = \frac{1}{q}.$$

**GD-dropout.** We consider the squared loss on the training set $(z_i, y_i)_{i=1}^m$. Writing $\mathbf{y} = (y_1, \dots, y_m) \in \mathbb{R}^m$, this leads, in presence of dropout, to the following objective function (which inherits randomness from $\eta$):

$$\mathcal{L}^{\eta}(\theta) = \frac{1}{2} \| \mathbf{f}^{\eta}(\theta) - \mathbf{y} \|_{2}^{2}.$$

The dynamics under investigation, which we refer to as *GD-dropout*, is gradient descent (GD) on the dropout loss with $\eta$ independently resampled at each iteration. Let $\mu_0 \in \mathcal{P}(\mathbb{R}^p)$ be an initial probability distribution and $\xi^1, \xi^2, \dots$ a sequence of independent samples from $\mu_0$. For a given width $n \in \mathbb{N}^*$, we define the (random) sequence $(\theta_k^n)_{k \in \mathbb{N}}$ recursively by $\theta_0^n = (\xi^1, \dots, \xi^n)$ and

$$\theta_{k+1}^n = \theta_k^n - \tau \cdot n \cdot \nabla \mathcal{L}^{\eta_k}(\theta_k^n)$$

where $\tau > 0$ is the master learning rate. The factor $n$ multiplying the gradient is introduced so that, in the absence of dropout, the limit is nondegenerate when $n \to \infty$ (this is the mean-field/μP learning rate, see references in Section 1). The dropout mask $\eta_k \in \mathbb{R}^n$ is independently sampled at each iteration $k$. In this setting, the sources of randomness are thus (i) the random initialization and (ii) the random dropout masks.

By the chain rule, we obtain the following dynamics in terms of the weights $\theta_k^n = (X_k^{1,n}, \dots, X_k^{n,n})$. From now on, we use capital letters for the parameters of the dynamics, to recall that we are dealing with random variables. For $i \in [1:n]$:

$$\left\{ \begin{array}{l} X_0^{i,n} = \xi^i, \\ X_k^{i,n} = X_{k-1}^{i,n} - \tau \cdot (1 + \eta_k^i) D\phi(X_{k-1}^{i,n})^\top \left( \frac{1}{n} \sum_{i=1}^n (1 + \eta_k^i) \phi(X_{k-1}^{i,n}) - \mathbf{y} \right), \quad k \in \mathbb{N}^*. \end{array} \right. \quad (3)$$

The rest of the paper is concerned with the analysis of this dynamical system.
# 3 Main result: characterization of the infinite-width limits

There are three hyperparameters (HP) involved in the definition the dynamics: width $n$, keep rate $q$, and learning rate $\tau$. Our goal is to describe how the dynamics behaves in the large-width limit $n \to \infty$ and as a function of the scaling of other HP with $n$.

## 3.1 Notation for infinite-width 2NN

To formulate the limit dynamics, we allow $f$ to be parameterized by a probability measure (thereby overloading notations). For $\mu \in \mathcal{P}(\mathbb{R}^p)$, let

$$f(\mu) = \int_{\mathbb{R}^p} \phi(x) d\mu(x).$$

A finite-width 2NN with weights $\theta = (x^1, \dots, x^n)$ is obtained as the particular case $f(\hat{\mu}^n)$ where $\hat{\mu}^n = \frac{1}{n} \sum_{i=1}^n \delta_{x^i}$ is the empirical distribution of weights. This parameterization by a measure leads to an objective function $\mathcal{L}: \mathcal{P}(\mathbb{R}^p) \to \mathbb{R}$

$$\mathcal{L}(\mu) = \frac{1}{2} \| \mathbf{f}(\mu) - \mathbf{y} \|_2^2.$$

Before stating our main result, let us introduce some useful definitions and notations.

1. The mean potential $V: \mathcal{P}(\mathbb{R}^p) \times \mathbb{R}^p \to \mathbb{R}$ is defined as $V[\mu](x) = \phi(x)^\top (\mathbf{f}(\mu) - \mathbf{y})$. This function represents the Fréchet derivative of $\mathcal{L}$, seen as a functional over finite signed measures, in the sense that for any finite signed measure $\sigma$ on $\mathbb{R}^p$, it holds as $\varepsilon \to 0$

$$\mathcal{L}(\mu + \varepsilon\sigma) = \mathcal{L}(\mu) + \varepsilon \int V[\mu](x)d\sigma(x) + o(\varepsilon).$$

In the following, we will frequently make use of the gradient of $V$ with respect to its second argument, that is

$$\nabla V : \mathcal{P}(\mathbb{R}^p) \times \mathbb{R}^p \to \mathbb{R}^p, \quad \nabla V[\mu](x) = D\phi(x)^\top (\mathbf{f}(\mu) - \mathbf{y}).$$

2. Given a probability measure $\mu \in \mathcal{P}(\mathbb{R}^p)$ and a measurable map $T: \mathbb{R}^p \to \mathbb{R}^p$, the pushforward measure $T_{\#}\mu$ is defined as $T_{\#}\mu(B) = \mu(T^{-1}(B))$ for all Borel sets $B \subset \mathbb{R}^p$. We also denote by $\nabla \cdot$ the divergence operator, which is the opposite of the adjoint of the gradient operator.

3. Given a Polish space $(\mathcal{X}, \mathrm{dist})$ (such as $\mathbb{R}^p$, see more details in Sec. 5.1), the $L^1$- Wasserstein distance on the space $\mathcal{P}_1(\mathcal{X})$ of probability measures with finite first moment is defined as

$$W_1(\mu, \nu) = \sup_{\psi \in \mathrm{1-Lip}(\mathcal{X})} \int \psi d(\mu - \nu) = \inf_{\pi \in \Pi(\mu, \nu)} \int \mathrm{dist}(x, y)d\pi(x, y)$$

where the supremum runs over $\mathrm{1-Lip}(\mathcal{X})$ the set of $\mathrm{1-Lipschitz}$ functions from $\mathcal{X}$ to $\mathbb{R}$, and the infimum runs over the set of transport plans $\Pi(\mu, \nu) = \{\pi \in \mathcal{P}(\mathcal{X} \times \mathcal{X}) ; \int \pi(\cdot, dy) = \mu$ and $\int \pi(dx, \cdot) = \nu\}$. The space $(\mathcal{P}_1(\mathcal{X}), W_1)$ is a Polish space [Villani, 2008, Chap. 6].
## 3.2 Main result

In this section, we describe the various limits for dropout-GD in the large-width limit $n \to \infty$. To this end, we consider two positive sequences of learning rates ($\tau_n$) and keep rates ($q_n$) indexed by the width $n$. Our main result deals with nondegenerate limits, which are obtained under the following conditions.

**Assumption B (Nondegenerate limits).** Assume that the sequences $\tau_n$, $q_n$, $\alpha_n = \frac{\tau_n}{q_n}$ and $\beta_n = \frac{1}{nq_n}$ admit finite limits, denoted $\tau, q, \alpha$ and $\beta$ respectively.

In addition to the sequences ($\tau_n$) and ($q_n$) already introduced, the two new sequences ($\alpha_n$) and ($\beta_n$) appearing in this definition are key to capture the asymptotic behavior of the dynamics:

* The scalar $\alpha_n$ represents the expected time between two updates of the same neuron, as measured in the natural time unit given by $t = k\tau_n$. Indeed, if $\eta_k^i > 0$ and $k' = \min\{\ell > k; \eta_\ell^i > 0\}$ then $k' - k$ is a geometric random variable with expectation $1/q_n$, which multiplied by $\tau_n$ gives $\alpha_n$. The scalar $\alpha_n$ is also equal to the effective learning rate, that is, the scalar factor multiplying the gradient of active units in (3). If $\alpha_n \to +\infty$ then the dynamics degenerates in two ways: (i) over a fixed time horizon $[0, T]$, only a vanishing fraction of the neurons are active and (ii) active neurons are subject to a diverging effective learning rate<fnref n="1" />. We thus assume that $\alpha = \lim_n \alpha_n < +\infty$.

* The scalar $\beta_n$ represents the inverse of the expected number of neurons which are updated at a given step. Indeed, this number follows a Bernoulli distribution with parameters $n$ and $q_n$, of expectation $nq_n$. When $\beta_n \to +\infty$, in the limit, almost all the GD steps are trivial with no active neuron<fnref n="2" />. We thus assume $\beta = \lim_n \beta_n < +\infty$.

With these key parameters in place, we are ready for our main result, which classifies all possible asymptotic behaviors of the empirical distribution of parameters $\hat{\mu}_k^n = \frac{1}{n} \sum_{i=1}^n \delta_{X_k^{i,n}}$ in the limit $n \to \infty$ (to the exception of the degenerate cases just mentioned).

This theorem is a corollary of more precise statements given in Section 5 that also describe the asymptotic dynamics of individual neurons. We refer to Figure 1 for a graphical illustration of the phase diagram.

**Theorem 1 (Main theorem).** Let Assumption A and B hold. Consider the sequence of GD-dropout dynamics (3) with $\mu_0 \in \mathcal{P}_1(\mathbb{R}^p)$. Then the limiting dynamics can be classified in four cases depending on the scaling of HP:

(I) (Discrete-time jump process limit) If $\tau_n \to \tau > 0$ (and hence $q, \alpha > 0$ and $\beta = 0$), then for all $k \in \mathbb{N}$ it holds $\mathbf{E}[W_1(\hat{\mu}_k^n, \rho_k)] \to 0$, where the limit sequence $(\rho_k)_{k \in \mathbb{N}}$ is defined recursively by $\rho_0 = \mu_0$ and

$$ \rho_{k+1} = q (\mathrm{Id} - \alpha \nabla V[\rho_k])_{\#} \rho_k + (1-q)\rho_k. $$

(II) (Wasserstein Gradient flow limit (penalized if $\beta > 0$)) If $\tau_n \to 0$ and $\alpha = 0$, then for any $t \ge 0$, $\mathbf{E}[W_1(\hat{\mu}_{[t/\tau_n]}^n, \rho_t)] \to 0$ where $\rho \in \mathcal{C}(\mathbb{R}_+; \mathcal{P}_1(\mathbb{R}^p))$ is the unique weak solution to

$$ \partial_t \rho_t = \nabla \cdot (\rho_t(\nabla V[\rho_t] + \beta \nabla P)), \quad \rho_0 = \mu_0 $$

<footnote n="1">¹The overall effect on the predictor of these two opposite trends – fewer but larger updates – depends in general on the growth of $\phi$; but under our assumption that $\phi$ is bounded, the limit dynamics when $\alpha_n \to \infty$ is stationary in predictor space because the effect of the updated units is overall negligible.</footnote>

<footnote n="2">²In this case, we can condition on having at least 1 neuron active – which means skipping trivial steps – which leads, up to time rescaling, to an equivalent dynamics with $\beta \le 1$. Therefore, it would not be restrictive to assume $\beta \le 1$ (we only assume $\beta < +\infty$).</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: Phase diagram of dropout in two-layer NN with mean-field scaling. For an HP scaling $(n, \tau_n^{-1}, q_n^{-1})$, the limit of $(\log n, -\log \tau_n, -\log q_n)/S_n$ where $S_n = \log n - \log \tau_n - \log q_n$ (when it exists) forms a point in the 2-simplex, which is represented in the triangle in barycentric coordinates. For instance, the central red point represents proportional limits while the blue line corresponds to $\tau_n^{-1}$ and $q_n^{-1}$ diverging proportionally while $n$ diverges faster. Red area: degenerate limits with either $\alpha = +\infty$ or $\beta = +\infty$. Grey area: Wasserstein gradient flow limit (the effect of dropout disappears in the limit). Orange vertex: discrete-time jump process limit (if $\tau, q > 0$). Blue line: continuous-time jump process limit. Red vertex: critical limit. Green line: penalized Wasserstein gradient flow limit.</figcaption></figure>

where $P(x) = \frac{1}{2}\|\phi(x)\|_2^2$ is the dropout penalty term.

(III) (Continuous-time jump process limit) If $\tau_n \to 0$, $\alpha > 0$ and $\beta = 0$, then for any $t \ge 0$, $\mathbf{E}[W_1(\hat{\mu}_{[t/\tau_n]}^n, \rho_t)] \to 0$ where $\rho \in \mathcal{C}(\mathbb{R}_+; \mathcal{P}_1(\mathbb{R}^p))$ is the unique weak solution to

$$ \partial_t \rho_t = \alpha^{-1}((\mathrm{Id} - \alpha \nabla V[\rho_t])_{\#} \rho_t - \rho_t), \quad \rho_0 = \mu_0. \quad (4) $$

(IV) (Critical limit) Otherwise $\tau_n \to 0$, $\alpha > 0$ and $\beta > 0$. This case is discussed below.

*Proof.* The cases (I), (II) and (III) are corollaries of Proposition 6, Proposition 7 and Proposition 14 respectively. Also, one can check that the dichotomy exhausts all cases under Assumption B (this is clear from the phase diagram on Figure 1). □

Regarding the case (IV) our ansatz is that the limit $(\rho_t)$ solves

$$ \partial_t \rho_t = \alpha^{-1} \left( \int K[\rho_t](\cdot|x) d\rho_t(x) - \rho_t \right), \quad \rho_0 = \mu_0 \quad (5) $$

where the transition kernel $K[\rho]$ is defined as

$$ K[\rho](\cdot|x) = \mathrm{Law}(x - \alpha \nabla V[\hat{M}(\rho)](x) - \alpha \beta \nabla P(x)). $$

Here $\hat{M}: \mathcal{P}_1(\mathbb{R}^p) \to \mathcal{P}_1(\mathbb{R}^p)$ is the random map defined as $\hat{M}(\rho) = \beta \sum_{i=1}^N \delta_{Z_i}$ where $N \sim \mathrm{Poisson}(1/\beta)$ and $Z_1, \dots, Z_N$ are independent samples from $\rho$. In the definition of $K[\rho]$, the law is over the randomness of $\hat{M}$, which is the only source of randomness in this expression. We formally justify this ansatz in Section 4.3, but we leave a rigorous analysis of
this case as an open question. Note that this critical limit corresponds to learning rates of order $\Theta(1/n)$ and is therefore of limited practical relevance.

The proof of Theorem 1 is deferred to Section 5 and a discussion of the various limits is the object of Section 4. For now, we can make the following short remarks:

1. When proving Theorem 1, we in fact obtain in cases (I) and (II) the stronger notion of convergence at the level of particle's trajectories. In mathematical terms, for any $m \in \mathbb{N}^*$ fixed, we obtain the convergence in distribution of the family of paths $(X^{1,n}, \dots, X^{m,n})$ towards $(Y^1, \dots, Y^m)$ where $(Y^i)$ are independent realizations of a limit process given in Section 5.

2. Interestingly, in all nondegenerate limits, $f(\rho_t)$ evolves deterministically.

3. Our results are proved for the squared loss, but it would not be difficult to extend them to any smooth loss when $\beta = 0$. When $\beta > 0$, the form of the limit would change because the dropout penalty does not have a simple closed form beyond the squared loss.

# 4 Qualitative properties of the various limits

In this section, we discuss some properties of the asymptotic processes presented previously.

## 4.1 The penalized Wasserstein gradient flow limit (II): equivalence with explicit penalization

The only two limits that exhibit the dropout penalty, that is when dropout noise remains asymptotically biased (see (2)), are cases (IV) and (II) with $\beta > 0$. These cases involve small learning rates $\tau_n = O(1/n)$ and keep rates $q_n = \Theta(1/n)$. Such learning rates are impractical in large-scale settings as they considerably slow down training.

If the penalty effect is indeed beneficial and desired, a more practical approach is to dispense with dropout altogether and instead add the deterministic penalty $P(\theta)$ directly to the objective. In the next proposition, we show that this procedure is asymptotically equivalent to dropout in regimes where the penalty emerges, focusing for simplicity on case (II), that is, with $\tau_n = o(1/n)$. However, unlike implementing the penalty via standard dropout, this approach does not require vanishing learning rates and is therefore more widely applicable, as well as being simple to implement.

**Proposition 2** (Explicit-implicit penalty equivalence). Let $\beta > 0$. Consider an initialization $(\xi^1, \dots, \xi^n)$ as in (3) and let:

(i) $(X_k^{i,n})$ the iterates of GD with dropout on the loss $\mathcal{L}^\eta$ with $q_n = \frac{1}{\beta n}$ and $\tau_n = o(\frac{1}{n})$;

(ii) $(\tilde{X}_k^{i,n})$ the iterates of GD without dropout ($q = 1$) on the penalized loss $\theta \mapsto \mathcal{L}(\theta) + \frac{\beta}{n} \sum_{i=1}^n P(X^{i,n})$ with the same $(\tau_n)$;

Then we have pathwise convergence over bounded time horizons: for any $T > 0$ and $m \in \mathbb{N}^*$, it holds

$$ \mathbf{E} \left[ \sup_{k \le T/\tau_n} \sup_{i \le \min\{m,n\}} \| \tilde{X}_k^{i,n} - X_k^{i,n} \|_2 \right] \xrightarrow{n \to \infty} 0. $$

In particular the respective distributional limits $(\rho_t)$ and $(\tilde{\rho}_t)$ exist and are equal $(\rho_t) = (\tilde{\rho}_t)$.

*Proof.* This follows from the fact that the proof of Proposition 7 applies to $X_k^{i,n}$ but also to $\tilde{X}_k^{i,n}$ (there are just fewer error terms in this case). So both dynamics converge pathwise to the limit process denoted $(Y^1, \dots, Y^m)$ in that statement. $\square$
## 4.2 The jump process limits (I) & (III)

In this section, we discuss properties of the jump process limits. Since the continuous-time jump process limit (III) is simply the small learning rate limit of the discrete-time jump process limit (I), we show properties for whichever of these two limits is more convenient; keeping in mind that an analogous property could be studied in the other limit.

### 4.2.1 Equivalence with the random metric (RaM) technique

First, let us show that in these asymptotics, dropout can be equivalently replaced by another, simpler training technique which consists in randomly removing units only after the gradient has been computed. In the two-layer case, this gives the following dynamics with $i \in [1:n]$ (compare with (3)),

$$ \left\{ \begin{aligned} \tilde{X}_0^{i,n} &= \xi^i, \\ \tilde{X}_k^{i,n} &= \tilde{X}_{k-1}^{i,n} - \tau_n \cdot (1 + \eta_k^i) D\phi(\tilde{X}_{k-1}^{i,n})^\top \left( \frac{1}{n} \sum_{i=1}^n \phi(\tilde{X}_{k-1}^{i,n}) - \mathbf{y} \right), \quad k \in \mathbb{N}^* \end{aligned} \right. \quad (6) $$

where $(\eta_k^i)$ is as in (3). This corresponds to GD but with a random diagonal (degenerate) metric (i.e. inverse preconditioner) with entries $(1 + \eta_k^i)^{-1}$, which takes the two values $\alpha_n^{-1}$ and $+\infty$. The updates are equivalent to random block-coordinate descent on the loss, an optimization technique that is usually exploited to leverage the lower computational complexity of sparse updates. Since the sparsity is not the motivation here and since the idea can be generalized to metrics with more than two values $(\alpha_n^{-1}$ and $+\infty)$, we prefer to refer to this as the random metric (RaM) technique.

Let us state the equivalence in the discrete-time case. This result is empirically illustrated on Figure 2.

**Proposition 3 (Dropout-RaM equivalence).** Let $\tau, q > 0$ and $\alpha = \tau/q$. Consider, with the same initialization $(\xi^1, \dots, \xi^n)$ as in (3):

(i) $(X_k^{i,n})$ the iterates of GD with dropout on the loss $\mathcal{L}^\eta$ with fixed dropout rate $q$ and learning rate $\tau$;

(ii) $(\tilde{X}_k^{i,n})$ the RaM iterates as in (6) with $\tau_n = \tau$ and the same dropout masks $(\eta_k^i)$;

Then we have pathwise convergence over bounded time horizons: for any $K > 0$ and $m \in \mathbb{N}^*$, it holds

$$ \mathbf{E} \left[ \sup_{k \le K} \sup_{i \le \min\{m,n\}} \| \tilde{X}_k^{i,n} - X_k^{i,n} \|_2 \right] \xrightarrow{n \to \infty} 0. $$

In particular the respective distributional limits $(\rho_k)$ and $(\tilde{\rho}_k)$ exist and are equal $(\rho_k) = (\tilde{\rho}_k)$.

*Proof.* This follows from the fact that the proof of Proposition 6 applies to $X_k^{i,n}$ but also to $\tilde{X}_k^{i,n}$ (there are just fewer error terms in this case). So both dynamics converge pathwise to the limit process denoted $(Y^1, \dots, Y^m)$ in that statement. $\square$

### 4.2.2 A benefit of RaM: large feature updates beyond the sharpness bound

RaM has an interesting consequence in the large learning rate regime: it makes it possible for certain neurons to take larger steps than what the sharpness of the loss would otherwise permit, as long as the average effective learning rate is small enough. To present this idea,
<figure><img src="image_4.png" /><figcaption>Figure 2: Illustration of the pathwise convergence between random metric (RaM) and dropout dynamics (Proposition 3). We train a width-5000 two-layer NN on a synthetic teacher-student task with the quadratic loss, either with GD (orange), GD with dropout (blue) or GD with RaM (green), with coupled randomness (of initialization and masks ($\eta_k^i$)). The paths shown correspond to two-dimensional projections of the trajectory of two randomly chosen neurons. The pathwise similarity between dropout and RaM illustrates Proposition 3 (the similarity degrades at large times because the width is finite). More details and additional plots are given in Appendix A.1.</figcaption></figure>

let us study the evolution of the loss along one random-metric GD step in the infinite-width limit.

For the purpose of this discussion, instead of parameterizing the neural network by a measure $\mu \in \mathcal{P}(\mathbb{R}^p)$, it is more convenient to parameterize it by a random variable $X \in L^2(\Omega; \mathbb{R}^p)$, in the spirit of Lions calculus [Cardaliaguet, 2010, Section 6.1] (here $\Omega$ is an abstract probability space). These two parameterizations lead to the same predictor if $\text{Law}(X) = \mu$. In this parameterization, for $\mathbb{R}^p$-valued $L^2$ random variables $X, Y_1, Y_2$, we have the following loss and $L^2$-derivatives [Cardaliaguet, 2010, Section 6.1]

$$ 
\begin{aligned}
\mathcal{L}(X) &= \frac{1}{2} \|\mathbf{E}[\phi(X)] - \mathbf{y}\|_2^2, \\
D\mathcal{L}(X)(Y_1) &= (\mathbf{E}[\phi(X)] - \mathbf{y})^\top \mathbf{E}[D\phi(X)[Y_1]] \\
D^2\mathcal{L}(X)(Y_1, Y_2) &= \mathbf{E}[D\phi(X)[Y_1]]^\top \mathbf{E}[D\phi(X)[Y_2]] + (\mathbf{E}[\phi(X)] - \mathbf{y})^\top \mathbf{E}[D^2\phi(X)[Y_1, Y_2]]
\end{aligned}
$$

From the expression of the first differential, we obtain that the gradient's expression is $\nabla \mathcal{L}(X) = D\phi(X)^\top (\mathbf{E}[\phi(X)] - \mathbf{y})$ (this is a bounded random variable under Assumption A). Now given $X$ and an independent nonnegative $L^2$ random variable $\eta$, the random geometry update is of the form $X_+ = X - \eta \nabla \mathcal{L}(X)$.

In these notations, the one-step evolution of the loss is given, up to second-order terms, by

$$ \mathcal{L}(X_+) \approx \mathcal{L}(X) - D\mathcal{L}(X)(\eta \nabla \mathcal{L}(X)) + \frac{1}{2} D^2\mathcal{L}(X)(\eta \nabla \mathcal{L}(X), \eta \nabla \mathcal{L}(X)) $$

Ignoring higher order terms, the decrease of the objective is guaranteed provided $S(X) \le 1$ with

$$ S(X) := \frac{D^2\mathcal{L}(\eta \nabla \mathcal{L}(X), \eta \nabla \mathcal{L}(X))}{2D\mathcal{L}(X)(\eta \nabla \mathcal{L}(X))}. $$

Using the independence of $\eta$, the denominator is given by

$$ 2D\mathcal{L}(X)(\eta \nabla \mathcal{L}(X)) = 2\mathbf{E}[\eta](\mathbf{E}[\phi(X)] - \mathbf{y})^\top \mathbf{E}[D\phi(X)D\phi(X)^\top](\mathbf{E}[\phi(X)] - \mathbf{y}) = A(X)\mathbf{E}[\eta] $$
for some $A(X) \ge 0$. The numerator is the sum of two terms (I) + (II) with

$$ 
\begin{aligned}
(I) &= \frac{1}{2} \mathbf{E}[D\phi(X)[\eta\nabla\mathcal{L}(X)]]^\top \mathbf{E}[D\phi(X)[\eta\nabla\mathcal{L}(X)]] \\
&= \frac{1}{2} \mathbf{E}[\eta]^2 \|\mathbf{E}[D\phi(X)D\phi(X)^\top(\mathbf{E}[\phi(X)] - \mathbf{y})]\|^2 \\
&= B(X)\mathbf{E}[\eta]^2
\end{aligned}
$$

for some $B(X) \ge 0$ and

$$ 
\begin{aligned}
(II) &= \frac{1}{2}(\mathbf{E}[\phi(X)] - \mathbf{y})^\top \mathbf{E}[D^2\phi(X)[\eta\nabla\mathcal{L}(X), \eta\nabla\mathcal{L}(X)]] \\
&= \frac{1}{2}\mathbf{E}[\eta^2](\mathbf{E}[\phi(X) - \mathbf{y}])^\top \mathbf{E}[D^2\phi(X)[D\phi(X)^\top(\mathbf{E}[\phi(X)] - \mathbf{y}), D\phi(X)^\top(\mathbf{E}[\phi(X)] - \mathbf{y})]] \\
&= C(X)\mathbf{E}[\eta^2]
\end{aligned}
$$

for some $C(X) \in \mathbb{R}$. It follows that

$$ S(X) = \frac{B(X)\mathbf{E}[\eta]^2 + C(X)\mathbf{E}[\eta^2]}{A(X)\mathbf{E}[\eta]} $$

If $\eta = \eta_0$ is deterministic, then the descent condition $S(X) \le 1$ is $\eta_0 \le \eta_{\max} := \frac{A(X)}{B(X)+C(X)}$ which is a local version of the usual smoothness-based descent lemma in optimization.

But when $\eta$ is random, the factor in front of the two terms of the Hessian is different. Moreover, it is a pervasive empirical observation in the context of neural networks [e.g. Wei et al., 2020] that it typically holds $B(X) \gg |C(X)|$. This means that provided one chooses $\mathbf{E}[\eta]$ a bit smaller than $\eta_{\max}$, there is ample room to increase the variance of $\eta$—and therefore allowing randomly chosen neurons to take large steps, beyond the sharpness threshold—while maintaining the descent property. We note that this discussion is valid as soon as $\mathbf{E}[\eta] > 0$ and $\mathbf{E}[\eta^2] < +\infty$, in particular the nonnegativity of $\eta$ is not essential.

### 4.2.3 Another asymptotic property: dynamical NTK averaging

Let us mention another interesting property of dropout in the jump process limits: the associated Neural Tangent Kernel (NTK) is smoothed in time via an exponential moving average. Let us present this phenomenon, which is an intuitive consequence of the fact that only a subset of the neurons are updated at each step, for the continuous-time dynamics (case (III)) where its expression is the cleanest.

As shown in Proposition 13, the limit dynamics ($\rho_t$) in case (III) can be characterized as the unique solution to the following integral equation

$$ \rho_t = \frac{1}{\alpha} \int_0^t e^{(s-t)/\alpha} T_\alpha(\rho_s) ds + e^{-t/\alpha} \mu_0 \quad (7) $$

where $T_\alpha(\mu) = (\mathrm{Id} - \alpha \nabla V[\mu])_\#\mu$ is the measure obtained from $\mu$ by taking a single GD step with learning rate $\alpha$ on the loss $\mathcal{L}$.

The expression has a counterpart in terms of the NTK. Indeed, the tangent kernel<fnref n="3" /> of a 2NN is a linear function of the measure $\mu \in \mathcal{P}(\mathbb{R}^p)$ that parameterizes it:

$$ \text{NTK}(\mu) = \int D\phi(w)D\phi(w)^\top d\mu(w). $$

<footnote n="3">³With a slight abuse of notation, we identify the NTK with its Gram matrix on the training set.</footnote>
After one gradient step with learning rate $\alpha$, the updated tangent kernel is NTK($T_\alpha(\mu)$). Hence, by (7), we see that the tangent kernel at time $t$ writes as the exponential moving average

$$ \text{NTK}(\rho_t) = \frac{1}{\alpha} \int_0^t e^{(s-t)/\alpha} \text{NTK}(T_\alpha(\rho_s))ds + e^{-t/\alpha} \text{NTK}(\rho_0). $$

### 4.3 Informal justification of the ansatz for the critical limit (IV)

In this section, we justify informally the ansatz given in (5) for the limit in the critical case $\tau_n \to 0$, $\alpha, \beta > 0$.

At each timestep, there is a probability $q_n$ that a neuron moves and $1 - q_n$ that it stays still. When $n$ is large, the distribution of neurons at iteration $k+1$ is therefore of the form $\hat{\mu}_{k+1}^n = q_n\bar{\mu}_{k+1}^n + (1 - q_n)\hat{\mu}_k^n$ where $\bar{\mu}_{k+1}^n$ is the distribution of the new position of the active neurons. Reorganizing this update as $\frac{1}{\tau_n}(\hat{\mu}_{k+1}^n - \hat{\mu}_k^n) = \frac{q_n}{\tau_n}(\bar{\mu}_{k+1}^n - \mu_k^n)$, this suggests a limit equation of the form

$$ \partial_t \rho_t = \frac{1}{\alpha}(\bar{\rho}_t - \rho_t) $$

where $\bar{\rho}_t$ is the distribution of the new positions of the neurons that have moved at time $t$ (this is the same form as the jump process limit (4)).

To derive the expression for $\bar{\rho}_t$, let us take a closer look at the update equation (3) for a single neuron. Conditioning on the neuron $i$ moving at iteration $k$ (that is, conditioning on $(1 + \eta_{k+1}^i) = 1/q_n$), its update can be decomposed as

$$ \begin{aligned} X_{k+1}^{i,n} - X_k^{i,n} &= -\tau_n \cdot (1 + \eta_{k+1}^i) D\phi(X_k^{i,n})^\top \left( \frac{1}{n} \sum_{j=1}^n (1 + \eta_{k+1}^j) \phi(X_k^{j,n}) - \mathbf{y} \right) \\ &= -\frac{\tau_n}{q_n} D\phi(X_k^{i,n})^\top \underbrace{\left( \frac{1}{n} \sum_{j \neq i} (1 + \eta_{k+1}^j) \phi(X_k^{j,n}) - \mathbf{y} \right)}_{(I)} - \underbrace{-\frac{\tau_n}{nq_n^2} D\phi(X_k^{i,n})^\top \phi(X_k^{i,n})}_{(II)} \end{aligned} $$

The second term (II) is simple: as $n \to \infty$, it converges to $-\alpha\beta\nabla P(X_k^{i,n})$. The first term (I) is more subtle because the number of terms in the sum follows a Binom($n-1, q_n$) distribution, whose expectation $(n-1)q_n \to 1/\beta$ does not diverge, so there is no limit theorem taking place to simplify this sum. Asymptotically, the sum behaves as drawing $N$ according to a Poisson distribution with parameter $1/\beta$ (this is the limit of Binom($n-1, q_n$) distributions) and drawing $N$ independent samples $Z_1, \dots, Z_N$ according to $\rho_t$ to form the sum $1/(nq_n) \sum_{j=1}^N \phi(Z_j)$. Hence the first term as a whole (I) is, in the limit, sampled via the formula $-\alpha\nabla V[\hat{M}(\rho_t)](X_k^{i,n})$ where $\hat{M}$ is the operator defined below (5).

Finally, the update of all the neurons that move at iteration $k$ are certainly not independent; however, in the limit, since $\tau_n \to 0$ there are infinitely many such moves over an infinitesimal time interval $[t, t+\delta t]$, and they are independent in time (considering that $\rho_t$ changes negligibly over this time interval). This leads to an averaging-over-time effect: the neurons located near $x$ and that move during this time interval are distributed, at time $t+\delta t$, according to

$$ K[\rho_t](\cdot|x) = \text{Law}(x - \alpha\nabla V[\hat{M}(\rho_t)](x) - \alpha\beta\nabla P(x)) $$

which is deterministic. This leads to the expression $\bar{\rho}_t = \int K[\rho_t](\cdot|x)d\rho_t(x)$ and concludes the justification of ansatz (5).
## 4.4 A synthesis: the three effects of dropout and experimental illustration

Our theory allows identifying three distinct effects of dropout, with each limit in the phase diagram characterized by the combination of one or several of these effects. This can be seen by suitably decomposing the update equation (3), which we recall is

$$X_{k+1}^{i,n} = X_k^{i,n} - \tau_n \cdot (1 + \eta_{k+1}^i) D\phi(X_k^{i,n})^\top \left( \frac{1}{n} \sum_{j=1}^n (1 + \eta_{k+1}^j) \phi(X_k^{j,n}) - \mathbf{y} \right).$$

Introducing an independent copy $\tilde{\eta}_{k+1}^j$ of $\eta_{k+1}^j$, we get

$$X_{k+1}^{i,n} = \underbrace{X_k^{i,n} - \tau_n D\phi(X_k^{i,n})^\top (\mathbf{f}(\hat{\mu}_k^n) - \mathbf{y})}_{\text{update without dropout}} - \underbrace{\frac{\tau_n(1 + \eta_{k+1}^i)}{n} \sum_{j=1}^n \tilde{\eta}_{k+1}^j D\phi(X_k^{i,n})^\top \phi(X_k^{j,n})}_{\text{PN}} \\ - \underbrace{\tau_n \eta_{k+1}^j D\phi(X_k^{i,n})^\top (\mathbf{f}(\hat{\mu}_k^n) - \mathbf{y})}_{\text{RaM}} + \underbrace{\tau_n \frac{\eta_{k+1}^i \tilde{\eta}_{k+1}^i - (\eta_{k+1}^i)^2}{n} D\phi(X_k^{i,n})^\top \phi(X_k^{j,n})}_{\text{Penalty}}.$$

Let us now give an interpretation of each of the terms:

1. **Propagation noise (PN) - active in case (IV).** This is the perturbation of the signal in the forward pass due to the random thinning and rescaling of the forward pass. It manifests itself in the randomness of the predictor. In our context, it is a centered perturbation of the update of units of variance $\Theta(\frac{\tau^2}{n^2}\mathbf{E}[(1+\eta)^2]\mathbf{E}[(\sum_i \eta^i)^2]) = \Theta(\tau\alpha\beta)$.

2. **Random metric (RaM) - active in cases (I), (III), and (IV).** This is, given a full gradient, the effect of its random thinning and rescaling. It is due to the mask applied in the backward pass. It leads to a centered perturbation of the update of all units of variance $\Theta(\tau^2\mathbf{E}[\eta^2]) = \Theta(\tau\alpha)$.

3. **Dropout penalty - active in cases (II) and (IV).** This is the bias of the dropout gradient, which is precisely the effect of the non-independence between the masks used in the forward and the backward pass. In the limit, the randomness of this term vanishes, and it leads to a deterministic perturbation of the update of all units of average magnitude $\Theta(\tau\beta)$ (see proof of Proposition 7 for details).

This decomposition can be used to intuit the phase diagram of Figure 1, by separating the cases according to whether $\alpha$ only is positive (case (I) with nonvanishing learning rate and case (III) with vanishing learning rate) or $\beta$ only is positive (case (II) with $\beta > 0$) or both (case (IV)) or neither (case (II) with $\beta = 0$).

To give experimental insight on the effect of each of these terms, we plot in Figure 3 the training logs for a width 4000 two-layer ReLU NN trained on binary classification on two classes of MNIST with the logistic loss (see details in Appendix A.2). The training configurations illustrate various combinations of the three effects presented above:

* "PN" refers to dropout applied in the forward pass only, and shows the isolated effect of propagation noise.

* "RaM" refers to dropout applied in the backward pass only, and shows the isolated effect of RaM.

* "PN + RaM" refers to dropout applied in both the forward and the backward passes, but with independent masks. This shows the combined effect of propagation noise and RaM. The difference with "dropout" is precisely the effect of the dropout penalty.
<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison of variants of dropout for the training of two-layer NNs on two classes of MNIST with logistic loss for several variants of dropout (see details in text), as a function of training steps (left) and of keep probability (right). More details are given in Appendix A.2.</figcaption></figure>

• “dropout” corresponds to standard dropout, which combines the effect of PN, RaM, and dropout penalty, while “no dropout” refers to vanilla SGD.

We observe that RaM almost achieves the performance of dropout. Adding dropout in the forward pass does not improve, or even degrades, performance, unless one uses the same masks in the forward and backward passes—that is standard dropout. In the early phase of training, the dropout, RaM and PN+RaM curves match, which suggests that the speed-up observed for dropout in this early phase is entirely due to RaM.

# 5 Proof of Theorem 1

Let us begin with some preliminary results which will be useful in all limits.

## 5.1 Preliminaries

The following consequence of Assumption A will be used repeatedly in the proofs, often without explicit reference.

**Lemma 4.** Under Assumption A, there exists $C > 0$ such that, for any $\mu, \mu' \in \mathcal{P}_1(\mathbb{R}^p)$ and $x, x' \in \mathbb{R}^p$ it holds

$$ \|\nabla V[\mu](x) - \nabla V[\mu'] (x')\|_2 \le C (\|x - x'\|_2 + W_1(\mu, \mu')), \quad (8) $$

$$ \text{and } \|\nabla V[\mu](x)\|_2 < C, \|\nabla P(x)\|_2 \le C, \mathcal{L}(\mu) < C. $$

*Proof.* We recall that Assumption A states that $\phi$ and $D\phi$ are uniformly bounded. The last three bounds directly follow from the expressions $f(\mu) = \int \phi d\mu$, $\nabla V[\mu](x) = D\phi(x)^{\top}(f(\mu) -
y), $\nabla P(x) = D\phi(x)^{\top}\phi(x)$ and $\mathcal{L}(\mu) = \frac{1}{2}\|\mathbf{f}(\mu) - \mathbf{y}\|_2^2$. To prove (8), we proceed as follows:

$$
\begin{align*}
\|\nabla V[\mu](x) - \nabla V[\mu](x')\|_2 &= \|D\phi(x)^{\top}(\mathbf{f}(\mu) - \mathbf{y}) - D\phi(x')^{\top}(\mathbf{f}(\mu') - \mathbf{y})\|_2 \\
&\le \|D\phi(x)^{\top}(\mathbf{f}(\mu) - \mathbf{y}) - D\phi(x')^{\top}(\mathbf{f}(\mu) - \mathbf{y})\|_2 \\
&\quad + \|D\phi(x')^{\top}(\mathbf{f}(\mu) - \mathbf{y}) - D\phi(x')^{\top}(\mathbf{f}(\mu') - \mathbf{y})\|_2 \\
&\le \|D\phi(x) - D\phi(x')\|_{2 \to 2} \cdot \|\mathbf{f}(\mu) - \mathbf{y}\|_2 \\
&\quad + \|D\phi(x')\|_{2 \to 2} \|\mathbf{f}(\mu) - \mathbf{f}(\mu')\|_2 \\
&\le \mathrm{Lip}(D\phi) \cdot \sqrt{2\mathcal{L}(\mu)} \cdot \|x - x'\|_2 + \mathrm{Lip}(\phi) \|\mathbf{f}(\mu) - \mathbf{f}(\mu')\|_2 \\
&\le \mathrm{Lip}(D\phi) \cdot (\|\phi\|_{\infty} + \|\mathbf{y}\|_2) \cdot \|x - x'\|_2 + \mathrm{Lip}(\phi)^2 W_1(\mu, \mu'). \quad \square
\end{align*}
$$

We note that the term $W_1(\mu, \mu')$ gives in general pessimistic estimates when it comes to quantitative convergence estimates to the mean-field limit, but this is not our goal in this paper so we adopt this bound for convenience. See Mei et al. [2019] or Chizat [2025] for the techniques to obtain tight error bounds.

We also recall a $W_1$ version of the law of large numbers for empirical distributions.

**Lemma 5.** Let $(\mathcal{X}, \mathrm{dist})$ be a Polish space (that is, a separable metric space that admits a complete metric equivalent to $\mathrm{dist}$), let $\mu \in \mathcal{P}_1(\mathcal{X})$, let $X^1, X^2, \dots$ be independent samples from $\mu$, and let $\hat{\mu}^n = \frac{1}{n} \sum_{i=1}^n \delta_{X^i}$ be the corresponding sequence of empirical measures. Then it holds $W_1(\mu, \hat{\mu}^n) \xrightarrow{a.s.} 0$ (hence, in particular $\mathbb{E}[W_1(\mu, \hat{\mu}^n)] \to 0$).

*Proof.* We provide a proof of this well-known fact for completeness. According to the law of large numbers for empirical measures [Varadarajan, 1958] $\hat{\mu}^n$ converges weakly (in duality with continuous and bounded functions) to $\mu$, almost surely. Moreover, posing $x_0 \in \mathcal{X}$ and $\psi(x) = \mathrm{dist}(x_0, x)$, we have $\int \psi(x) d\hat{\mu}^n = \frac{1}{n} \sum_{i=1}^n \psi(X^i)$ which is a sum of i.i.d. random variables with finite first moment because $\mu \in \mathcal{P}_1(\mathcal{X})$. So, by the strong law of large numbers, $\int \mathrm{dist}(x_0, x) d\hat{\mu}^n(x)$ converges to $\int \mathrm{dist}(x_0, x) d\mu(x)$ almost surely. It follows, by Villani [2008, Theorem 6.9] that $W_1(\hat{\mu}^n, \mu) \xrightarrow{n \to \infty} 0$ almost surely. $\square$

We will consider the following Polish spaces:

* The usual Euclidean space $(\mathbb{R}^p)^K, \|\cdot\|_2$;

* The space of continuous paths $(\mathcal{C}([0, t]; \mathbb{R}^p), d_t)$ where $d_t(x, x') = \sup_{s \le t} \|x(s) - x'(s)\|_2$;

* The Skorokhod space $(\mathcal{D}([0, t], \mathbb{R}^p), d'_t)$ of càdlàg functions (right-continuous functions with left limits at all points) where $d'_t$ is the Skorokhod metric [Billingsley, 2013, Chap. 3]:

$$d'_{t}(x, x') = \inf_{\lambda \in \Lambda} \sup_{s \in [0, t]} |\lambda(s) - s| + \sup_{s \in [0, t]} \|x(\lambda(s)) - x'(s)\|_{2}$$

where $\Lambda$ is the set of strictly increasing continuous maps from $[0, t]$ to $[0, t]$.

The topologies above are those that we associate to these spaces throughout.

## 5.2 Case (I) : the discrete-time jump process limit

In this section we assume that $n \to \infty$ and that $q_n \to q \in [0, 1]$ and $\tau_n \to \tau > 0$. This is the simplest limit to study, in particular because it remains discrete-in-time in the limit.
**Limit dynamics.** Consider the (discrete-time) stochastic process $(Y_k)_{k \in \mathbb{N}}$ defined recursively by $Y_0 \sim \mu_0$ and

$$ Y_{k+1} = \begin{cases} Y_k & \text{with probability } 1-q, \\ Y_k - \alpha \nabla V[\rho_k](Y_k) & \text{with probability } q, \end{cases} \quad \rho_k = \text{Law}(Y_k). \quad (9) $$

By construction, it holds

$$ \rho_0 = \mu_0, \quad \rho_{k+1} = q(\text{Id} - \alpha \nabla V[\rho_k]) \# \rho_k + (1-q)\rho_k, \forall k \ge 1. $$

Then the discrete-time jump process limit in the main theorem (Theorem 1(I)) is a consequence of the following, more precise statement:

**Proposition 6.** Let Assumption A hold and assume that $q_n \to q \in [0, 1]$, that $\tau_n \to \tau > 0$. For any $k \in \mathbb{N}$, it holds $\mathbf{E}[W_1(\hat{\mu}_k^n, \rho_k)] \xrightarrow{n \to \infty} 0$. Moreover, for any fixed $m \in \mathbb{N}^*$ and time horizon $K \in \mathbb{N}^*$, it holds

$$ (X^{1,n}, \dots, X^{m,n}) \xrightarrow{n \to \infty} (Y^1, \dots, Y^m) $$

where $(Y^1, \dots, Y^m)$ are independent copies of the $(\mathbb{R}^p)^{K+1}$-valued random variable $Y$.

*Proof.* Let us consider constant sequences $\tau_n = \tau$ and $q_n = q$, as the general case can be proved along the same lines.

**Preparation.** The update equation (3) can be decomposed, for $k \in \mathbb{N}$, as

$$ \begin{aligned} X_{k+1}^{i,n} &= X_k^{i,n} - \tau \cdot (1 + \eta_{k+1}^i) D\phi(X_k^{i,n})^\top \left( \frac{1}{n} \sum_{j=1}^n (1 + \eta_{k+1}^j) \phi(X_k^{j,n}) - \mathbf{y} \right) \\ &= X_k^{i,n} - \tau \cdot (1 + \eta_{k+1}^i) \left( \nabla V[\hat{\mu}_k^n](X_k^{i,n}) + \gamma_{k+1}^{i,n} \right) \end{aligned} $$

where we have gathered errors terms in $\gamma_k^{i,n}$ defined as

$$ \gamma_k^{i,n} = \frac{1}{n} \sum_{j=1}^{n} \eta_k^j D\phi(X_{k-1}^{i,n})^\top \phi(X_{k-1}^{j,n}). $$

We have under Assumption A

$$ \mathbf{E} \| \gamma_k^{i,n} \|_2 \le \sqrt{\mathbf{E} \| \gamma_k^{i,n} \|_2^2} = \sqrt{O\left(\frac{n\mathbf{E}[\eta^2]}{n^2}\right)} = O(\sqrt{\beta_n}) \xrightarrow{n \to \infty} 0. \quad (10) $$

**Coupling argument.** Consider $n$ independent samples from $Y$, labelled $(Y^1, \dots, Y^n)$ which we couple with $(X^{1,n}, \dots, X^{n,n})$ via the initialization and the random dropout masks as follows: $Y_0^i = X_0^{i,n} = \xi^i$ and

$$ Y_{k+1}^i = Y_k^i - \tau(1 + \eta_{k+1}^i) \nabla V[\rho_k](Y_k^i), $$

where we recall that $\rho_k$ is defined in (9). Observe that the definition of $Y_{k+1}^i$ involves the same dropout masks as that of $X_{k+1}^i$. Let us also define $\hat{\rho}_\ell^n = \frac{1}{n} \sum_{i=1}^n \delta_{Y_\ell^i}$.
By construction, it holds $\|X_0^{i,n} - Y_0^i\|_2 = 0$ and for $k \ge 1$, using Lemma 4 and the fact that $|1 + \eta| \le \frac{1}{q}$,

$$ \begin{aligned} \|X_{k+1}^{i,n} - Y_{k+1}^i\|_2 &\le \|X_k^{i,n} - Y_k^i\|_2 + \frac{\tau}{q} \|\nabla V[\hat{\mu}_k^n](X_k^{i,n}) - \nabla V[\rho_k](Y_k^i)\|_2 + \frac{\tau}{q} \|\gamma_{k+1}^{i,n}\|_2 \\ &\le (1 + C\alpha) \|X_k^{i,n} - Y_k^i\|_2 + C\alpha W_1(\hat{\mu}_k^n, \rho_k) + \alpha \|\gamma_{k+1}^{i,n}\|_2 \end{aligned} \quad (11) $$

$$ \le \sum_{\ell=0}^{k} (1 + \alpha C)^{k-\ell} (\alpha C W_1(\hat{\mu}_{\ell}^n, \rho_{\ell}) + \alpha \|\gamma_{\ell+1}^{i,n}\|_2). \quad (12) $$

Denoting $h_k^n = \frac{1}{n} \sum_{i=1}^n \|X_k^{i,n} - Y_k^i\|_2$, it holds $h_0^n = 0$, and we deduce from (11) and from the inequality $W_1(\hat{\mu}_k^n, \hat{\rho}_k^n) \le h_k^n$ that

$$ \begin{aligned} h_{k+1}^n &\le (1 + C\alpha)h_k^n + C\alpha W_1(\hat{\mu}_k^n, \hat{\rho}_k^n) + C\alpha W_1(\hat{\rho}_k^n, \rho_k) + \frac{\alpha}{n} \sum_{i=1}^n \|\gamma_{k+1}^{i,n}\|_2 \\ &\le (1 + 2C\alpha)h_k^n + E_k^n \end{aligned} $$

where $E_k^n := C\alpha W_1(\hat{\rho}_k^n, \rho_k) + \frac{\alpha}{n} \sum_{i=1}^n \|\gamma_{k+1}^{i,n}\|_2$. Applying this inequality recursively (a discrete-time Grönwall argument), we get

$$ h_k^n \le \sum_{\ell=0}^{k-1} (1 + 2C\alpha)^{k-1-\ell} E_\ell^n \le (1 + 2C\alpha)^k \sum_{\ell=0}^{k-1} E_\ell^n. $$

Since $\mathbf{E}[W_1(\hat{\rho}_k^n, \rho_k)] \to 0, \forall \ell \in \mathbb{N}$, (Lemma 5), and $\mathbf{E}[\|\gamma_{k+1}^{i,n}\|_2] \to 0$ (by (10)), we have $\mathbf{E}[E_\ell^n] \to 0$. We deduce that

$$ \mathbf{E}[h_k^n] \le (1 + 2C\alpha)^k \sum_{\ell=0}^{k-1} \mathbf{E}[E_\ell^n] \xrightarrow{n \to \infty} 0, \forall k \in \mathbb{N}. $$

It follows

$$ \mathbf{E}[W_1(\hat{\mu}_k^n, \rho_k)] \le \mathbf{E}[W_1(\hat{\mu}_k^n, \hat{\rho}_k^n)] + \mathbf{E}[W_1(\hat{\rho}_k^n, \rho_k)] \le \mathbf{E}[h_k^n] + \mathbf{E}[W_1(\hat{\rho}_k^n, \rho_k)] \to 0 $$

which is the first claim. For the second claim, we have by (12) that for $m \in \mathbb{N}^*$,

$$ \mathbf{E} \left[ \max_{k \in [1:K]} \max_{i \in [1:m]} \|X_k^{i,n} - Y_k^i\|_2 \right] \le \sum_{\ell=0}^{K-1} (1 + \alpha C)^{K-1-\ell} \mathbf{E} \left[ \alpha C W_1(\hat{\mu}_k^n, \rho_\ell) + \alpha \max_{i \in [1:m]} \|\gamma_{\ell+1}^{i,n}\|_2 \right] \xrightarrow{n \to \infty} 0. $$

It follows that $(X^{1,n}, \dots, X^{m,n})$ converges in $L^1$ to $(Y^1, \dots, Y^m)$. This implies the second claim. $\square$

## 5.3 Case (II) : the Wasserstein gradient flow limit

In this section we assume that $n \to \infty$, that $\tau_n \to 0$ and that $\alpha_n = \frac{\tau_n}{q_n} \to 0$. Because $\alpha_n$, which represents the average time period between two updates of a neuron – vanishes in the limit, we get a flow-like dynamics in the limit, with a continuous trajectory for each neuron.
### 5.3.1 Precise statement and proof

**Limit dynamics.** Consider the continuous-time stochastic process Y that satisfies $Y_0 \sim \mu_0$ and

$$dY_t = -\nabla V[\rho_t](Y_t)dt - \beta \nabla P(Y_t)dt, \quad \rho_t = \text{Law}(Y_t)$$

where we recall that $P(x) = \frac{1}{2}\|\phi(x)\|_2^2$ (hence $\nabla P(x) = D\phi(x)^\top\phi(x)$). This is a MacKean-Vlasov process of a simple type because it is noiseless and the only source of randomness is the initialization. Under Assumption A, the existence and uniqueness of strong solutions is guaranteed (see for instance [Bauer and Meyer-Brandis, 2019] which deals with much weaker spatial regularity). Also we have for any time horizon $T$, $\text{Law}(Y) \in \mathcal{P}_1(\mathcal{C}([0, T]; \mathbb{R}^p))$. It is moreover not difficult to show that $\rho = (\rho_t)_t$ belongs to $\mathcal{C}(\mathbb{R}_+; \mathcal{P}_1(\mathbb{R}^p))$ and that it is the unique solution in the weak sense<fnref n="4" /> of

$$\partial_t \rho_t = \nabla \cdot (\rho_t \nabla [V[\rho_t] + \beta P]), \quad \rho_0 = \mu_0.$$

If $\mu_0 \in \mathcal{P}_2(\mathbb{R}^p)$ (that is, if $\mu_0$ has finite second moments), then this dynamics can be interpreted as the Wasserstein ($W_2$) gradient flow [Santambrogio, 2015, Chap. 8] of the penalized risk

$$\mu \mapsto \frac{1}{2} \|\boldsymbol{f}(\mu) - \boldsymbol{y}\|_2^2 + \beta \int P(x)d\mu(x).$$

**Interpolated dynamics.** Our statement will also involve the following continuous-time process. We denote by $(\bar{X}^{i,n})$ the piecewise affine interpolation of the discrete-time dynamics, defined by

$$\bar{X}_t^{i,n} = X_k^{i,n} + \frac{t - k\tau_n}{\tau_n}(X_{k+1}^{i,n} - X_k^{i,n}), \quad \text{for } t \in [k\tau_n, (k+1)\tau_n[$.

The gradient flow limit of the main theorem (Theorem 1(II)) is then a consequence of the following, more precise statement:

**Proposition 7.** Assume that $\tau_n \to 0$, that $\alpha_n = \frac{\tau_n}{q_n} \to 0$, and Assumption A. Then for any $t \ge 0$, it holds $\mathbf{E}[W_1(\hat{\mu}_{[t/\tau_n]}^n, \rho_t)] \xrightarrow{n \to \infty} 0$. Moreover, for any fixed $m \in \mathbb{N}^*$ and time horizon $T > 0$, it holds (in $\mathcal{P}(\mathcal{C}([0, T]; \mathbb{R}^p)^m)$)

$$(\bar{X}^{1,n}, \dots, \bar{X}^{m,n}) \xrightarrow{\text{distr.}} (Y^1, \dots, Y^m)$$

where $(Y^1, \dots, Y^m)$ are independent copies of $Y$.

**Proof.** In order to better highlight the structure of the argument, we obtain this result as a corollary of a general and stronger mean-field stochastic approximation result, Proposition 8, stated and proved the next subsection. Let us here show that the assumptions of that proposition are satisfied, that is, that we can indeed decompose the update as in (13), with $\gamma_k^{i,n}$ a Martingale difference sequence (MDS) satisfying $\max_{k \le T/\tau_n, i \le n} \mathbb{E}\|\gamma_k^{i,n}\|_2^2 = o(\tau_n)$.

By suitably decomposing the update of (3), which we recall is

$$X_{k+1}^{i,n} = X_k^{i,n} - \tau_n \cdot (1 + \eta_{k+1}^i) D\phi(X_k^{i,n})^\top \left( \frac{1}{n} \sum_{j=1}^n (1 + \eta_{k+1}^j)^\top \phi(X_k^{j,n}) - \boldsymbol{y} \right),$$

<footnote n="4">⁴By “weak solution”, we mean that $t \mapsto \rho_t$ is absolutely continuous for $W_1$ and that for any $\psi \in \mathcal{C}^1(\mathbb{R}^p)$ and a.e. $t \in \mathbb{R}_+$, it holds $\frac{d}{dt} \int \psi d\rho_t = - \int \nabla \psi \cdot \nabla [V[\rho_t] + \beta P] d\rho_t$, see e.g. [Santambrogio, 2015, Sec. 4.1.2].</footnote>
we get $X_{k+1}^{i,n} = X_k^{i,n} + \tau_n b_n(\hat{\mu}_k^n, X_k^{i,n}) + \gamma_{k+1}^{i,n}$ with $b_n$ defined by

$$ 
\begin{aligned}
b_n(\mu, x) &= -D\phi(x)^{\top}(\mathbf{f}(\mu) - \mathbf{y}) - \frac{\mathbf{E}[(\eta_{k+1}^i)^2]}{n} D\phi(x)^{\top}\phi(x) \\
&= -D\phi(x)^{\top}(\mathbf{f}(\mu) - \mathbf{y}) - \frac{(1 - q_n)}{nq_n} D\phi(x)^{\top}\phi(x)
\end{aligned}
$$

and

$$ 
\begin{aligned}
\gamma_{k+1}^{i,n} &= -\underbrace{\tau_n \eta_{k+1}^i D\phi(X_k^{i,n})^{\top} (\mathbf{f}(\hat{\mu}_k^n) - \mathbf{y})}_{(I)} - \underbrace{\tau_n \frac{(\eta_{k+1}^i)^2 - \mathbf{E}(\eta_{k+1}^i)^2}{n} D\phi(X_k^{i,n})^{\top} \phi(X_k^{i,n})}_{(II)} \\
&\quad - \underbrace{\tau_n \frac{n}{\eta_{k+1}^i} \sum_{j=1}^n \eta_{k+1}^j D\phi(X_k^{i,n})^{\top} \phi(X_k^{j,n})}_{(III)} - \underbrace{\tau_n \frac{(\eta_{k+1}^i)}{n} \sum_{\substack{j=1 \\ j \neq i}}^n \eta_{k+1}^j D\phi(X_k^{i,n})^{\top} \phi(X_k^{j,n})}_{(IV)}.
\end{aligned}
$$

* Under Assumption A, it is clear that $(b_n)$ converges uniformly to $b : (\mu, x) \mapsto -D\phi(x)^{\top}(\mathbf{f}(\mu) - \mathbf{y}) - \beta D\phi(x)^{\top}\phi(x) = -\nabla V[\mu](x) - \beta \nabla P(x)$ and that the latter satisfies the required regularity assumptions by Lemma 4.

* By construction, $(\gamma_k^{i,n})_k$ is a $(\mathbb{R}^p)^n$-valued MDS. It remains to check that it satisfies $\max_{k \le T/\tau_n, i \le n} \|\gamma_k^{i,n}\|_2^2 = o(\tau_n)$. Since by Cauchy-Schwartz's inequality, a family of centered random variables $(A_i)_{i=1}^k$ satisfies $\mathbf{E}[(\sum_{i=1}^k A_i)^2] \le \mathbf{E}[k \sum_{i=1}^k A_i^2] = k \sum_{i=1}^k \mathbf{E}[A_i^2]$, it is sufficient to check this property for the terms (I) to (IV) separately. Under Assumption A it holds

1. $\mathbf{E}\|(I)\|_2^2 = O(\mathbf{E}[(\tau_n \eta^2)]) = O(\tau_n^2 \frac{1-q_n}{q_n}) = O(\alpha_n \tau_n) = o(\tau_n)$ since $\alpha_n \to 0$;

2. $\mathbf{E}\|(II)\|_2^2 \le O(\tau_n^2 \frac{1}{n^2} \mathbf{E}[(\eta^2 - \mathbf{E}\eta^2)^2])$. Note that $\mathbf{E}\eta^4 = O(1/q_n^3)$ and $\mathbf{E}\eta^2 = O(1/q_n)$ so $\mathbf{E}[(\eta^2 - \mathbf{E}\eta^2)^2] = \mathbf{E}[\eta^4] - (\mathbf{E}[\eta^2])^2 = O(1/q_n^3)$. Overall $\mathbf{E}\|(II)\|_2^2 = O(\tau_n^2 \frac{1}{n^2 q_n^3}) = O(\tau_n \alpha_n \beta_n^2) = o(\tau_n)$, again because $\alpha_n \to 0$;

3. $\mathbf{E}\|(III)\|_2^2 = O(\tau_n^2 \frac{1}{n^2} \mathbf{E}(\sum_j \eta^j)^2) = O(\tau_n^2 \sum_j \mathbf{E}(\eta_j)^2) = O(\tau_n^2 \frac{n}{q_n}) = O(\tau_n \alpha_n) = o(\tau_n)$;

4. $\mathbf{E}\|(IV)\|_2^2 = O(\tau_n^2 \frac{1}{n^2} \mathbf{E}[(\eta_i \sum_{j \neq i} \eta^j)^2]) = O(\tau_n^2 \frac{1}{n^2} \mathbf{E}[(\eta_i)^2] \mathbf{E}[(\sum_{j \neq i} \eta^j)^2]) = O(\tau_n^2 \frac{n}{n^2} (n-1)(\mathbf{E}\eta^2)^2) = O(\tau_n^2 \frac{1}{nq_n^2}) = O(\tau_n \alpha_n \beta_n) = o(\tau_n)$.

Since all these terms are $o(\tau_n)$, this proves that Proposition 8 applies, and thus concludes the proof. $\square$

### 5.3.2 A mean-field stochastic approximation result

In this section, we prove the general mean-field stochastic approximation result used in the previous section. Consider a sequence of step-sizes $(\tau_n)$ going to 0 and a measure-dependent drift vector field $b : \mathcal{P}(\mathbb{R}^p) \times \mathbb{R}^p \to \mathbb{R}^p$ that is uniformly bounded and such that there exists $C > 0$ such that

$$ \|b(\mu, x) - b(\mu', x')\|_2 \le C(W_1(\mu, \mu') + \|x - x'\|_2), \quad \forall \mu, \mu' \in \mathcal{P}_1(\mathbb{R}^p), x, x' \in \mathbb{R}^p. $$

Consider also a sequence of velocity fields $b_n$ that converges uniformly on $\mathcal{P}_1(\mathbb{R}^p) \times \mathbb{R}^p$ to $b$. Let $\mu_0 \in \mathcal{P}(\mathbb{R}^p)$ and let $\xi^1, \xi^2, \dots$ a sequence of i.i.d. samples from $\mu_0$. Consider $\mathbb{R}^p$-valued martingale difference sequences (MDS) $(\gamma_k^{i,n})_{k \in \mathbb{N}^*}$ for $n \in \mathbb{N}^*$ and $i \in [1:n]$. This means that we assume for all $n \in \mathbb{N}^*$, $i \in [1:n]$, and $k \in \mathbb{N}^*$:
1. $\mathbf{E}[\|\gamma_k^{i,n}\|_2] < +\infty$, and

2. $\mathbf{E}[\gamma_k^{i,n}|\mathcal{F}_{k-1}^n] = 0$

where the filtration is $\mathcal{F}_k^n = \sigma(\{\xi^i\}_{i=1}^n, \{\gamma_1^{i,n}\}_{i=1}^n, \dots, \{\gamma_k^{i,n}\}_{i=1}^n)$ and $\mathcal{F}_0^n = \sigma(\{\xi^i\}_{i=1}^n)$.

**Discrete and limit dynamics.** We consider the sequence of $\mathbb{R}^p$-valued stochastic processes $(X^{1,n}, \dots, X^{n,n})$ indexed by $n$ and defined for all $n \in \mathbb{N}^*, i \in [1:n]$ and $k \in \mathbb{N}^*$ by

$$X_{k+1}^{i,n} = X_k^{i,n} + \tau_n b_n(\hat{\mu}_k^n, X_k^{i,n}) + \gamma_{k+1}^{i,n}, \quad X_0^{i,n} = \xi^i, \quad \hat{\mu}_k^n = \frac{1}{n} \sum_{i=1}^n \delta_{X_k^{i,n}}. \quad (13)$$

Our goal is to prove, under suitable assumptions, that this process converges in some sense to the continuous-time $\mathbb{R}^p$-valued stochastic process solving

$$dY_t = b(\rho_t, Y_t)dt, \quad Y_0 \sim \mu_0, \quad \rho_t = \text{Law}(Y_t).$$

The classical theory of MacKean-Vlasov equations guarantees that there is a unique strong solution $Y$ for this system with $\rho = (\rho_t)_t \in \mathcal{C}(\mathbb{R}_+; \mathcal{P}_1(\mathbb{R}^p))$. Moreover, for any time horizon $T > 0$, we have $\text{Law}(Y) \in \mathcal{P}_1(\mathcal{C}([0, T]; \mathbb{R}^p))$ (e.g. [Bauer and Meyer-Brandis, 2019]).

**Interpolated dynamics.** In order to formalize the sense in which $(X^{i,n})$ converges to $Y$, let us map the discrete-time processes $(X^{i,n})$ into continuous-time ones via a piecewise affine interpolation of the trajectories. For $k \in \mathbb{N}$, let $t_k = \tau_n k$, and for $t \in [t_k, t_{k+1}]$ let

$$\begin{aligned} \bar{X}_t^{i,n} &= X_k^{i,n} + \frac{t-t_k}{\tau_n}(X_{k+1}^{i,n} - X_k^{i,n}) \\ &= X_k^{i,n} + (t-t_k)b_n(\hat{\mu}_k, X_k^{i,n}) + \frac{(t-t_k)}{\tau_n}\gamma_{k+1}^{i,n} \end{aligned}$$

Finally, we define $\bar{\mu}_t^n = \frac{1}{n} \sum_{i=1}^n \delta_{\bar{X}_t^{i,n}}$ which interpolates the sequence of discrete measures $(\hat{\mu}_k^n)_k$. The next result guarantees that the discrete-time finite-particle dynamics converges to the continuous-time mean-field dynamics, provided the MDS terms are small enough.

**Proposition 8.** Fix a time horizon $T > 0$. As $n \to \infty$, assume that $\max_{k \le T/\tau_n, i \le n} \mathbb{E}\|\gamma_k^{i,n}\|_2^2 = o(\tau_n)$. Then $\mathbf{E}[\sup_{t \le T} W_1(\bar{\mu}_t^n, \mu_t)] \to 0$. Moreover, for any fixed $m \in \mathbb{N}^*$, it holds

$$(\bar{X}^{1,n}, \dots, \bar{X}^{m,n}) \xrightarrow{n \to \infty} (Y^1, \dots, Y^m),$$

where $(Y^1, \dots, Y^m)$ is a family of $m$ independent copies of $Y$.

*Proof.* The piecewise affine path $\bar{X}^{i,n}$ is differentiable over each interval $]t_k, t_{k+1}]$ with

$$\begin{aligned} \frac{d}{dt} \bar{X}_t^{i,n} &= b_n(\hat{\mu}_k^n, X_k^{i,n}) + \frac{1}{\tau_n} \gamma_{k+1}^{i,n} \\ &= b(\bar{\mu}_t^n, \bar{X}_t^{i,n}) + (b(\hat{\mu}_k^n, X_k^{i,n}) - b(\bar{\mu}_t^n, \bar{X}_t^{i,n})) + (b_n(\hat{\mu}_k^n, X_k^{i,n}) - b(\hat{\mu}_k^n, X_k^{i,n})) + \frac{1}{\tau_n} \gamma_{k+1}^{i,n}. \end{aligned}$$

Now consider $(Y^1, Y^2, \dots)$ a sequence of independent copies of $Y$ which is such that $Y_0^i = \xi^i$ and let us define

$$e_t^{i,n} = \bar{X}_t^{i,n} - Y_t^i, \quad h_t^n = \frac{1}{n} \sum_{i=1}^n \|e_t^{i,n}\|_2 = \frac{1}{n} \sum_{i=1}^n \|\bar{X}_t^{i,n} - Y_t^i\|_2. \quad (14)$$
We recall that $\rho_t = \text{Law}(Y_t)$ and $\hat{\rho}_t^n = \frac{1}{n} \sum_{i=1}^n \delta_{Y_t^i}$. Clearly, $t \mapsto e_t^{i,n}$ is differentiable over each interval $]t_k, t_{k+1}]$ with

$$ \frac{d}{dt} e_t^{i,n} = \underbrace{(b_n(\hat{\mu}_k^n, X_k^{i,n}) - b(\hat{\mu}_k^n, X_k^{i,n}))}_{o(1)} + \underbrace{(b(\hat{\mu}_k, X_k^{i,n}) - b(\bar{\mu}_t^n, \bar{X}_t^{i,n}))}_{I_t^{i,n}} + \underbrace{(b(\bar{\mu}_t^n, \bar{X}_t^{i,n}) - b(\rho_t, Y_t^i))}_{J_t^{i,n}} + \frac{1}{\tau_n} \gamma_{k+1}^{i,n}. $$

Since $e_t^{i,n}$ is absolutely continuous and $e_0^{i,n} = 0$, it holds $e_t^{i,n} = \int_0^t \left(\frac{d}{ds} e_s^{i,n}\right) ds$ and thus

$$ \|e_t^{i,n}\|_2 \le \left\| \int_0^t \frac{d}{ds} e_s^{i,n} ds \right\|_2 \le o(1) + \int_0^t (\|I_s^{i,n}\|_2 + \|J_s^{i,n}\|_2) ds + \sup_{s \le t} \|M_s^{i,n}\|_2 \quad (15) $$

where we have defined $M_t^{i,n} = \sum_{k=1}^{\lfloor t/\tau_n \rfloor+1} \gamma_k^{i,n}$. Note that we make these computations over a bounded time domain $[0, T]$, hence we write $to(1) = o(1)$. The last term comes from bounding the integral of the martingale term which is

$$ \sum_{k=0}^{\lfloor t/\tau_n \rfloor - 1} \int_{t_k}^{t_{k+1}} \frac{1}{\tau_n} \gamma_{k+1}^{i,n} ds + \int_{\lfloor t/\tau_n \rfloor \tau_n}^{t} \frac{1}{\tau_n} \gamma_{\lfloor t/\tau_n \rfloor + 1}^{i,n} ds \in [M_{t-\tau_n}^{i,n}, M_t^{i,n}] $$

and the last inclusion is understood as inclusion in a line segment in $\mathbb{R}^p$. Let us now derive estimates for $I_s^{i,n}$ and $J_s^{i,n}$. It holds for $t \in [t_k, t_{k+1}]$

$$ \|I_t^{i,n}\|_2 = \|b(\hat{\mu}_k^n, X_k^{i,n}) - b(\bar{\mu}_t^n, \bar{X}_t^{i,n})\|_2 \le C(W_1(\hat{\mu}_k^n, \bar{\mu}_t^n) + \|X_k^{i,n} - \bar{X}_t^{i,n}\|_2) $$

and using our assumption that $b$ is uniformly bounded (and therefore $b_n$ as well),

$$ \|X_k^{i,n} - \bar{X}_t^{i,n}\|_2 \le \tau_n \|b_n(\hat{\mu}_k^n, X_k^{i,n})\|_2 + \|\gamma_{k+1}^{i,n}\|_2 = O(\tau_n) + \|\gamma_{k+1}^{i,n}\|_2 $$

and

$$ W_1(\hat{\mu}_k^n, \bar{\mu}_t^n) \le \frac{1}{n} \sum_{i=1}^n \|X_k^{i,n} - \bar{X}_t^{i,n}\|_2 \le O(\tau_n) + \frac{1}{n} \sum_{i=1}^n \|\gamma_{k+1}^{i,n}\|_2. $$

Plugging these in the previous estimate yields

$$ \|I_t^{i,n}\|_2 \le C(O(\tau_n) + \|\gamma_{k+1}^{i,n}\|_2 + \frac{1}{n} \sum_{j=1}^n \|\gamma_{k+1}^{j,n}\|_2) $$

therefore

$$ \frac{1}{n} \sum_{i=1}^{n} \int_{0}^{t} \|I_{s}^{i,n}\|_{2} ds = O(\tau_{n}) + 2C \sum_{k=1}^{\lfloor t/\tau_{n} \rfloor+1} \frac{1}{n} \sum_{i=1}^{n} \tau_{n} \|\gamma_{k}^{i,n}\|_{2}. $$

Regarding the second error term $J_t^{i,n}$, we have

$$ \begin{aligned} \|J_t^{i,n}\|_2 &= \|b(\bar{\mu}_t^n, \bar{X}_t^{i,n}) - b(\rho_t, Y_t^i)\|_2 \\ &\le C(W_1(\bar{\mu}_t^n, \bar{\rho}_t^n) + W_1(\bar{\rho}_t^n, \rho_t) + \|\bar{X}_t^{i,n} - Y_t^i\|_2) \\ &\le C(h_t^n + W_1(\bar{\rho}_t^n, \rho_t) + \|\bar{X}_t^{i,n} - Y_t^i\|_2) \end{aligned} $$
where we have used the quantity $h_t^n$ defined in (14). Returning to (15) with the estimates collected so far yields

$$ \|e_t^{i,n}\|_2 \le o(1) + C \int_0^t \|e_s^{i,n}\|_2 ds + C \int_0^t h_s^n ds + C \int_0^t W_1(\hat{\rho}_s^n, \rho_s) ds \\ + O(\tau_n) + C \sum_{k=1}^{\lfloor t/\tau_n \rfloor+1} \tau_n \left( \|\gamma_k^{i,n}\|_2 + \frac{1}{n} \sum_{j=1}^n \|\gamma_k^{j,n}\|_2 \right) + \sup_{s \le t} \|M_s^{i,n}\|_2. \quad (16) $$

By averaging over all particles we get for any $t \in [0, T]$

$$ h_t^n \le 2C \int_0^t h_s^n ds + g_t^n $$

where

$$ g_t^n = o(1) + O(\tau_n) + C \int_0^t W_1(\hat{\rho}_s^n, \rho_s) ds \\ + 2C \sum_{k=1}^{\lfloor t/\tau_n \rfloor+1} \frac{1}{n} \sum_{i=1}^n \tau_n \|\gamma_k^{i,n}\|_2 + \frac{1}{n} \sum_{i=1}^n \sup_{s \le t} \|M_s^{i,n}\|_2. $$

Here we may choose $t \mapsto g_t^n$ to be nonnegative and nondecreasing. Therefore by Grönwall's lemma integral form (recalled in Lemma 10) it holds

$$ \sup_{t \le T} h_t^n \le e^{2CT} g_T^n \quad (17) $$

We now show that $\mathbf{E}[g_T^n] \to 0$. For this purpose, let investigate the terms involving $W_1(\hat{\rho}_s^n, \rho_s)$, $\gamma_k^{i,n}$ and $M_s^{i,n}$ in its definition. By Lemma 5, it holds $\mathbf{E}[W_1(\frac{1}{n}\sum_{i=1}^n \delta_{Y^i}, \text{Law}(Y))] \to 0$ in $\mathcal{P}_1(\mathcal{C}([0, T]; \mathbb{R}^p))$ endowed with the distance $\|x - x'\|_T = \sup_{s \le T} \|x(s) - x'(s)\|_2$. It follows (using [Lacker, 2018, Eq. (3.7)] to upper-bound the supremum of marginal $W_1$ distances by the $W_1$ distance in path-space):

$$ \mathbf{E} \left[ \int_0^T W_1(\hat{\rho}_s^n, \rho_s) ds \right] \le T \mathbf{E} \left[ \sup_{s \le T} W_1(\hat{\rho}_s^n, \rho_s) \right] \le T \mathbf{E} \left[ W_1 \left( \frac{1}{n} \sum_{i=1}^n \delta_{Y^i}, \text{Law}(Y) \right) \right] \underset{n \to \infty}{\longrightarrow} 0. \quad (18) $$

For the martingale term, we have by Doob's inequality (recalled in Lemma 9 below)

$$ \begin{aligned} \mathbf{E} \left[ \max_{s \le T} \|M_s^{i,n}\|_2 \right] &\stackrel{\text{Jensen}}{\le} \left( \mathbf{E} \left[ \max_{s \le T} \|M_s^{i,n}\|_2^2 \right] \right)^{\frac{1}{2}} \\ &\stackrel{\text{Lemma 9}}{\le} 2 \left( \mathbf{E} \left[ \sum_{1 \le k \le \lfloor T/\tau_n \rfloor + 1} \|\gamma_k^{i,n}\|_2^2 \right] \right)^{\frac{1}{2}} \\ &\stackrel{\text{assumpt.}}{=} 2 \left( \sum_{1 \le k \le \lfloor T/\tau_n \rfloor + 1} o(\tau_n) \right)^{\frac{1}{2}} = o(1). \end{aligned} $$

Finally, we have

$$ \mathbf{E} \left[ \sum_{k=1}^{\lfloor T/\tau_n \rfloor+1} \tau_n \|\gamma_k^{i,n}\|_2 \right] \le \sum_{k=1}^{\lfloor T/\tau_n \rfloor+1} \tau_n \left( \mathbf{E} \left[ \|\gamma_k^{i,n}\|_2^2 \right] \right)^{\frac{1}{2}} \le \frac{T}{\tau_n} \tau_n \sqrt{o(\tau_n)} = o(\tau_n^{1/2}). $$
All in all, and using our assumption that $\tau_n = o(1)$, we have shown that $\mathbf{E}[g_T^n] \to 0$. Plugging into (17), this gives

$$ \mathbf{E}[\sup_{s \le T} h_s^n] \le e^{2C} \mathbf{E}[g_T^n] \xrightarrow{n \to \infty} 0 $$

By the triangle inequality and (18), it follows

$$ \mathbf{E}[\sup_{t \in [0,T]} W_1(\bar{\mu}_t^n, \rho_t)] \le \mathbf{E}[\sup_{t \in [0,T]} W_1(\bar{\mu}_t^n, \hat{\rho}_t^n) + \sup_{t \in [0,T]} W_1(\hat{\rho}_t^n, \rho_t)] = \mathbf{E}[\sup_{t \in [0,T]} h_t^n] + o(1) \xrightarrow{n \to \infty} 0. $$

This proves the first claim.

As for the second claim of convergence in distribution in path space, we have, using (16) and the error estimates that follow:

$$ \|e_t^{i,n}\|_2 \le C \int_0^t \|e_s^{i,n}\|_2 ds + \tilde{g}_t^{i,n} $$

for some nonnegative $\tilde{g}_t^{i,n}$ that satisfies $\mathbf{E}[\tilde{g}_T^{i,n}] \to 0$. Therefore, by Grönwall's lemma in integral form (Lemma 10) it holds

$$ \mathbf{E}[\sup_{t \le T} \|e_t^{i,n}\|_2] \le e^{Ct} \mathbf{E}[\tilde{g}_T^{i,n}] \xrightarrow{n \to \infty} 0. $$

For a fixed $m \in \mathbb{N}^*$, this implies

$$ \mathbf{E}[\max_{i=\lceil 1:m \rceil} \sup_{s \le T} \|\bar{X}_s^{i,n} - Y_s^i\|_2] \le \sum_{i=1}^m \mathbf{E}[\sup_{s \le T} \|e_s^{i,n}\|_2] \xrightarrow{n \to \infty} 0. $$

This shows that $(\bar{X}^{i,n}, \dots, \bar{X}^{m,n})$ converges in $L^1$ to $(Y^1, \dots, Y^m)$ and the convergence in distribution follows. $\square$

For convenience, we recall below two classical results which are used in the proof: Doob's inequality (see, e.g. [Acciaio et al., 2013, Theorem 1.1]) and Grönwall's lemma.

**Lemma 9 (Doob's inequality).** Let $(M_k, \mathcal{F}_k)_{k \ge 1}$ be a real-valued martingale. Then for any $k \in \mathbb{N}^*$,

$$ \mathbf{E}[\max_{1 \le l \le k} \|M_l\|_2^2] \le 4\mathbf{E}\|M_k\|_2^2 = 4\mathbf{E}\left[\sum_{l=2}^k \|M_l - M_{l-1}\|_2^2\right]. $$

**Lemma 10 (Grönwall lemma (integral form)).** Suppose that $h$ and $g$ are nonnegative functions on $[0, T]$ with $h$ continuous and $g$ nondecreasing, and $C > 0$. If for all $t \in [0, T]$, $h(t) \le g(t) + C \int_0^t h(s)ds$ then for all $t \in [0, T]$, $h(t) \le e^{Ct}\alpha(t)$.

## 5.4 Case (III) : the continuous-time jump process limit

In this section we assume that $\tau_n \to 0$, that $\alpha_n \to \alpha > 0$ and that $\beta_n \to 0$. Because $\alpha_n$ (representing the average elapsed time between two updates of a given neuron) does not converge to zero, we get a continuous-time jump process in the limit.

### 5.4.1 The limit process: definition and well-posedness

Let $\Delta T_1, \Delta T_2, \dots$ be independent exponential random variables of expectation $\alpha$ (usually denoted $\text{Exp}(1/\alpha)$ where $1/\alpha$ is the rate), and let

$$ T_0 = 0, \quad T_j = \sum_{\ell=1}^{j} \Delta T_{\ell}, \quad N_t = \max\{j; T_j \le t\}. $$
Remark that $T_j$ follows a Gamma(j, 1/α) distribution and $N_t$ follows a Poisson(t/α) distribution (and more precisely, $t \mapsto N_t$ is a Poisson process). We consider Y the $\mathbb{R}^p$-valued continuous-time stochastic process that solves $Y_0 \sim \mu_0$ and

$$Y_t = Y_0 - \alpha \sum_{\ell=1}^{N_t} \nabla V[\rho_{(T_\ell)^-}](Y_{(T_\ell)^-}), \quad \rho_t = \text{Law}(Y_t) \quad (19)$$

where $Y_{t-}$ denotes the left limit of Y at $t$ and analogously for $\rho_{(T_\ell)^-}$. Here $T_j$ represents the time at which the particle makes its j-th jump and $N_t$ represents the total number of jumps before time $t$. It is not obvious at first that this self-referential equation has a solution and whether the solution is unique, so we first study its well-posedness. We recall the notation $\mathcal{D}(\mathbb{R}_+; \mathbb{R}^p)$ for càdlag paths in $\mathbb{R}^p$.

**Proposition 11.** There exists a unique couple $(Y, \rho)$ with $Y$ a $\mathcal{D}(\mathbb{R}_+; \mathbb{R}^p)$-valued stochastic process and $\rho \in \mathcal{C}(\mathbb{R}_+; \mathcal{P}_1(\mathbb{R}^p))$ that solves (19). For any $T > 0$, it holds $\text{Law}(Y) \in \mathcal{P}_1(\mathcal{D}([0, T]; \mathbb{R}^p))$. Moreover, there exists $L > 0$ such that for all $0 \le s, t \le T$, $\mathbf{E}[\|Y_t - Y_s\|_2] \le L\|t - s\|_2$. In particular, $t \mapsto \rho_t$ is L-Lipschitz continuous for $W_1$.

*Proof.* For $t \ge 0$, consider the space $S_t = \mathcal{C}([0, t]; \mathcal{P}_1(\mathbb{R}^p))$ of continuous functions from $[0, t]$ to the Wasserstein-1 space $\mathcal{P}_1(\mathbb{R}^p)$. Endowed with the distance $d_t(\mu, \mu') = \sup_{s \in [0, t]} W_1(\mu_s, \mu'_s)$ the space $S_t$ is a complete metric space (because $(\mathcal{P}_1(\mathbb{R}^p), W_1)$ itself is complete).

Given $t > 0$ and the random variables $(\Delta T_k)$ and $Y_0$, notice that for any fixed $\mu \in S_t$, there is easily existence of $(Y^\mu)$, the unique $\mathcal{D}([0, t]; \mathbb{R}^p)$-valued stochastic process that is a strong solution of

$$Y_t^\mu = Y_0 - \alpha \sum_{\ell=1}^{N_t} \nabla V[\mu_{T_\ell^{-}}](Y_{T_\ell^{-}}^\mu) = Y_0 - \alpha \sum_{\ell=1}^{N_t} \nabla V[\mu_{T_\ell}](Y_{T_\ell^{-1}}^\mu)$$

where we have replaced $\mu_{T_\ell^{-}}$ by $\mu_{T_\ell}$ in the last equation using the continuity of $s \mapsto \nabla V[\mu_s]$ for $\mu \in S_T$, and we have replaced $Y_{T_\ell^{-}}^\mu$ by $Y_{T_\ell^{-1}}^\mu$ which are equal, since $Y^\mu$ is constant on the intervals $[T_i, T_{i+1}]$. The map $\Phi: S_t \to S_t$ defined by $\Phi(\mu) = (\text{Law}(Y_s^\mu))_{s \in [0, t]}$ is well-posed since, for $s \le t$,

$$\mathbf{E}\|Y_t^\mu - Y_s^\mu\|_2 = \alpha \mathbf{E} \left\| \sum_{\ell=N_s+1}^{N_t} \nabla V[\mu_{T_\ell}](Y_{T_\ell-1}^\mu) \right\|_2 = O(1) \mathbf{E}[|N_s - N_t|] = O(1)|t-s|. \quad (20)$$

We are looking to prove existence and uniqueness of a fixed point for $\Phi$, at least if $t$ is small enough.

Let $\mu, \nu \in S_t$. Our first step is to control the distance between the corresponding strong solutions $Y^\mu$ and $Y^\nu$ in terms of the distance between $\mu$ and $\nu$. On the one hand, it holds

$$\|Y_t^\mu - Y_t^\nu\|_2 = \alpha \left\| \sum_{\ell=1}^{N_t} \nabla V[\mu_{T_\ell}](Y_{T_{\ell-1}}^\mu) - \sum_{\ell=1}^{N_t} \nabla V[\nu_{T_\ell}](Y_{T_{\ell-1}}^\nu) \right\|_2 \\ \le C' \sum_{\ell=1}^{N_t} \left( W_1(\mu_{T_\ell}, \nu_{T_\ell}) + \|Y_{T_{\ell-1}}^\mu - Y_{T_{\ell-1}}^\nu\|_2 \right)$$

where $C' = \alpha C$ and $C$ is the constant from (8).
On the other hand, we have by recursion that

$$
\begin{aligned}
\|Y_{T_\ell}^\mu - Y_{T_\ell}^\nu\|_2 &\le \|Y_{T_{\ell-1}}^\mu - Y_{T_{\ell-1}}^\nu\|_2 + \alpha \|\nabla V[\mu_{T_\ell}](Y_{T_{\ell-1}}^\mu) - \nabla V[\nu_{T_\ell}](Y_{T_{\ell-1}}^\mu)\|_2 \\
&\le \|Y_{T_{\ell-1}}^\mu - Y_{T_{\ell-1}}^\nu\|_2 + C' (W_1(\mu_{T_\ell}, \nu_{T_\ell}) + \|Y_{T_{\ell-1}}^\mu - Y_{T_{\ell-1}}^\nu\|_2) \\
&\le \sum_{k=1}^\ell C'(1+C')^{\ell-k} W_1(\mu_{T_k}, \nu_{T_k}).
\end{aligned}
$$

Since $\|Y_t^\mu - Y_t^\nu\|_2 = \|Y_{T_{N_t}}^\mu - Y_{T_{N_t}}^\nu\|_2$ because $Y^\mu$ and $Y^\nu$ are constant on the intervals $[T_i, T_{i+1}]$, this recursion gives

$$
\sup_{s \in [0,t]} \|Y_s^\mu - Y_s^\nu\|_2 \le C'(1+C')^{N_t-1} \sum_{k=1}^{N_t} W_1(\mu_{T_k}, \nu_{T_k}).
$$

Now by the law of total expectation, it holds

$$
\mathbf{E}[\sup_{s \in [0,t]} \|Y_s^\mu - Y_s^\nu\|_2] \le C' \sum_{N=0}^{+\infty} \mathbf{P}[N_t = N](1+C')^{N-1} \mathbf{E}\left[\sum_{k=1}^{N} W_1(\mu_{T_k}, \nu_{T_k}) \mid N_t = N\right].
$$

We recall that $N_t$ follows a Poisson distribution with intensity $t/\alpha$. By the memoryless property of the Poisson distribution, conditioned on $N_t = N$, the random variables $(T_1, \dots, T_N)$ are just a re-ordering (in increasing order) of $N$ independent uniform samples on $[0, t]$. Hence

$$
\begin{aligned}
\mathbf{E}[\sup_{s \in [0,t]} \|Y_s^\mu - Y_s^\nu\|_2] &\le C' \sum_{N=0}^{+\infty} \mathbf{P}[N_t = N](1+C')^{N-1} \frac{N}{t} \int_0^t W_1(\mu_s, \nu_s) \mathrm{d}s \\
&= \frac{C'}{t} e^{-t/\alpha} \left( \sum_{N=0}^{+\infty} \frac{(t/\alpha)^N}{N!} (1+C')^{N-1} N \right) \int_0^t W_1(\mu_s, \nu_s) \mathrm{d}s \\
&= \frac{C't}{\alpha t} e^{tC'/\alpha} \int_0^t W_1(\mu_s, \nu_s) \mathrm{d}s
\end{aligned}
$$

Since $(Y_t^\mu, Y_t^\nu)$ is a coupling of $\Phi(\mu)_t$ and $\Phi(\nu)_t$, it follows

$$
d_t(\Phi(\mu), \Phi(\nu)) \le \frac{C'}{\alpha} e^{tC'/\alpha} \int_0^t W_1(\mu_s, \nu_s) \mathrm{d}s \le \frac{C't}{\alpha} e^{tC'/\alpha} d_t(\mu, \nu).
$$

Hence by taking $t_0$ small enough, the map $\Phi$ is a contraction on $S_{t_0}$. By Banach's fixed point theorem, $\Phi$ has a unique fixed-point, which proves existence and uniqueness of solutions to (19) on a time interval $[0, t_0]$. Since the constant $C'$ above does not depend on the initial distribution $\mu_0$, we can repeat the same argument starting from $t_0$ and recursively, we get existence and uniqueness over any time interval $[0, t]$ for any $t > 0$ and hence, on $\mathbb{R}_+$.

Then, from (20) we deduce that there exists $L > 0$ such that for $0 \le s, t \le T$ $\mathbf{E}[\|Y_t - Y_s\|_2] \le L|t - s|$. Since $W_1(\mu_s, \mu_t) \le \mathbf{E}\|Y_s - Y_t\|_2$, this implies that $t \mapsto \mu_t$ is $L$-Lipschitz continuous for $W_1$. $\square$

5.4.2 Evolution equation in the space of probability measures

In this subsection, we study the equation satisfied by $\rho_t = \text{Law}(Y_t)$ which, as will be proven, is given by the following evolution equation

$$
\partial_t \rho_t = \frac{1}{\alpha} ((\text{Id} - \alpha \nabla V[\rho_t]) \# \rho_t - \rho_t), \qquad \rho_0 = \mu_0. \tag{21}
$$
We say that $\rho$ solves this equation if $\rho \in C(\mathbb{R}_+; \mathcal{P}_1(\mathbb{R}^p))$ and (21) is satisfied in the weak sense, that is

$$ \frac{d}{dt} \int \psi(x) d\rho_t(x) = \frac{1}{\alpha} \int \left( \psi(x - \alpha \nabla V[\rho_t](x)) - \psi(x) \right) d\rho_t(x), \quad \forall \psi \in \mathcal{C}^b(\mathbb{R}^p) $$

where $\mathcal{C}^b(\mathbb{R}^p)$ is the set of bounded and continuous functions on $\mathbb{R}^p$. Here again, the fact that this equation is well-posed is not immediate. We first prove that the law of $(Y_t)$ where $Y$ the solution to (19) is a solution – thereby proving existence of solutions – and separately, we will prove uniqueness of solutions below.

**Proposition 12.** If $(Y, \rho)$ is the unique solution to (19), then $\rho$ is a solution to (21).

*Proof.* Let $\psi \in \mathcal{C}^b(\mathbb{R}^p)$ be a test function. For any $t, \Delta t > 0$, by partitioning the expectation according to the number of jumps $\Delta N = N_{t+\Delta t} - N_t$ between $t$ and $t+\Delta t$ (which is a Poisson random variable of intensity $\Delta t/\alpha$) we have, by continuity of $t \mapsto \nabla V[\rho_t]$,

$$ \begin{aligned} \mathbf{E}[\psi(Y_{t+\Delta t})] &= \sum_{N=0}^{\infty} \mathbf{P}(\Delta N = N) \mathbf{E} \left[ \psi \left( Y_t - \alpha \sum_{\ell=1}^{N} \nabla V[\rho_{T_{N_t+\ell}}](Y_{T_{N_t+\ell}}^{-}) \right) \mid \Delta N = N \right] \\ &= e^{-\Delta t/\alpha} \mathbf{E} \left[ \psi(Y_t) \right] + \frac{\Delta t}{\alpha} e^{-\Delta t/\alpha} \mathbf{E} \left[ \psi(Y_t - \alpha \nabla V[\rho_{T_{N_t+1}}](Y_{T_{N_t+1}}^{-})) \mid \Delta N = 1 \right] \\ &\quad + \sum_{N=2}^{\infty} \frac{(\Delta t/\alpha)^N}{N!} e^{-\Delta t/\alpha} \mathbf{E} \left[ \psi \left( Y_t - \alpha \sum_{\ell=1}^{N} \nabla V[\rho_{T_{N_t+\ell}}](Y_{T_{N_t+\ell}}^{-}) \right) \mid \Delta N = N \right]. \end{aligned} $$

Using that $Y_{T_{N_t+1}}^{-} = Y_t$, it follows

$$ \begin{aligned} & \frac{\mathbf{E}[\psi(Y_{t+\Delta t})] - \mathbf{E}[\psi(Y_t)]}{\Delta t} \\ &= \frac{e^{-\Delta t/\alpha} - 1}{\Delta t} \mathbf{E}[\psi(Y_t)] + \frac{1}{\alpha} e^{-\Delta t/\alpha} \mathbf{E}[\psi(Y_t - \alpha \nabla V[\rho_{T_{N_t+1}}](Y_t)) \mid \Delta N = 1] \\ &\quad + \frac{1}{\Delta t} \sum_{N=2}^{\infty} \frac{(\Delta t/\alpha)^N}{N!} e^{-\Delta t/\alpha} \mathbf{E}[\psi(Y_t - \alpha \sum_{\ell=1}^{N} \nabla V[\rho_{T_{N_t+\ell}}](Y_{T_{N_t+\ell}}^{-})) \mid \Delta N = N]. \end{aligned} $$

We then get

$$ \frac{\mathbf{E}[\psi(Y_{t+\Delta t})] - \mathbf{E}[\psi(Y_t)]}{\Delta t} \underset{\Delta t \to 0}{\longrightarrow} \frac{-1}{\alpha} \mathbf{E}[\psi(Y_t)] + \frac{1}{\alpha} \mathbf{E}[\psi(Y_t - \alpha \nabla V[\rho_t](Y_t))] + 0, $$

by the continuity of $t \mapsto \nabla V[\rho_t]$ in the supremum norm for the second term (indeed, if $\Delta N = 1$ then $T_{N_1+1} \in [t, t+\Delta t]$), and the boundedness of $\psi$ to obtain that the last expression converges to 0. This shows that $t \mapsto \mathbf{E}[\psi(Y_t)]$ is differentiable, with derivative

$$ \frac{d}{dt} \mathbf{E}[\psi(Y_t)] = \frac{1}{\alpha} \mathbf{E}[\psi(Y_t - \alpha \nabla V[\rho_t](Y_t)) - \psi(Y_t)] $$

which, complemented with the initial condition $\mathbf{E}[\psi(Y_0)] = \int \psi d\mu_0$, shows that $\rho$ is a weak solution to (21). $\square$

Regarding uniqueness of solutions of the evolution equation (21), we will use a fixed-point argument that directly proves existence and uniqueness, independently of the stochastic process solving (19).
**Proposition 13.** Given $\mu_0 \in \mathcal{P}_1(\mathbb{R}^p)$, there exists a unique $\rho \in \mathcal{C}(\mathbb{R}_+; \mathcal{P}_1(\mathbb{R}^p))$ that solves (21). It is also the unique solution to the integral equation:

$$ \rho_t = e^{-t/\alpha}\mu_0 + \frac{1}{\alpha}\int_0^t e^{(s-t)/\alpha}(\mathrm{Id} - \alpha\nabla V[\rho_s])_{\#}\rho_s ds. $$

*Proof.* Consider the operator $\Phi : \mathcal{C}(\mathbb{R}_+; \mathcal{P}_1(\mathbb{R}^p)) \to \mathcal{C}(\mathbb{R}_+; \mathcal{P}_1(\mathbb{R}^p))$ defined as

$$ \Phi(\rho)_t = e^{-t/\alpha}\mu_0 + \frac{1}{\alpha}\int_0^t e^{(s-t)/\alpha}(\mathrm{Id} - \alpha\nabla V[\rho_s])_{\#}\rho_s ds. $$

Let us check that this operator indeed maps into $\mathcal{C}(\mathbb{R}_+; \mathcal{P}_1(\mathbb{R}^p))$. First, the right-hand side of this equation should be understood in duality with bounded continuous functions, that is $\Phi(\rho)_t$ is characterized by

$$ \int_{\mathbb{R}^p} \psi d\Phi(\rho)_t = e^{-t/\alpha} \int_{\mathbb{R}^p} \psi d\mu_0 + \frac{1}{\alpha} \int_0^t e^{(s-t)/\alpha} \int_{\mathbb{R}^p} \psi(x - \alpha \nabla V[\rho_s](x)) d\rho_s(x) ds. $$

Since the right-hand side of this equation defines a continuous linear map over $\mathcal{C}^b(\mathbb{R}^p)$ that is nonnegative and maps the constant function equal to 1 to 1, this indeed defines a probability measure. Moreover, for any $\psi \in \mathcal{C}^b(\mathbb{R}^p)$, the map $t \mapsto \int \psi\Phi(\rho)_t$ is differentiable and

$$ \begin{aligned} \frac{d}{dt} \int \psi d\Phi(\rho)_t &= -\frac{e^{-t/\alpha}}{\alpha} \int \psi d\mu_0 - \frac{1}{\alpha^2} \int_0^t e^{(s-t)/\alpha} \int \psi(x - \alpha \nabla V[\rho_s](x)) d\rho_s(x) ds \\ &\quad + \frac{1}{\alpha} \int \psi(x - \alpha \nabla V[\rho_t](x)) d\rho_t(x) \\ &= \frac{1}{\alpha} \left[ \int \psi(x - \alpha \nabla V[\rho_t](x)) d\rho_t(x) - \int \psi d\Phi(\rho_t) \right]. \end{aligned} $$

Clearly, if $\Phi(\rho) = \rho$ then $\rho$ solves (21). On the other hand assume that $\rho$ solves (21); then by the previous computation, for $\psi \in \mathcal{C}^b(\mathbb{R}^p)$ it holds for almost every $t \in \mathbb{R}_+$,

$$ \frac{d}{dt} \int \psi d\Phi(\rho)_t = \frac{d}{dt} \int \psi d\rho_t + \frac{1}{\alpha} \int \psi d\rho_t - \frac{1}{\alpha} \int \psi d\Phi(\rho)_t. $$

So if we let $h(t) = \int \psi d(\Phi(\rho)_t - \rho_t)$ this function is absolutely continuous with $h(0) = 0$ and for almost every $t \in \mathbb{R}_+$

$$ h'(t) = -\frac{1}{\alpha}h(t) $$

and so $h(t) = 0$. This implies that $\rho$ is a fixed point of $\Phi$. So $\rho$ solves (21) iff $\Phi(\rho) = \rho$. We are thus looking for fixed points of $\Phi$.

Recall that $S_t = \mathcal{C}([0, t]; \mathcal{P}_1(\mathbb{R}^p))$ is endowed with the metric $d_t(\rho, \rho') = \sup_{s \in [0, t]} W_1(\rho_s, \rho'_s)$. We will next show that $\Phi$ is a contraction of the complete metric space $(S_t, d_t)$ for $t$ small enough. By the subadditivity of $W_1$ ([Villani, 2008, Theorem 4.8]), we have

$$ W_1(\Phi(\rho)_t, \Phi(\rho')_t) \le \frac{1}{\alpha} \int_0^t e^{(s-t)/\alpha} W_1((I - \alpha \nabla V[\rho_s])_{\#}\rho_s, (I - \alpha \nabla V[\rho'_s])_{\#}\rho'_s) ds. $$

Now for all Lipschitz maps $T, T' : \mathbb{R}^p \to \mathbb{R}^p$ it is not difficult to show that

$$ \begin{aligned} W_1((\mathrm{Id} + T)_{\#}\rho, (\mathrm{Id} + T')_{\#}\rho') &\le W_1((\mathrm{Id} + T)_{\#}\rho, (\mathrm{Id} + T)_{\#}\rho') + W_1((\mathrm{Id} + T)_{\#}\rho', (\mathrm{Id} + T')_{\#}\rho') \\ &\le \mathrm{Lip}(T)W_1(\rho, \rho') + \|T - T'\|_{L^1(\rho')}. \end{aligned} $$
Coming back to our previous formula, we deduce

$$W_1(\Phi(\rho)_t, \Phi(\rho')_t) \leq \frac{1}{\alpha} \int_0^t e^{(s-t)/\alpha} (LW_1(\rho_s, \rho'_s) + \alpha \|\nabla V[\rho_s] - \nabla V[\rho'_s]\|_{\infty}) ds$$

where $L$ here refers to a bound on the Lipschitz constant of any map of the form $x \mapsto -\alpha \nabla V[\rho_s](x)$ which exists under Assumption A. Since the last term is also bounded by $CW_1(\rho_s, \rho'_s)$ by Lemma 4, one has that for some $C' > 0$,

$$W_1(\Phi(\rho)_t, \Phi(\rho')_t) \leq C' \int_0^t e^{(s-t)/\alpha} W_1(\rho_s, \rho'_s) ds.$$

This implies

$$d_t(\Phi(\rho), \Phi(\rho')) \leq C' \int_0^t e^{(s-t)/\alpha} W_1(\rho_s, \rho'_s) ds \leq C' \alpha (1 - e^{-t/\alpha}) d_t(\rho, \rho').$$

Hence for $t$ small enough, $\Phi$ is a contraction and admits a unique fixed point. Since $C'$ does not depend on $\mu_0$, we can repeat the argument and get existence and uniqueness on $\mathbb{R}_+$. $\square$

### 5.4.3 Proof of the mean-field limit

**Piecewise constant interpolation.** In order to formalize the sense in which $(X^{i,n})$ converges to $Y$, let us map the discrete-time processes $(X^{i,n})$ into continuous-time ones via a piecewise constant interpolation of the trajectories. For $k \in \mathbb{N}$, let $t_k = \tau k$, and for $t \in [t_k, t_{k+1}]$ [let

$$\bar{X}_t^{i,n} = X_k^{i,n}, \qquad \bar{\mu}_t^n = \frac{1}{n} \sum_{i=1}^n \delta_{\bar{X}_t^{i,n}}.$$

Claim (III) of the main theorem (Theorem 1) is then a consequence of the following statement:

**Proposition 14.** Assume that $\tau_n \to 0$, that $\alpha_n \to \alpha > 0$ and that $\beta_n \to 0$. For any $t \geq 0$, it holds $\mathbf{E}[W_1(\hat{\mu}_{[t/\tau_n]}, \rho_t)] \xrightarrow{n \to \infty} 0$ where $\rho$ is the unique solution to (21).

*Proof.* For this proof, we will build a coupling between the discrete and the limit dynamics and estimate the errors in an asynchronous way. This requires first to reformulate the discrete dynamics in a form analogous to the limit one.

**Step 1: Reformulation of the discrete dynamics.** The first step of this proof is identical to the beginning of the proof of Proposition 6. The update equation (3) can be decomposed, for $k \in \mathbb{N}$, as

$$X_{k+1}^{i,n} = X_k^{i,n} - \tau_n \cdot (1 + \eta_{k+1}^i) D\phi(X_k^{i,n})^\top \left( \frac{1}{n} \sum_{j=1}^n (1 + \eta_{k+1}^j) \phi(X_k^{j,n}) - \mathbf{y} \right) \\ = X_k^{i,n} - \tau_n \cdot (1 + \eta_{k+1}^i) \left( \nabla V[\hat{\mu}_k^n](X_k^{i,n}) + \gamma_{k+1}^{i,n} \right)$$

where we have gathered errors terms in $\gamma_k^{i,n}$ defined as

$$\gamma_k^{i,n} = \frac{1}{n} \sum_{j=1}^{n} \eta_k^j D\phi(X_{k-1}^{i,n})^\top \phi(X_{k-1}^{j,n})$$
We have under Assumption A

$$ \mathbf{E}\|\gamma_k^{i,n}\|_2 \le \sqrt{\mathbf{E}\|\gamma_k^{i,n}\|_2^2} = \sqrt{O\left(\frac{n\mathbf{E}[\eta^2]}{n^2}\right)} = O(\sqrt{\beta_n}) \to 0, \quad (22) $$

which suggests that these quantities indeed are error terms, as it will be confirmed in the rest of the proof.

Let $\Delta K_{\ell}^{i}$ be the number of steps between jump number $\ell$ and $\ell+1$ of the $i$-th neuron, that is if $\eta_{k}^{i} > 0$ is the $\ell$-th nonnegative term in the sequence $(\eta_{k}^{i})_{k}$, we set $\Delta K_{\ell}^{i} = \min\{k'; \eta_{k+k'}^{i} > 0\}$. By construction, the variables $\Delta K_{1}^{i}, \Delta K_{2}^{i}, \dots$ are independent $\text{Geom}(q)$ random variables, with support $\{1, 2, \dots\}$. Also let

$$ K_0^i = 0, \quad K_j^i = K_0^i + \sum_{\ell=1}^j \Delta K_\ell^i, \quad M_k^i = \max\{j; K_j^i \le k\} $$

Here $(K_j^i)$ for $j \ge 1$ is the index of the (landing of the) $j$-th jump and corresponds to the $j$-th time $(\eta_k^i)_k$ takes a positive value, while $M_k^i$ is the number of jumps landing before or on step $k$.

**Step 2: Coupling with the limit dynamics.** Now consider $(Y^1, Y^2, \dots)$ independent copies of $Y$ (the solution to (19)). Let us couple these copies with the discrete dynamics as follows: (i) we take $Y_0^i = \xi^i$ and (ii) we take the coupling $(\Delta T_j^i, \Delta K_j^i)$ given by Lemma 16. For $k \in [K_{\ell}^i : K_{\ell+1}^i - 1]$ (hence it holds $\ell = M_k^i$), we define

$$ e_k^{i,n} = X_k^{i,n} - Y_{T_{\ell}^i} $$

Since $Y^i$ is constant on the interval $[T_{\ell}^i, T_{\ell+1}^i[$, we could have replaced in this definition $Y_{T_{\ell}^i}$ by $Y_t$ for any $t$ in this interval without consequences. In words, $e_k^{i,n}$ compares the coupled processes $X^{i,n}$ and $Y^i$ after the same number of jumps. Note that, in expectation, this quantity behaves similarly as the difference between these processes at the same time since

$$ \mathbf{E}\|X_k^{i,n} - Y_{\tau_n k}^i\|_2 \le \mathbf{E}\|X_k^{i,n} - Y_{T_k^i}^{i,n}\|_2 + \mathbf{E}\|Y_{T_k^i}^{i,n} - Y_{\tau_n k}^i\|_2 \le \mathbf{E}\|e_k^{i,n}\|_2 + O(\tau_n + |\alpha_n - \alpha|) \quad (23) $$

where the estimate $\mathbf{E}\|Y_{T_k^i}^i - Y_{\tau_n k}^i\|_2 = \mathbf{E}\|Y_{T_k^i}^i - Y_{T_{\ell+1}^i}^i\|_2 = O(\tau_n + |\alpha_n - \alpha|)$ is proved in Lemma 19. However, the comparator $e_k^{i,n}$ will prove simpler to work with in the following.

If $\eta_{k+1}^i = -1$, then it holds $e_{k+1}^{i,n} = e_k^{i,n}$ (there is no jump), otherwise $\eta_{k+1}^i = (1 - q_n)/q_n$ and $k+1 = K_{\ell+1}^i$ for $\ell = M_k^i$ so it holds

$$ e_{k+1}^{i,n} = e_k^{i,n} - \frac{\tau_n}{q_n} \left( \nabla V[\hat{\mu}_k^n](X_k^{i,n}) + \gamma_{k+1}^{i,n} \right) + \alpha \nabla V[\rho_{T_{\ell+1}^i}^i](Y_{T_{\ell}^i}^i) $$

hence

$$ \|e_{k+1}^{i,n}\|_2 \le \|e_k^{i,n}\|_2 + \alpha C \left( W_1(\hat{\mu}_k^n, \rho_{T_{\ell+1}^i}) + \|X_k^{i,n} - Y_{T_{\ell}^i}^i\|_2 + \|\gamma_{k+1}^{i,n}\|_2 \right) \\ + (C + \|\gamma_{k+1}^{i,n}\|_2)|\alpha_n - \alpha|. $$

Moreover, the Wasserstein distance can be decomposed as

$$ W_1(\hat{\mu}_k^n, \rho_{T_{\ell+1}^i}) \le W_1(\hat{\mu}_k^n, \hat{\rho}_{k\tau_n}^n) + W_1(\hat{\rho}_{k\tau_n}^n, \rho_{k\tau_n}) + W_1(\rho_{k\tau_n}, \rho_{T_{\ell+1}^i}). $$
Also, by the Lipschitz property of Proposition 11, it holds

$$W_1(\rho_{k\tau_n}, \rho_{T_{\ell+1}^i}) \le L|k\tau_n - T_{\ell+1}^i| \le L \cdot |\tau_n(K_{\ell+1}^i - 1) - T_{\ell+1}^i| \le L \cdot |\tau_n K_{\ell+1}^i - T_{\ell+1}^i| + o(1).$$

Hence, still when $k+1 = K_{\ell+1}^i$, it holds (assuming $|\alpha_n - \alpha| \le \alpha C$, which holds eventually):

$$\|e_{k+1}^{i,n}\|_2 \le \|e_k^{i,n}\|_2 + \alpha C \left( W_1(\hat{\mu}_k^n, \hat{\rho}_{k\tau_n}^n) + W_1(\hat{\rho}_{k\tau_n}^n, \rho_{k\tau_n}) \right. \\ \left. \qquad + L \cdot \sup_{\ell \le M_k^i} |\tau_n K_{\ell+1}^i - T_{\ell+1}^i| + \|e_k^{i,n}\|_2 + 2\|\gamma_{k+1}^{i,n}\|_2 \right) + o(1)$$

Now, for a fixed (deterministic) $k$, accounting for the probability $q_n$ of jumps, it follows taking expectation and using (22) that

$$\mathbf{E}\|e_{k+1}^{i,n}\|_2 \le \mathbf{E}\|e_k^{i,n}\|_2 \\ + \alpha q_n C \mathbf{E} \left( W_1(\hat{\mu}_k^n, \hat{\rho}_{k\tau_n}^n) + W_1(\hat{\rho}_{k\tau_n}^n, \rho_{k\tau_n}) + L \cdot |\tau_n K_{\ell+1}^i - T_{\ell+1}^i| + \|e_k^{i,n}\|_2 \right) + o(q_n) \\ \le (1 + \alpha C q_n) \mathbf{E}\|e_k^{i,n}\|_2 + \alpha q_n C \mathbf{E}W_1(\hat{\mu}_k^n, \hat{\rho}_{k\tau_n}^n) + o(q_n)$$

where the second inequality uses Corollary 18 (that shows $\mathbf{E} \sup_{\ell \le M_k^i} |\tau_n K_{\ell+1}^i - T_{\ell+1}^i| = o(1)$) and Lemma 15 (that shows $\mathbf{E}[\sup_{s \le t} W_1(\hat{\rho}_s^n, \rho_s)] \xrightarrow{n \to \infty} 0$). We note that the term $o(q_n)$ is uniform for $k \in [0, t/\tau_n]$.

Consider now the quantity $h_k^n = \frac{1}{n} \sum_{i=1}^n \mathbf{E}\|e_k^{i,n}\|_2$ which satisfies, by (23),

$$\mathbf{E}W_1(\hat{\mu}_k^n, \hat{\rho}_{k\tau_n}^n) \le \frac{1}{n} \sum_{i=1}^n \mathbf{E}\|X_k^{i,n} - Y_{k\tau_n}^i\|_2 \le h_k^n + o(1).$$

From the previous inequality, we deduce

$$h_{k+1}^n \le (1 + 2\alpha C q_n)h_k^n + o(q_n).$$

Since $h_0^n = 0$, it follows that

$$\sup_{k \le t/\tau_n} h_k^n \le o(q_n) \sum_{i=0}^{k} (1 + 2\alpha C q_n)^i \le \frac{o(q_n)}{\tau_n} (1 + 2C\tau_n)^{t/\tau_n} = o\left(\frac{1}{\alpha_n}\right) \xrightarrow{n \to \infty} 0.$$

This shows that for any $t \ge 0$,

$$\mathbf{E}W_1(\hat{\mu}_{[t/\tau_n]}^n, \rho_t) \le \mathbf{E}W_1(\hat{\mu}_{[t/\tau_n]}^n, \hat{\rho}_{\tau_n[t/\tau_n]}^n) + \mathbf{E}W_1(\hat{\rho}_{\tau_n[t/\tau_n]}^n, \rho_{\tau_n[t/\tau_n]}) + \mathbf{E}W_1(\rho_{\tau_n[t/\tau_n]}, \rho_t) \\ \le h_{[t/\tau_n]}^n + o(1) + L|t - \tau_n[t/\tau_n]| \xrightarrow{n \to \infty} 0. \quad \square$$

### 5.4.4 Auxiliary lemmas

In the following lemma, we consider $\rho$ and $\hat{\rho}^n$ as defined in the proof of Proposition 14.

**Lemma 15.** For any $t > 0$, it holds $\mathbf{E}[\sup_{s \le t} W_1(\hat{\rho}_s^n, \rho_s)] \to 0$.

*Proof.* We have seen in (20) that $t \mapsto \rho_t$ is Lipschitz continuous for $W_1$, and the same holds for $t \mapsto \hat{\rho}_t^n$ via the same argument. Let us denote $L$ a common Lipschitz constant. For $\varepsilon > 0$, we thus have

$$\sup_{s \le t} W_1(\hat{\rho}_s^n, \rho_s) \le \varepsilon + \sup_{k \in [0: \lfloor Lt/\varepsilon \rfloor]} W_1(\hat{\rho}_{s_k}^n, \rho_{s_k})$$
where $s_k = k\epsilon/L$. Taking expectations, it follows

$$ \mathbf{E} \left[ \sup_{s \le t} W_1(\hat{\rho}_s^n, \rho_s) \right] \le \epsilon + \sum_{k=0}^{\lfloor Lt/\epsilon \rfloor} \mathbf{E} [W_1(\hat{\rho}_{s_k}^n, \rho_{s_k})] \xrightarrow{n \to \infty} \epsilon $$

by Lemma 5. Since $\epsilon$ is arbitrary, this proves the claim. □

**Lemma 16 (Jump times coupling).** Let $\tau_n, q_n > 0$ be such that $\tau_n \to 0$ and $\alpha_n = \frac{\tau_n}{q_n} \to \alpha > 0$. If $X \sim \text{Exp}(1/\alpha)$ and $Y_n/\tau_n \sim \text{Geom}(q_n)$ then there exists a sequence of couplings $(X, Y_n)$ such that $\lim_{n \to \infty} \mathbf{E}|Y_n - X| = 0$. Moreover, for any $x > 0$, these couplings satisfy $\lim_{n \to \infty} \text{ess-sup} [|Y_n - X| | X \le x] = 0$.

*Proof.* For $\tau > 0$, the random variable $Z = \lfloor X/\tau \rfloor$ has distribution $\text{Geom}(1 - e^{-\tau/\alpha})$. Indeed, for $k = 1, 2, \dots$ we have

$$ P(Z \ge k) = P(X > (k-1)\tau) = \exp(-(k-1)\tau/\alpha) = (e^{-\tau/\alpha})^{k-1} $$

Recall that a random variable $\zeta$ has distribution $\text{Geom}(\pi)$ iff for $k = 1, 2, \dots$ it holds that

$$ P(\zeta \ge k) = \sum_{j=k}^{\infty} (1-\pi)^{j-1}\pi = (1-\pi)^{k-1}\pi \sum_{j=0}^{\infty} (1-\pi)^j = (1-\pi)^{k-1}. $$

Hence, $Z \sim \text{Geom}(1 - e^{-\tau/\alpha})$. Now consider $X'_n = -X \cdot \tau_n/(\alpha \log(1 - q_n))$ and $Y_n = \tau_n[X'_n/\tau_n]$. By construction, $X'_n \sim \text{Exp}(-\log(1 - q_n)/\tau_n)$ and $Y_n/\tau_n \sim \text{Geom}(q_n)$. Moreover, we have

$$ \begin{aligned} |X - Y_n| &\le |X - X'_n| + |X'_n - \tau_n[X'_n/\tau_n]| \\ &\le \left|1 - \frac{\tau_n}{-\alpha \log(1 - q_n)}\right| |X| + \tau_n \\ &= o(1)|X| + \tau_n. \end{aligned} $$

The first claim follows by taking the expectation, the second by taking the essential-supremum conditioned on $X \le x$. □

**Lemma 17.** If $(\Delta T_l, \Delta K_l)$ follows the coupling in Lemma 16 the following inequalities hold for large enough $n$

$$ \Delta T_l \frac{\alpha_n}{(1+q_n)\alpha} \le \tau_n \Delta K_l < \Delta T_l \frac{\alpha_n}{\alpha} + \tau_n. $$

Moreover, for any $j \le N_T$ and large enough $n$ we have

$$ T_j \frac{\alpha_n}{(1+q_n)\alpha} \le \tau_n K_j < T_j \frac{\alpha_n}{\alpha} + j\tau_n. $$

In particular, for $n$ large enough, and any $k \le T/\tau_n$,

$$ M_k = \sup\{j; K_j \le k\} \le N_{2\tau_n k} = \sup\{j; T_j \le 2\tau_n k\}. $$

*Proof.* From the Taylor expansion $-\log(1 - q_n) = q_n + q_n^2/2 + O(q_n^3)$ it follows that for large enough $n$

$$ \frac{1}{\alpha(q_n + q_n^2)} \le \frac{1}{-\alpha \log(1 - q_n)} < \frac{1}{\alpha q_n}. $$
Now recalling that $\tau_n \Delta K_l = \tau_n [-\Delta T_l / (\alpha \log(1 - q_n))]$ we get

$$ \tau_n \Delta K_l \le \tau_n \left[ \frac{\Delta T_l}{\alpha_n q_n} \cdot \frac{\alpha_n}{\alpha} \right] \le \tau_n \left( \frac{\Delta T_l}{\alpha_n q_n} \cdot \frac{\alpha_n}{\alpha} + 1 \right) = \tau_n + \Delta T_l \frac{\alpha_n}{\alpha} $$

and

$$ \tau_n \Delta K_l \ge \tau_n \left[ \frac{\Delta T_l}{\alpha(q_n + q_n^2)} \right] \ge \Delta T_l \frac{\alpha_n}{(1 + q_n)\alpha} $$

For the last expression, we use

$$ M_k^n = \sup\{j; \tau_n K_j^n \le \tau_n k\} \le \sup\{j; T_j \frac{\alpha_n}{(1+q_n)\alpha} \le \tau_n k\} = N_{\tau_n k \frac{(1+q_n)\alpha}{\alpha_n}} \le N_{2\tau_n k} $$

for $n$ large enough.

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

**Corollary 18.** It holds

$$ \mathbf{E} \sup_{l \le M_{T/\tau_n}} |\tau_n K_{l+1} - T_{l+1}| = o(1). $$

*Proof.* By Lemma 17 we have

$$ |\tau_n K_{l+1} - T_{l+1}| \le \max \left\{ T_{l+1} \left| \frac{\alpha_n}{\alpha(1+q_n)} - 1 \right|, (l+1)\tau_n + T_{l+1} |\alpha_n/\alpha - 1| \right\}. $$

Taking supremum and expectation we get

$$ \mathbf{E} \sup_{l \le M_{T/\tau_n}} |\tau_n K_{l+1} - T_{l+1}| \le \mathbf{E} \max \left\{ o(1) T_{M_{\lfloor T/\tau_n \rfloor}+1}, \tau_n (M_{\lfloor T/\tau_n \rfloor} + 1) + o(1) T_{M_{\lfloor T/\tau_n \rfloor}+1} \right\} = o(1). $$

The last identity holds by observing that, by Lemma 17, for $k := \lfloor T/\tau_n \rfloor$, $M_k \le N_{2\tau_n k}$, hence

$$ \mathbf{E}(M_k + 1) \le \mathbf{E}(N_{2\tau_n k} + 1) = \frac{2\tau_n k}{\alpha} + 1 = O(1), $$

and

$$ \mathbf{E}(T_{M_{\lfloor T/\tau_n \rfloor}+1}) \le \mathbf{E}(T_{N_{2\tau_n k}+1}) \le \frac{2\tau_n k + 1}{\alpha} = O(1). $$

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

**Lemma 19.** We have

$$ \sup_{k \le T/\tau_n} \mathbf{E} \| Y_{T_{M_k}} - Y_{T_{N_{\tau_k}}} \| \le C(\tau_n + |\alpha_n - \alpha|) $$

for some constant $C > 0$ independent of $n$.

*Proof.* Take some $k \le T/\tau_n$. First we show that $|N_{\tau_n k} - M_k|$ is dominated by the number of points $T_j$ that land in a neighborhood of $\tau_n k$ of size $\tau_n O(N_T)$. Let $r$ be a natural number, $l = N_{\tau_n k}$, and observe the following sequence of implications

$$ \begin{align*} N_{\tau_n k} - M_k \ge r &\implies T_l \le \tau_n k < T_{l+1} \text{ and } K_{l-r+1} > k \\ &\implies T_{l-r+1} \le \dots \le T_l \le \tau_n k \le \tau_n K_{l-r+1} \\ &\implies T_{l-r+1}, \dots, T_l \in [\tau_n k - \varepsilon, \tau_n k + \varepsilon], \end{align*} $$
where $\varepsilon = \max\{T|\frac{\alpha_n}{\alpha(1+q_n)} - 1|, \tau_n N_T + T|\alpha_n/\alpha - 1|\}$ and the last implication uses the bound $|T_{l-r+1} - \tau_n K_{l-r+1}| \le \varepsilon$ from Lemma 17. Hence, if $N_{\tau_n k} - M_k \ge r$, then at least $r$ points among $T_1, \dots, T_{N_T}$ fall into the interval $[\tau_n k - \varepsilon, \tau_n k + \varepsilon]$. The same way we can show that $M_k - N_{\tau_n k} \ge r$ implies that at least $r$ points fall into the interval $[\tau_n k - \varepsilon, \tau_n k + \varepsilon]$. Thus we have

$$|N_{\tau_n k} - M_k| \le \sum_{j=1}^{N_T} 1\{T_j \in [\tau_n k - \varepsilon, \tau_n k + \varepsilon]\}$$

Recall that $T_j|N_T = m$ has the same distribution as $j$-th order statistic $u_{(j)}$ of i.i.d. uniform rvs $u_1, \dots, u_m \sim U([0, T])$. Before taking the conditional expectation, let us for convenience bound $\varepsilon$ by $(N_T + 1)\varepsilon'$ where

$$\varepsilon' := \max \left\{ T \middle| \frac{\alpha_n}{\alpha(1+q_n)} - 1 \middle|, \tau_n + T \middle| \frac{\alpha_n}{\alpha} - 1 \middle| \right\}. \quad (24)$$

Hence,

$$\mathbf{E}(|N_{\tau_n k} - M_k| | N_T = m) \le \mathbf{E} \left( \sum_{j=1}^{m} 1\{u_j \in [\tau_n k - (m+1)\varepsilon', \tau_n k + (m+1)\varepsilon']\} \right) = \frac{2\varepsilon'm(m+1)}{T}.$$

Now we can finally derive the upper bound for $\mathbf{E}\|Y_{T_{N_{\tau_n k}}} - Y_{T_{M_k}}\|$. By construction we have for any $\ell, \ell' > 0$ that $\|Y_{T_\ell} - Y_{T_{\ell'}}\| \le |\ell - \ell'| \max \|\alpha \nabla V[\mu_{T_j}](T_j)\| \le \alpha C |\ell - \ell'|$ for some $C > 0$. Finally, using the law of total expectation we get

$$\begin{align*}
\mathbf{E}\|Y_{T_{N_{\tau_n k}}} - Y_{T_{M_k}}\| &\le \sum_{m=1}^{\infty} \alpha \max_{\ell \le m} \|\nabla V[\mu_{T_\ell}](T_\ell)\| \mathbf{E}(|N_{\tau_n k} - M_k| | N_T = m) \mathbf{P}(N_T = m) \\
&\le \alpha C \sum_{m=1}^{\infty} \frac{2\varepsilon'm(m+1)}{T} \cdot e^{-T/\alpha} \frac{(T/\alpha)^m}{m!} \\
&\le C\varepsilon',
\end{align*}$$

where $C$ is redefined on the last line. One can verify by inspection that $C$ is independent of $k$ and $n$. Furthermore, from its definition (24), we see that $\varepsilon'$ is bounded by a constant times $\tau_n + |\alpha_n - \alpha|$, which concludes the proof. $\square$
# A Experimental details

All experiments are coded in JAX [Bradbury et al., 2018]. Our code is provided at https://github.com/PierreMarion23/dropout-phase-diagram.

## A.1 Teacher-student experiment

We train a two-layer ReLU neural networks with *d*-dimensional input and *n* neurons. The neurons are initialized with i.i.d. Gaussians with variance 1/*d* on the inner layer and variance 1 on the outer layer. The input data is standard Gaussian, while the target is generated from a teacher network with the same architecture. The network is trained with (full-batch) gradient descent on the mean squared error. We use the feature learning scaling, meaning that the output of the network is multiplied by 1/*n*, while the learning rate is *n*η₀, where η₀ is a base learning rate. The numerical values of the parameters are given in Table 1.

We compare four optimization algorithms: plain gradient descent (GD), GD with dropout, GD with random geometry (RaM), and GD with propagation noise and random geometry (PN + RaM), see Section 4.4 for details. For the first three of these algorithms, Figure 2 in the main text shows a realization of the weight trajectory for two arbitrary weights. For completeness, we show in Figure 4 the same realization for the four algorithms.

To provide metrics for the comparison between these algorithms, we compute the RMS distance (ℓ₂ distance normalized by $\sqrt{n}$) in parameter space between the dropout version and each of the other three versions. The results are presented in Figure 5a as a function of the number of gradient steps for two different widths, and in Figure 5b as a function of the

<table><caption>Table 1: Parameters for the teacher-student experiment.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Data dimension <i>d</i></td><td>20</td></tr><tr><td>Width <i>n</i></td><td>{200; 500; 1,000; 2,000; 5,000}</td></tr><tr><td>Number of gradient steps</td><td>1,000</td></tr><tr><td>Dropout probability <i>p</i></td><td>0.3</td></tr><tr><td>Training set size</td><td>500</td></tr><tr><td>Width of the teacher</td><td>15</td></tr><tr><td>Base learning rate η<sub>0</sub></td><td>0.5</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 4: Illustration of the pathwise convergence between random metric (RaM) and dropout dynamics (Proposition 3) for the teacher-student experiment. Same realization as in Figure 2, adding GD with propagation noise and random geometry (PN + RaM).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>(a) As a function of the training time, for two widths.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) As a function of the width, at a fixed training time (after 100 steps).</figcaption></figure>

Figure 5: RMS distance in parameter space between the dropout algorithm and the three other tested algorithms (plain GD, GD with random geometry and GD with propagation noise and random geometry) for the teacher-student experiment.

<table><caption>Table 2: Parameters for the MNIST experiment.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Width <i>n</i></td><td>4,000</td></tr><tr><td>Number of gradient steps</td><td>40,000</td></tr><tr><td>Dropout probability <i>p</i></td><td>{0.1, 0.2, 0.3, 0.4}</td></tr><tr><td>Base learning rate <i>η</i><sub>0</sub></td><td>0.25</td></tr><tr><td>Batch size</td><td>64</td></tr></tbody></table>

width for a fixed number of steps. Confidence intervals over 10 repetitions are plotted. We observe that dropout is closer to RaM than to plain GD and PN+RaM, especially for large width and small training time, which is expected from our theory (Proposition 3).

We observe in Figure 6 that, in terms of performance metrics, in this setting, RaM and plain GD both outperform dropout, and perform very similarly to each other. In other words, in this setting the propagation noise from the forward pass is harmful for performance, even more so for the PN+RaM algorithm.

The experiment takes about 30 minutes to run on a consumer laptop.

## A.2 MNIST experiment

We train a two-layer ReLU network on the binary classification of digits 4 and 7 of MNIST. The images are flattened into vectors of dimension 784, and pixel values are mapped into [0, 1]. The loss is the logistic loss. We use the standard data split from MNIST, removing at random 2,000 samples from the training set to serve as validation set. The initialization and scaling of the neural network are as described in Section A.1.

We compare five optimization algorithms: plain mini-batch stochastic gradient descent (SGD, no dropout), SGD with dropout, SGD with random geometry (RaM), SGD with propagation noise (PN), SGD with propagation noise and random geometry (PN + RaM), see Section 4.4 for details. The numerical values of the parameters are given in Table 2.

Confidence intervals over 10 repetitions are plotted. In Figure 3b, we log the loss value over the course of training every 100 steps, and report the minimum logged test loss.

The experiment takes about 3 hours to run on an Nvidia GTX 1080Ti GPU.
<figure><img src="image_11.png" /><figcaption>(a) Train loss for width 200.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Test loss for width 200.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) Train loss for width 5,000.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) Test loss for width 5,000.</figcaption></figure>

Figure 6: Evolution of losses during training for the teacher-student experiment with the four tested optimization algorithms.

## Acknowledgments

P.M. is supported by a Google PhD Fellowship. This work was done in part while P.M. was visiting the Simons Institute for the Theory of Computing.
References

Beatrice Acciaio, Mathias Beiglböck, Friedrich Penkner, Walter Schachermayer, and Johannes Temme. A trajectorial interpretation of Doob's martingale inequalities. *The Annals of Applied Probability*, pages 1494–1505, 2013.

Raman Arora, Peter Bartlett, Poorya Mianjy, and Nathan Srebro. Dropout: Explicit forms and capacity control. In *International Conference on Machine Learning*, pages 351–361. PMLR, 2021.

Martin Bauer and Thilo Meyer-Brandis. Existence and regularity of solutions to multi-dimensional mean-field stochastic differential equations with irregular drift. arXiv preprint arXiv:1912.05932, 2019.

Patrick Billingsley. *Convergence of probability measures*. John Wiley & Sons, 2013.

James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL http://github.com/google/jax.

Pierre Cardaliaguet. Notes on mean field games. Technical report, Université Paris-Dauphine, 2010.

Lénaïc Chizat. The hidden width of deep ResNets: Tight error bounds and phase diagrams. arXiv preprint arXiv:2509.10167, 2025.

Lénaïc Chizat and Francis Bach. On the global convergence of gradient descent for over-parameterized models using optimal transport. *Advances in Neural Information Processing Systems*, 31, 2018.

Lénaïc Chizat and Francis Bach. Implicit bias of gradient descent for wide two-layer neural networks trained with the logistic loss. In *Conference on Learning Theory*, pages 1305–1338. PMLR, 2020.

Lénaïc Chizat, Edouard Oyallon, and Francis Bach. On lazy training in differentiable programming. *Advances in Neural Information Processing Systems*, 32, 2019.

Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. BERT: Pre-training of deep bidirectional transformers for language understanding. In *Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)*, pages 4171–4186. Association for Computational Linguistics, 2019.

Nolan Dey, Bin Claire Zhang, Lorenzo Noci, Mufan Li, Blake Bordelon, Shane Bergsma, Cengiz Pehlevan, Boris Hanin, and Joel Hestness. Don't be lazy: CompleteP enables compute-efficient deep transformers. arXiv preprint arXiv:2505.01618, 2025.

Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby. An image is worth 16x16 words: Transformers for image recognition at scale. In *International Conference on Learning Representations*, 2021.

Mario Geiger, Stefano Spigler, Arthur Jacot, and Matthieu Wyart. Disentangling feature and lazy training in deep neural networks. *Journal of Statistical Mechanics: Theory and Experiment*, 2020(11):113301, 2020.
David P Helmbold and Philip M Long. On the inductive bias of dropout. *The Journal of Machine Learning Research*, 16(1):3403–3454, 2015.

Edward J Hu, Yelong Shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, Lu Wang, and Weizhu Chen. Lora: Low-rank adaptation of large language models. *arXiv preprint arXiv:2106.09685*, 2021.

Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kernel: Convergence and generalization in neural networks. *Advances in Neural Information Processing Systems*, 31, 2018.

Daniel Lacker. Mean field games and interacting particle systems. *Lecture notes*, 2018.

Zhuang Liu, Zhiqiu Xu, Joseph Jin, Zhiqiang Shen, and Trevor Darrell. Dropout reduces underfitting. In *International Conference on Machine Learning*, pages 22233–22248. PMLR, 2023.

Henry P. McKean. Propagation of chaos for a class of non-linear parabolic equations. *Stochastic Differential Equations (Lecture Series in Differential Equations, Session 7, Catholic Univ., 1967)*, pages 41–57, 1967.

Song Mei, Andrea Montanari, and Phan-Minh Nguyen. A mean field view of the landscape of two-layer neural networks. *Proceedings of the National Academy of Sciences*, 115(33): E7665–E7671, 2018.

Song Mei, Theodor Misiakiewicz, and Andrea Montanari. Mean-field theory of two-layers neural networks: dimension-free bounds and kernel limit. In *Conference on Learning Theory*, pages 2388–2464. PMLR, 2019.

Poorya Mianjy and Raman Arora. On convergence and generalization of dropout training. *Advances in Neural Information Processing Systems*, 33:21151–21161, 2020.

Wenlong Mou, Yuchen Zhou, Jun Gao, and Liwei Wang. Dropout training, data-dependent regularization, and generalization bounds. In *International Conference on Machine Learning*, pages 3645–3653. PMLR, 2018.

Long Ouyang, Jeffrey Wu, Xu Jiang, Diogo Almeida, Carroll Wainwright, Pamela Mishkin, Chong Zhang, Sandhini Agarwal, Katarina Slama, Alex Ray, et al. Training language models to follow instructions with human feedback. *Advances in Neural Information Processing Systems*, 35, 2022.

Alec Radford, Karthik Narasimhan, Tim Salimans, and Ilya Sutskever. Improving language understanding by generative pre-training. Technical report, OpenAI, 2018.

Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen. Hierarchical text-conditional image generation with clip latents. *arXiv preprint arXiv:2204.06125*, 1(2): 3, 2022.

Grant Rotskoff and Eric Vanden-Eijnden. Trainability and accuracy of artificial neural networks: An interacting particle system approach. *Communications on Pure and Applied Mathematics*, 75(9):1889–1935, 2022.

Filippo Santambrogio. *Optimal transport for applied mathematicians*, volume 87. Springer, 2015.
Justin Sirignano and Konstantinos Spiliopoulos. Mean field analysis of neural networks: A law of large numbers. *SIAM Journal on Applied Mathematics*, 80(2):725-752, 2020.

Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. Dropout: a simple way to prevent neural networks from overfitting. *The Journal of Machine Learning Research*, 15(1):1929-1958, 2014.

Alain-Sol Sznitman. Topics in propagation of chaos. *Ecole d'été de probabilités de Saint-Flour XIX—1989*, 1464:165–251, 1991.

Ross Taylor, Marcin Kardas, Guillem Cucurull, Thomas Scialom, Anthony Hartshorn, Elvis Saravia, Andrew Poulton, Viktor Kerkez, and Robert Stojnic. Galactica: A large language model for science. *arXiv preprint arXiv:2211.09085*, 2022.

Veeravalli S Varadarajan. On the convergence of sample probability distributions. *Sankhyā: The Indian Journal of Statistics (1933-1960)*, 19(1/2):23–26, 1958.

Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. *Advances in Neural Information Processing Systems*, 30, 2017.

Cédric Villani. *Optimal transport: old and new*, volume 338. Springer, 2008.

Nikhil Vyas, Alexander Atanasov, Blake Bordelon, Depen Morwani, Sabarish Sainathan, and Cengiz Pehlevan. Feature-learning networks are consistent across widths at realistic scales. *Advances in Neural Information Processing Systems*, 36, 2023.

Stefan Wager, Sida Wang, and Percy S Liang. Dropout training as adaptive regularization. *Advances in Neural Information Processing Systems*, 26, 2013.

Colin Wei, Sham Kakade, and Tengyu Ma. The implicit and explicit regularization effects of dropout. In *International Conference on Machine Learning*, pages 10181–10192. PMLR, 2020.

Stephan Wojtowytsch. On the convergence of gradient descent training for two-layer relu-networks in the mean field regime. *arXiv preprint arXiv:2005.13530*, 2020.

Fuzhao Xue, Yao Fu, Wangchunshu Zhou, Zangwei Zheng, and Yang You. To repeat or not to repeat: Insights from scaling LLM under token-crisis. In *Advances in Neural Information Processing Systems*, volume 36, 2023.

Greg Yang and Edward J Hu. Tensor programs IV: Feature learning in infinite-width neural networks. In *International Conference on Machine Learning*, pages 11727–11737. PMLR, 2021.

Zhongwang Zhang and Zhi-Qin John Xu. Implicit regularization of dropout. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 2024.

Zhongwang Zhang, Yuqing Li, Tao Luo, and Zhi-Qin John Xu. Stochastic modified equations and dynamics of dropout algorithm. In *International Conference on Learning Representations*, 2024.