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

HAL
open science

# Particle method and quantization-based schemes for the simulation of the McKean-Vlasov equation

Yating Liu

► To cite this version:

Yating Liu. Particle method and quantization-based schemes for the simulation of the McKean-Vlasov equation. 2022. hal-03915182

HAL Id: hal-03915182

https://hal.science/hal-03915182

Preprint submitted on 29 Dec 2022

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.
# Particle method and quantization-based schemes for the simulation of the McKean-Vlasov equation

Yating Liu *

29th December, 2022

## Abstract

In this paper, we study three numerical schemes for the McKean-Vlasov equation

$$ \begin{cases} dX_t = b(t, X_t, \mu_t) dt + \sigma(t, X_t, \mu_t) dB_t, \\ \forall t \in [0, T], \mu_t \text{ is the probability distribution of } X_t, \end{cases} $$

where $X_0: (\Omega, \mathcal{F}, \mathbb{P}) \to (\mathbb{R}^d, \mathcal{B}(\mathbb{R}^d))$ is a known random variable. Under the assumption on the Lipschitz continuity of the coefficients $b$ and $\sigma$, our first result proves the convergence rate of the particle method with respect to the Wasserstein distance, which extends a previous work [BT97] established in one dimensional setting. In the second part, we present and analyse two quantization-based schemes, including the recursive quantization scheme (deterministic scheme) in the Vlasov setting, and the hybrid particle-quantization scheme (random scheme, inspired by the K-means clustering). Two examples are simulated at the end of this paper: the Burger's equation introduced in [BT97] and the network of FitzHugh-Nagumo neurons (see [BFFT12] and [BFT15]) in dimension 3.

**Keywords:** McKean-Vlasov equation, Mean-field limits, Numerical analysis and simulation, Optimal quantization, Particle method, Quantization-based scheme.

## 1 Introduction

The McKean-Vlasov equation was originally introduced in [McK67] as a stochastic model naturally associated to a class of non-linear PDEs. Nowadays, it refers to the whole family of stochastic differential equations whose coefficients not only depend on the position of the process $X_t$ at time $t$ but also depend on its probability distribution $\mu_t$. This distribution dependent structure of the McKean-Vlasov equation is widely used to model phenomena in Statistical Physics (see e.g. [MA01], [BJR22]), in mathematical biology (see e.g. [BFFT12], [BFT15]), also in social sciences and in quantitative finance, often motivated by the development of the Mean-Field games and the interacting diffusion models (see e.g. [LL18], [CL18] and [CD18]). Besides the original paper [McK67], an excellent reference of the general theory of the McKean-Vlasov equation and the propagation of chaos is [Szn91]. See also the lecture note [Lac18].

In this paper, we consider a filtered probability space $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P})$ satisfying the usual conditions and an $(\mathcal{F}_t)$-standard Brownian motion $(B_t)_{t \ge 0}$ valued in $\mathbb{R}^q$. For a separable Banach space $(E, \|\cdot\|_E)$, we denote the set of probability distributions on $E$ by $\mathcal{P}(E)$ and denote by $\mathcal{P}_p(E)$ the set of probability distributions on $E$ having $p$-th finite moment, $p \ge 1$. The McKean-Vlasov equation writes

$$ \begin{cases} dX_t = b(t, X_t, \mu_t) dt + \sigma(t, X_t, \mu_t) dB_t, \\ \mu_t \text{ is the probability distribution of } X_t, \ t \in [0, T], \end{cases} \quad (1.1) $$

<footnote>*CEREMADE, CNRS, UMR 7534, Université Paris-Dauphine, PSL University, 75016 Paris, France, liu@ceremade.dauphine.fr. This paper is the Chapter 7 of the author's Ph.D. thesis [Liu19].</footnote>
where $X_0 : (\Omega, \mathcal{F}, \mathbb{P}) \to (\mathbb{R}^d, \mathcal{B}(\mathbb{R}^d))$ is a known random variable independent of the Brownian motion $(B_t)_{t \ge 0}$ and $b, \sigma$ are Borel functions defined on $[0, T] \times \mathbb{R}^d \times \mathcal{P}(\mathbb{R}^d)$ having respective values in $\mathbb{R}^d$ and $M_{d,q}(\mathbb{R})$, the space of real matrices with $d$ rows and $q$ columns.

This paper aims to show three implementable numerical methods for the simulation of the McKean-Vlasov equation (1.1), including the particle method and two quantization-based schemes, accompanied by a quantitative analysis of the corresponding simulation error.

In the following discussion, the notation $C_{p_1, \dots, p_n}$ means a positive constant depending on parameters $p_1, \dots, p_n$, $|\cdot|_d$ denotes the Euclidean norm on $\mathbb{R}^d$ (we drop the subscript $d$ when there is no ambiguity) and $\langle \cdot | \cdot \rangle$ denotes the Euclidean inner product. We endow $M_{d,q}(\mathbb{R})$ with an operator norm $\|\cdot\|$ defined by $\|A\| := \sup_{|z|_q \le 1} |Az|_d$. For a random variable $X$, we write $P_X$ or $\mathcal{L}(X)$ for its probability distribution and write $\|X\|_p$ for its $L^p$-norm defined by $\|X\|_p := [\mathbb{E}|X|^p]^{1/p}$, $p \ge 1$. Moreover, let $\mathcal{W}_p$ denote the $L^p$-Wasserstein distance on $\mathcal{P}_p(\mathbb{R}^d)$, $p \ge 1$, defined by

$$ \begin{aligned} \forall \mu, \nu \in \mathcal{P}_p(\mathbb{R}^d), \mathcal{W}_p(\mu, \nu) &:= \left( \inf_{\pi \in \Pi(\mu, \nu)} \int_{\mathbb{R}^d \times \mathbb{R}^d} |x - y|^p \pi(dx, dy) \right)^{\frac{1}{p}} \\ &= \inf \left\{ [\mathbb{E} |X - Y|^p]^{\frac{1}{p}}, X, Y : (\Omega, \mathcal{F}, \mathbb{P}) \to (\mathbb{R}^d, \mathcal{B}(\mathbb{R}^d)) \text{ with } P_X = \mu, P_Y = \nu \right\}, \quad (1.2) \end{aligned} $$

where in the first line of (1.2), $\Pi(\mu, \nu)$ denotes the set of all probability measures on $(\mathbb{R}^d \times \mathbb{R}^d, \mathcal{B}(\mathbb{R}^d)^{\otimes 2})$ with marginals $\mu$ and $\nu$. We denote by $\delta_a$ the Dirac mass at a point $a \in \mathbb{R}^d$.

Throughout this paper, we make the following assumption on the coefficients $b, \sigma$ and on the initial random variable $X_0$. Remark that the following assumption depends on an index $p \in [2, +\infty)$.

**Assumption I.** There exists $p \in [2, +\infty)$ such that the following conditions (1) and (2) hold.

(1) The initial random variable $X_0$ satisfies $\|X_0\|_p < +\infty$.

(2) The coefficient functions $b$ and $\sigma$ are $\gamma$-Hölder continuous in $t$, $\gamma \in (0, 1]$, and Lipschitz continuous in $x$ and in $\mu$ in the following sense: there exists a constant $L > 0$ such that

$$ \begin{aligned} \forall (x, \mu) \in \mathbb{R}^d \times \mathcal{P}_p(\mathbb{R}^d), \forall s, t \in [0, T] \text{ with } s \le t, \\ |b(t, x, \mu) - b(s, x, \mu)| \vee \|\sigma(t, x, \mu) - \sigma(s, x, \mu)\| \le L(1 + |x| + \mathcal{W}_p(\mu, \delta_0))(t - s)^\gamma, \quad (1.3) \end{aligned} $$

and such that

$$ \begin{aligned} \forall t \in [0, T], \forall x, y \in \mathbb{R}^d \text{ and } \forall \mu, \nu \in \mathcal{P}_p(\mathbb{R}^d), \\ |b(t, x, \mu) - b(t, y, \nu)| \vee \|\sigma(t, x, \mu) - \sigma(t, y, \nu)\| \le L(|x - y| + \mathcal{W}_p(\mu, \nu)). \quad (1.4) \end{aligned} $$

Assumption I is a classical assumption for existence and uniqueness of a strong solution $X = (X_t)_{t \in [0, T]}$ to the McKean-Vlasov equation (1.1) (see e.g. [Lac18], [Liu19, Chapter 5]) and the convergence of the Euler scheme of (1.1), defined further in (1.6) (see [LP22, Proposition 2.1]). Besides, in the Vlasov case, that is, there exist $\beta : [0, T] \times \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d$ and $a : [0, T] \times \mathbb{R}^d \times \mathbb{R}^d \to M_{d,q}(\mathbb{R})$ such that

$$ b(t, x, \mu) = \int_{\mathbb{R}^d} \beta(t, x, u)\mu(du) \quad \text{and} \quad \sigma(t, x, \mu) = \int_{\mathbb{R}^d} a(t, x, u)\mu(du), \quad (1.5) $$

a sufficient condition for Assumption I is Assumption $\Pi_V$ (see below), which follows from the Kantorovich-Rubinstein dual representation of the Wasserstein distance $\mathcal{W}_1$:

$$ \mathcal{W}_1(\mu, \nu) = \sup \left\{ \int_{\mathbb{R}^d} f(\xi)\mu(d\xi) - \int_{\mathbb{R}^d} f(\xi)\nu(d\xi) \middle| f \text{ Lipschitz continuous with } [f]_{\text{Lip}} := \sup_{x \ne y} \frac{|f(x) - f(y)|}{|x - y|} \le 1 \right\} $$

and follows from the fact that for every $p \ge 1$, $\mathcal{W}_1(\mu, \nu) \le \mathcal{W}_p(\mu, \nu)$ (see e.g. [Edw11] and [Vil09, Chapter 6]).

**Assumption $\Pi_V$.** The functions $\beta$ and $a$ in (1.5) are $\gamma$-Hölder continuous in $t$, $\gamma \in (0, 1]$, Lipschitz continuous in $(x, u)$ uniformly with respect to $t \in [0, T]$ in the following sense: there exists a constant
$L > 0$ such that

$$ \forall x, y \in \mathbb{R}^d, \forall s, t \in [0, T] \text{ with } s \le t, \\ |\beta(t, x, y) - \beta(s, x, y)| \vee \|a(t, x, y) - a(s, x, y)\| \le L(1 + |x| + |y|)(t - s)^\gamma, $$

and such that

$$ \forall t \in [0, T], \forall x_1, x_2, y_1, y_2 \in \mathbb{R}^d, \\ |\beta(t, x_1, y_1) - \beta(t, x_2, y_2)| \vee \|a(t, x_1, y_1) - a(t, x_2, y_2)\| \le L(|x_1 - x_2| + |y_1 - y_2|). $$

We then turn to the construction of three numerical schemes, which contains two parts: the temporal discretization on $[0, T]$ and the spatial discretization on $\mathbb{R}^d$.

### (1) TEMPORAL DISCRETIZATION BY THEORETICAL EULER SCHEME.

We fix $M \in \mathbb{N}^*$ and define $h = \frac{T}{M}$ as the time step. Let $t_m = t_m^M := m \cdot h$, $0 \le m \le M$ and let $Z_{m+1} := \frac{1}{\sqrt{h}}(B_{t_{m+1}} - B_{t_m})$, $0 \le m \le M-1$. The random variables $Z_m$, $1 \le m \le M$, have standard normal distribution $N(0, I_q)$, where $I_q$ denotes the identity matrix of size $q \times q$. The theoretical Euler scheme of (1.1) is defined as follows,

$$ \bar{X}_{t_{m+1}}^M = \bar{X}_{t_m}^M + h \cdot b(t_m^M, \bar{X}_{t_m}^M, \bar{\mu}_{t_m}^M) + \sqrt{h} \sigma(t_m^M, \bar{X}_{t_m}^M, \bar{\mu}_{t_m}^M) Z_{m+1}, \quad \bar{X}_0^M = X_0, \quad (1.6) $$

where for every $m \in \{0, ..., M\}$, $\bar{\mu}_{t_m}^M$ denotes the probability distribution of $\bar{X}_{t_m}^M$. We also define the continuous expansion of $(\bar{X}_{t_0}^M, ..., \bar{X}_{t_M}^M)$, denoted by $(\bar{X}_t^M)_{t \in [0, T]}$, by the following continuous Euler scheme

$$ \forall m = 0, ..., M-1, \forall t \in (t_m, t_{m+1}], \quad \bar{X}_t := \bar{X}_{t_m}^M + (t - t_m)b(t_m^M, \bar{X}_{t_m}^M, \bar{\mu}_{t_m}^M) + \sigma(t_m^M, \bar{X}_{t_m}^M, \bar{\mu}_{t_m}^M)(B_t - B_{t_m}) \quad (1.7) $$

with the same $\bar{X}_0 = X_0$. When there is no ambiguity, we will omit the superscript $M$ in (1.6) and in (1.7).

Under Assumption I, [LP22, Proposition 2.1] gives the following convergence rate of the theoretical Euler scheme.

**Proposition 1.1.** ([LP22, Proposition 2.1]) Assume that Assumption I holds for an index $p \in [2, +\infty)$. Let $X = (X_t)_{t \in [0,T]}$ be the unique strong solution of (1.1) and let $\bar{X}^M = (\bar{X}_t^M)_{t \in [0,T]}$ denote the process defined by the continuous Euler scheme (1.7). Then there exists a constant $C_{p,d,T,L,\gamma, \|X_0\|_p}$ such that
$$ \| \sup_{t \in [0,T]} |X_t - \bar{X}_t^M| \|_p \le C_{p,d,T,L,\gamma, \|X_0\|_p} h^{\frac{1}{2}\wedge\gamma}. $$

### (2) SPATIAL DISCRETIZATION.

We call (1.6) theoretical Euler scheme since it is not directly implementable, contrary to the Euler scheme for a standard diffusion $dX_t = b(t, X_t)dt + \sigma(t, X_t)dB_t$. The reason is that the scheme (1.6) does not directly indicate how to simulate $\bar{\mu}_{t_m}^M$ in the coefficient functions $b$ and $\sigma$. To do this, we need a further spatial discretization on $\mathbb{R}^d$, a key point of this paper, to construct a discrete approximation of $\bar{\mu}_{t_m}^M$.

#### (2.1) PARTICLE METHOD.
A first way to perform the spatial discretization is called the *particle method*, inspired by the propagation of chaos property of the McKean-Vlasov equation (see e.g. [Szn91], [G88], [Lac18] and [CST22]). In one dimensional setting, the convergence rate of the distribution function and the density function of the particle method has been established in [BT97]. In this paper, we establish the convergence rate with respect to the Wasserstein distance, which also holds for a higher dimensional setting ($d \ge 2$).

We consider the same temporal discretization number $M$ and the same time step $h$ as in (1.6). For the simplicity of notation, we will omit the superscript $M$ in the following discussion. Let $\bar{X}_0^{1,N}$, ..., $\bar{X}_0^{N,N}$ be i.i.d copies of $X_0$ in (1.1) and let $B^n := (B_t^n)_{t \in [0,T]}$, $1 \le n \le N$, be $N$ independent standard Brownian motions valued in $\mathbb{R}^q$, independent of the Brownian motion $(B_t)_{t \in [0,T]}$ in the initial McKean-Vlasov equation (1.1) and of $(X_0, \bar{X}_0^{1,N}, ..., \bar{X}_0^{N,N})$. The main idea of the *particle method* is to construct an
N-particle system $(\bar{X}_{t_m}^{1,N}, ..., \bar{X}_{t_m}^{N,N})_{0\le m\le M}$ by computing for every $1 \le n \le N, 0 \le m \le M - 1$,

$$ egin{cases} 
\bar{X}_{t_{m+1}}^{n,N} = \bar{X}_{t_m}^{n,N} + hb(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) + \sqrt{h}\sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N)Z_{m+1}^n, \\ 
\bar{\mu}_{t_m}^N := \frac{1}{N}\sum_{n=1}^N \delta_{\bar{X}_{t_m}^{n,N}}, \quad Z_{m+1}^n = \frac{1}{\sqrt{h}}(B_{t_{m+1}}^n - B_{t_m}^n) \stackrel{\text{i.i.d}}{\sim} \mathcal{N}(0, \mathbf{I}_q), 
\end{cases} 
\quad (1.8) $$

and then to use $\bar{\mu}_{t_m}^N$ defined in (1.8) as an estimator of $\bar{\mu}_{t_m}$ in (1.6) at each time step $t_m, 1 \le m \le M$. Remark that $\bar{\mu}_{t_m}^N$ in (1.8) is a random probability distribution, that is, it depends on $\omega \in \Omega$.

Our first main result is the convergence rate of $W_p(\bar{\mu}_{t_m}^N, \bar{\mu}_{t_m})$ in $L^p$ for every $0 \le m \le M$. Let $(C([0, T], \mathbb{R}^d), \|\cdot\|_{\text{sup}})$ denote the space of $\mathbb{R}^d$-valued continuous applications defined on $[0, T]$, equipped with the uniform norm $\|\alpha\|_{\text{sup}} := \sup_{t \in [0, T]} |\alpha_t|$. Let $W_p$ denote the Wasserstein distance on

$$ \mathcal{P}_p(C([0, T], \mathbb{R}^d)) := \left\{ \mu \text{ probability distribution on } C([0, T], \mathbb{R}^d) \text{ s.t. } \int_{C([0, T], \mathbb{R}^d)} \|\alpha\|_{\text{sup}}^p \mu(d\alpha) < +\infty \right\}, \quad (1.9) $$

which is defined for every $\mu, \nu \in \mathcal{P}_p(C([0, T], \mathbb{R}^d))$ by

$$ \begin{aligned} W_p(\mu, \nu) &:= \left( \inf_{\pi \in \Pi(\mu, \nu)} \int_{C([0,T], \mathbb{R}^d) \times C([0,T], \mathbb{R}^d)} \|\alpha - \tilde{\alpha}\|_{\text{sup}} \pi(d\alpha, d\tilde{\alpha}) \right)^{\frac{1}{p}} \\ &= \inf \left\{ \left[ \mathbb{E} \|X - Y\|_{\text{sup}}^p \right]^{\frac{1}{p}}, X, Y : (\Omega, \mathcal{F}, \mathbb{P}) \to (C([0,T], \mathbb{R}^d), \|\cdot\|_{\text{sup}}) \text{ with } P_X = \mu, P_Y = \nu \right\}, \quad (1.10) \end{aligned} $$

where in the first line of (1.10), $\bar{\Pi}(\mu, \nu)$ denotes the set of all probability measures on $C([0, T], \mathbb{R}^d) \times C([0, T], \mathbb{R}^d)$ (equipped with the Borel $\sigma$-algebra generated by $\|\cdot\|_{\text{sup}}$) having marginals $\mu$ and $\nu$.

**Theorem 1.2.** Assume that Assumption I is in force for an index $p \in [2, +\infty)$. Set $M \in \mathbb{N}^*$ and $h = \frac{T}{M}$ for the temporal discretization. Let $(\bar{\mu}_{t_m})_{1 \le m \le M}$, $(\bar{\mu}_{t_m}^N)_{1 \le m \le M}$ be probability distributions respectively defined by the theoretical Euler scheme (1.6) and the particle method (1.8).

(i) We have

$$ \left\| \sup_{0 \le m \le M} W_p(\bar{\mu}_{t_m}^N, \bar{\mu}_{t_m}) \right\|_p \le C_{d,p,L,T} \|W_p(\bar{\mu}, \nu^N)\|_p, $$

where $\bar{\mu}$ is the probability distribution of the process $\bar{X} = (\bar{X}_t)_{t \in [0, T]}$ defined by (1.7) and $\nu^N := \frac{1}{N} \sum_{i=1}^N \delta_{Y^i}$ is the empirical measure on the i.i.d. copies $(Y^1, ..., Y^N)$ of $\bar{X}$.

(ii) If, in addition, $\|X_0\|_{p+\varepsilon} < +\infty$ for some $\varepsilon > 0$, we have the following inequality

$$ \left\| \sup_{0 \le m \le M} W_p(\bar{\mu}_{t_m}^N, \bar{\mu}_m) \right\|_p \le \tilde{C} \times \begin{cases} N^{-\frac{1}{2p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p > d/2 \text{ and } \varepsilon \ne p, \\ N^{-\frac{1}{2p}} \left[ \log(1+N) \right]^{\frac{1}{p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p = d/2 \text{ and } \varepsilon \ne p, \\ N^{-\frac{1}{d}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p \in (0, d/2) \text{ and } p + \varepsilon \ne \frac{d}{(d-p)}, \end{cases} $$

where $\tilde{C}$ is a constant depending on $p, \varepsilon, d, b, \sigma, L, T$ and $\|X_0\|_{p+\varepsilon}$.

Theorem 1.2 can be considered as a strong error of the particle method. The weak error, that is, the upper-bound of $|\Phi(\bar{\mu}_{t_m}) - \mathbb{E}\Phi(\bar{\mu}_{t_m}^N)|$ for a function $\Phi: \mathcal{P}_p(\mathbb{R}^d) \to \mathbb{R}^d$ with appropriate differentiability, can be established in the future by applying similar techniques as [CST22]. Moreover, we refer to [AKH02] and [HL22] for the studies of the density simulation of $\bar{\mu}_{t_m}$ based on the particle method (1.8).

(2.2) **QUANTIZATION-BASED SCHEME.** A second way to implement the spatial discretization is to use the (optimal) vector quantization, also known as K-means clustering in unsupervised learning.

Consider a probability distribution $\mu$ on $\mathbb{R}^d$. The main idea of the vector quantization is to use the projection of $\mu$ on a fixed quantizer (called also quantization grid in the literature) $x = (x_1, ..., x_K) \in (\mathbb{R}^d)^K$ as an approximation of $\mu$. More specifically, for a fixed quantization level $K \in \mathbb{N}^*$, for a quantizer
$x = (x_1, ..., x_K)$ satisfying $x_i \neq x_j$ if $i \neq j$, this projection, denoted by $\hat{\mu}^x$, is defined by

$$ \hat{\mu}^x := \sum_{k=1}^{K} \mu(V_k(x)) \cdot \delta_{x_k} = \mu \circ \mathrm{Proj}_x^{-1} \qquad (1.11) $$

where $(V_k(x))_{1 \leq k \leq K}$ denotes a Voronoï partition generated by $x = (x_1, ..., x_K)$, which is a Borel partition on $\mathbb{R}^d$ satisfying

$$ V_k(x) \subset \{y \in \mathbb{R}^d \mid |y - x_k| = \min_{1 \leq j \leq K} |y - x_j|\}, \quad 1 \leq k \leq K, \qquad (1.12) $$

(see Figure 1) and where the projection function $\mathrm{Proj}_x$ is defined by

$$ \xi \in \mathbb{R}^d \longmapsto \mathrm{Proj}_x(\xi) := \sum_{k=1}^{K} x_k \mathbb{1}_{V_k(x)}(\xi) \in \{x_1, ..., x_K\} \qquad (1.13) $$

based on the chosen Voronoï partition $(V_k(x))_{1 \leq k \leq K}$. (1)

The approximation $\hat{\mu}^x$ defined in (1.11) is deterministic. Moreover, if $\mu \in \mathcal{P}_p(\mathbb{R}^d)$, $p \geq 1$, there exists (at least) an optimal quantizer $x^* = (x_1^*, ..., x_K^*) \in (\mathbb{R}^d)^K$ in the sense of

$$ \mathcal{W}_p(\mu, \hat{\mu}^{x^*}) = \inf \{ \mathcal{W}_p(\mu, \nu) \mid \nu \in \mathcal{P}(\mathbb{R}^d) \text{ s.t. the support of } \nu \text{ has at most } K \text{ points.} \}. \qquad (1.14) $$

The equality (1.14) shows that $\hat{\mu}^{x^*}$ is the closest probability measure to $\mu$ with respect to the Wasserstein distance $\mathcal{W}_p$, among all probability distributions having a support of at most $K$ points ([GL00, Lemma 3.4]). Moreover, if $\mu \in \mathcal{P}_{p+\varepsilon}(\mathbb{R}^d)$ for some $\varepsilon > 0$, the optimal error $\mathcal{W}_p(\mu, \hat{\mu}^{x^*})$ has the following upper bound with convergence rate $K^{-1/d}$ (Zador's Theorem, see e.g. [LP08] and [Pag18, Theorem 5.2])

$$ \mathcal{W}_p(\mu, \hat{\mu}^{x^*}) \leq C_{d,p,\varepsilon} \mathcal{M}_{p+\varepsilon}(\mu) K^{-1/d} \quad \text{with} \quad \mathcal{M}_{p+\varepsilon}(\mu) := \inf_{a \in \mathbb{R}^d} \left[ \int_{\mathbb{R}^d} |\xi - a|^{p+\varepsilon} \mu(d\xi) \right]^{\frac{1}{p+\varepsilon}}. $$

In the quadratic setting ($p=2$), the optimal quantizer $x^*$ can be found by several numerical methods such as Lloyd's algorithm (see further Algorithm 1 and [Llo82], [Kie82], [PY16]) or the CLVQ algorithm (see e.g. [Pag15, Section 3.2]). Figure 2 is an illustration of the quadratic optimal quantization at level $K=60$ for the standard normal distribution $\mu = \mathcal{N}(0, \mathbf{I}_2)$, where the blue points represent an optimal quantizer and the colour in red represents the weight $\mu(V_k(x))$ of each Voronoï cell (the darker the heavier).

<figure><img src="image_2.png" /><figcaption>Figure 1: A Voronoï partition.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: The quadratic optimal quantization for the standard normal distribution $\mathcal{N}(0, \mathbf{I}_2)$.</figcaption></figure>

In the framework of unsupervised learning with unlabeled dataset $\{\xi_1, ..., \xi_n\} \subset \mathbb{R}^d$, the K-means

<footnote n="1">¹The Voronoï partition $(V_k(x))_{1 \leq k \leq K}$ generated by a fixed quantizer $x = (x_1, ..., x_K)$ is not unique since we can place points on the hyperplane $H_{i,j} := \{\xi \in \mathbb{R}^d \mid |x_i - \xi| = |x_j - \xi|\}$ in either the Voronoï cell $V_i(x)$ or $V_j(x)$. However, the choice of the Voronoï partition does not impact on the quantization error of the quantizer $x$ (see further discussion in (3.1)-(3.2)).</footnote>
clustering is to compute the (quadratic) optimal quantizer of the empirical measure $\nu = \frac{1}{n} \sum_{i=1}^{n} \delta_{\xi_i}$ on the dataset.

In this paper, we proposed the following two quantization-based schemes for the simulation of the McKean-Vlasov equation (2):

(a) Recursive quantization scheme in the Vlasov setting (1.5),

(b) Hybrid particle-quantization scheme.

Different from the particle method, the recursive quantization scheme is deterministic, that is, the simulation result does not depend on $\omega \in \Omega$.

(a) **RECURSIVE QUANTIZATION SCHEME.** The most natural idea of applying optimal quantization to the simulation of the McKean-Vlasov equation is to replace $\tilde{X}_{t_m}$ and $\tilde{\mu}_{t_m}$ in (1.6) by their respective quantized projection. More precisely, we consider a quantizer sequence for each time step $t_m$

$$x^{(m)} = (x_1^{(m)}, \dots, x_K^{(m)}) \in (\mathbb{R}^d)^K, \quad 0 \le m \le M \qquad (1.15)$$

and define the quantized random variable $\tilde{X}_{t_m}$ and its probability distribution $\hat{\mu}_{t_m}$ by

$$\left\{ \begin{array}{l} \tilde{X}_0 = X_0, \quad \tilde{X}_0 = \operatorname{Proj}_{x^{(0)}}(\tilde{X}_0), \quad \hat{\mu}_0 = P_{\tilde{X}_0}, \\ \tilde{X}_{t_{m+1}} = \tilde{X}_{t_m} + h \cdot b(t_m, \tilde{X}_{t_m}, \hat{\mu}_{t_m}) + \sqrt{h} \sigma(t_m, \tilde{X}_{t_m}, \hat{\mu}_{t_m}) Z_{m+1}, \\ \tilde{X}_{t_{m+1}} = \operatorname{Proj}_{x^{(m+1)}}(\tilde{X}_{t_{m+1}}), \quad \hat{\mu}_{t_{m+1}} = P_{\tilde{X}_{t_{m+1}}}, \quad m = 0, \dots, M-1. \end{array} \right. \qquad (1.16)$$

In (1.16), at each time step $t_m$, we add a quantization procedure for $\tilde{X}_{t_m}, \tilde{\mu}_{t_m} = \mathcal{L}(\tilde{X}_{t_m})$ and inject their quantization projection $\tilde{X}_{t_m} = \operatorname{Proj}_{x^{(m)}}(\tilde{X}_{t_m})$ as well as $\hat{\mu}_{t_m} = \mathcal{L}(\tilde{X}_{t_m}) = \tilde{\mu}_{t_m} \circ \operatorname{Proj}_{x^{(m)}}^{-1}$ into the time step $t_{m+1}$. This scheme (1.16) is not directly implementable because of $\hat{\mu}_{t_m}$ for the same reason as the theoretical Euler scheme (1.6), so we call (1.16) *theoretical quantization-based scheme*. However, in the Vlasov setting (1.5), we can circumvent this issue by using the recursive quantization method, which is firstly introduced in [PS15] for the diffusion equation $dX_t = b(t, X_t)dt + \sigma(t, X_t)dB_t$.

The main idea of the recursive quantization method is to construct the Markovian transitions of $(\tilde{X}_{t_m}, \hat{\mu}_{t_m})$ based on (1.16). To be more specific, for each time step $t_m$, by considering Voronoï partitions $(V_k(x^{(m)}))_{1 \le k \le K}$ generated by $x^{(m)}, 0 \le m \le M$, and by applying (1.11), we can rewrite $\hat{\mu}_{t_m}$ in (1.16) as follows

$$\hat{\mu}_{t_m} = \sum_{k=1}^{K} p_k^{(m)} \delta_{x_k^{(m)}} \qquad (1.17)$$

where $p_k^{(m)} = \tilde{\mu}_{t_m}(V_k(x^{(m)})) = \mathbb{P}(\tilde{X}_{t_m} \in V_k(x^{(m)})) = \mathbb{P}(\tilde{X}_{t_m} = x_k^{(m)}), 1 \le k \le K$. The expression (1.17) shows that for every $m \in \{0, \dots, M\}$, the quantized distribution $\hat{\mu}_{t_m}$ is determined by the quantizer $x^{(m)} = (x_1^{(m)}, \dots, x_K^{(m)})$ and its weight vector $p^{(m)} = (p_1^{(m)}, \dots, p_K^{(m)})$. In the Vlasov setting (1.5), the transition step of (1.16), with respect to $\hat{\mu}_{t_m}$ in (1.17), writes

$$\tilde{X}_{t_{m+1}} = \hat{X}_{t_m} + h \sum_{k=1}^{K} p_k^{(m)} \beta(t_m, \hat{X}_{t_m}, x_k^{(m)}) + \sqrt{h} \sum_{k=1}^{K} p_k^{(m)} a(t_m, \hat{X}_{t_m}, x_k^{(m)}) Z_{m+1}.$$

Consequently, given $\hat{X}_{t_m}$ and $p^{(m)}$, the transition probability is

$$\begin{align} \mathbb{P}(\hat{X}_{t_{m+1}} = x_j^{(m+1)} | \hat{X}_{t_m} = x_i^{(m)}, p^{(m)}) \tag{1.18} \\ = \mathbb{P}\left[ \left(x_i^{(m)} + h \sum_{k=1}^{K} p_k^{(m)} \beta(t_m, x_i^{(m)}, x_k^{(m)}) + \sqrt{h} \sum_{k=1}^{K} p_k^{(m)} a(t_m, x_i^{(m)}, x_k^{(m)}) Z_{m+1}\right) \in V_j(x^{(m+1)}) \right] \nonumber \end{align}$$

<footnote n="2">²This idea of applying optimal quantization to simulate the McKean-Vlasov equation was firstly introduced in [GPPP05][Section 4] in a different framework. Besides, another quantization-based scheme is proposed in [Liu19, Section 7.4], called *doubly quantized scheme*, in which we implement the quantized Gaussian random variables instead of $(Z_m)_{1 \le m \le M}$ in (1.6).</footnote>
so that given $p^{(m)}$, we can compute $p_j^{(m+1)} = \mathbb{P}(\hat{X}_{t_{m+1}} = x_j^{(m+1)}), 1 \le j \le K$, by

$$p_j^{(m+1)} = \mathbb{P}(\hat{X}_{t_{m+1}} = x_j^{(m+1)} | p^{(m)}) = \sum_{i=1}^{K} \mathbb{P}(\hat{X}_{t_{m+1}} = x_j^{(m+1)} | \hat{X}_{t_m} = x_i^{(m)}, p^{(m)}) \cdot p_i^{(m)}. \quad (1.19)$$

A detailed proof of the above equalities is provided in Section 3.2, where we also explain how to combine this recursive quantization method with Lloyd's algorithm to optimally compute the quantizer $x^{(m)}$ at each time step. The following theorem, whose proof is in Section 3.1, shows the error analysis of the quantization based scheme (1.16).

**Theorem 1.3 (Error analysis of the quantization-based scheme).** Assume that Assumption I is satisfied with $p = 2$. Set $M \in \mathbb{N}^*$ and $h = \frac{T}{M}$ for the temporal discretization. Let $(\hat{X}_{t_m})_{0 \le m \le M}$ be random variables defined by the Euler scheme (1.6). Consider a fixed K-level quantizer sequence $x^{(m)} = (x_1^{(m)}, \dots, x_K^{(m)})$, $0 \le m \le M$ and define $(\hat{X}_{t_m})_{0 \le m \le M}$, $(\hat{X}_{t_m})_{0 \le m \le M}$ by the quantization-based scheme (1.16).

(1) For every $m \in \{1, \dots, M\}$, we have

$$\| \bar{X}_{t_m} - \hat{X}_{t_m} \|_2 \le \sum_{j=0}^{m} [1 + 2Lh(1 + Lh + Lq)]^j \cdot \Xi_{m-j}, \quad (1.20)$$

where for every $j \in \{0, \dots, m\}$, $\Xi_j$ denotes the quadratic quantization error of the j-th step, that is, $\Xi_j = \|\tilde{X}_{t_j} - \hat{X}_{t_j}\|_2$.

(2) If moreover, there exists $\varepsilon > 0$ such that $\|X_0\|_{2+\varepsilon} < +\infty$ and if for every $0 \le m \le M$, $x^{(m)}$ is a quadratic optimal quantizer of $\tilde{X}_{t_m}$, we have

$$\| \bar{X}_{t_m} - \hat{X}_{t_m} \|_2 \le C_{b, \sigma, L, T, d, \varepsilon, q, \|X_0\|_{2+\varepsilon}} \cdot K^{-1/d} \left( \sum_{j=0}^{m} [1 + 2Lh(1 + Lh + Lq)]^j \right). \quad (1.21)$$

(b) HYBRID PARTICLE-QUANTIZATION SCHEME. Another way to simulate the McKean-Vlasov equation is to apply the optimal quantization method to the particle system (1.8) and subsequently devise the following *hybrid particle-quantization scheme*. Consider the same initial random variables $\bar{X}_0^{1,N}$, ..., $\bar{X}_0^{N,N}$ and the same Brownian motions $B^n$, $1 \le n \le N$ as (1.8) and define

$$\left\{ \begin{array}{l} \tilde{X}_{t_{m+1}}^{n,N} = \tilde{X}_{t_m}^{n,N} + h \cdot b(t_m, \tilde{X}_{t_m}^{n,N}, \hat{\mu}_{t_m}^K) + \sqrt{h} \sigma(t_m, \tilde{X}_{t_m}^{n,N}, \hat{\mu}_{t_m}^K) Z_{m+1}^n, \quad 1 \le n \le N, \\ \hat{\mu}_{t_m}^K = \left( \frac{1}{N} \sum_{n=1}^N \delta_{\tilde{X}_{t_m}^{n,N}} \right) \circ \mathrm{Proj}_{x^{(m)}}^{-1} = \sum_{k=1}^K \left[ \delta_{x_k^{(m)}} \cdot \sum_{n=1}^N \mathbb{1}_{V_k(x^{(m)})(\tilde{X}_{t_m}^{n,N})} \right], \quad (1.22) \\ \tilde{X}_{t_0}^{n,N} = \bar{X}_0^{n,N} \quad \text{i.i.d} \quad \tilde{X}_0, \quad Z_m^n = \frac{1}{\sqrt{h}} (B_{t_{m+1}}^n - B_{t_m}^n) \quad \text{i.i.d} \quad \mathcal{N}(0, \mathbf{I}_q). \end{array} \right. $$

This scheme adds a quantization step of the empirical measure $\frac{1}{N} \sum_{i=1}^N \delta_{\bar{X}_{t_m}^{n,N}}$ at time $t_m$ and injects the quantized measure $\hat{\mu}_{t_m}^K$ as an input to the simulation for $\tilde{X}_{t_{m+1}}^{n,N}$. This scheme is inspired by the consistency of the optimal quantization established in [LP20]. Namely, we can envisage that an optimal quantizer of the particle system at time $t_m$ is quasi-optimal for the measure $\bar{\mu}_{t_m}$ of (1.6), as a consequence of the convergence of the particle method obtained in Theorem 1.2.

The error analysis of the hybrid scheme (1.22) is established in the following proposition.

**Proposition 1.4.** Set the same temporal discretization as Theorem 1.2 and Theorem 1.3. Assume that Assumption I holds true with $p = 2$. For any $m \in \{1, \dots, M\}$, let $\bar{\mu}_{t_m}^N$ be the empirical measure on the particles $(\bar{X}_{t_m}^{0,N}, \dots, \bar{X}_{t_m}^{N,N})$, defined by the particle system (1.8) and let $\hat{\mu}_{t_m}^K$ be the quantized measure defined in (1.22). Then,

$$\forall 1 \le m \le M, \quad \mathbb{E} \left[ \mathcal{W}_2(\hat{\mu}_{t_m}^K, \bar{\mu}_{t_m}^N) \right] \le C_2 \sum_{j=0}^{m-1} C_1^j \sqrt{\mathbb{E} \hat{\Xi}_{m-1-j}^2} + \mathbb{E} \hat{\Xi}_m. \quad (1.23)$$
where $\hat{\Xi}_m = W_2(\hat{\mu}_{t_m}^K, \frac{1}{N} \sum_{i=1}^N \delta_{\bar{X}_{t_m}^{n,N}})$ denotes the quadratic quantization error at time $t_m$ and $C_1, C_2$ are positive constants depending on $L, q$ and $T$.

**Remark 1.5.** We would like to highlight that the upper-bounds in (1.20), (1.21) and (1.23) are not numerically optimal since the two sums $\sum_{j=0}^M [1+2Lh(1+Lh+Lq)]^j$ and $\sum_{j=0}^{M-1} C_1^j \sqrt{E\hat{\Xi}_{m-1-j}^2}$ at the right hand side converge to infinity as the temporal discretization number $M \to +\infty$. These are common error bounds when dealing with Markovian quantization based scheme (see e.g. [PS15]). However, we do not observe on numerical experiments (see further Section 4) a large simulation error, even when we consider a large values of $M$.

This paper is organised as follows.

In Section 2, we prove Theorem 1.2, the convergence rate of the particle method. The proof shares the same idea of the well-known propagation of chaos (see e.g. [Lac18, Theorem 3.3]). We compare the particle system (1.8) with another particle system without interaction, that is, a system composed by $N$ i.i.d. Itô processes $(Y_t^1, ..., Y_t^N)_{t \in [0,T]}$ simulated by the continuous Euler scheme (1.7). Then the error of the particle method can be bounded by, up to a constant multiplier, the Wasserstein distance between the probability measure of $\bar{X}$ defined by (1.7) and its empirical measure defined on its i.i.d. copies $(Y^1, ..., Y^N)$.

Section 3 is devoted to the quantization-based schemes (1.16), (1.18) and (1.22) and their error analyses. We first give a review of the optimal quantization method in Subsection 3.1 and prove Theorem 1.3, the $L^2$-error analysis of the theoretical Euler scheme (1.16). Then in Subsection 3.2, we show how to obtain the transition probability (1.18) for the recursive quantization scheme in the Vlasov setting, as well as the integration of Lloyd's algorithm. Finally, Subsection 3.3 contains the proof of Proposition 1.4, the $L^2$-error analysis of the hybrid particle-quantization scheme (1.22).

In Section 4, we illustrate two simulation examples by using the above numerical methods. The first simulation is for a one-dimensional Burgers equation, introduced in [Szn91] and [BT97]. This equation has a closed form solution so we can compute and compare the accuracy of different methods. The second example is a 3-dimensional model, the network of FitzHugh-Nagumo neurons, firstly introduced and simulated in [BFFT12] (see also some corrections in [BFT15]) and also simulated in [dRES22].

Finally, for the reader's convenience, all numerical methods presented in this paper and their connection are displayed in Appendix B-Figure 18, in which we also briefly mention their convergence rates.

**Remark 1.6.** (Comments on the particle method and on the hybrid particle-quantization scheme) In many research areas and applications related to the McKean-Vlasov equation such as mean field games (see e.g. [CD18]) and opinion dynamics (see e.g. [HK02]), the main research target is an $N$-particle system $(X_t^{1,N}, ..., X_t^{N,N})_{t \in [0,T]}$ satisfying $(X_0^{1,N}, ..., X_0^{N,N}) \stackrel{\text{i.i.d.}}{\sim} X_0$ and

$$ \begin{cases} dX_t^{n,N} = b(t, X_t^{n,N}, \mu_t^N)dt + \sigma(t, X_t^{n,N}, \mu_t^N)dB_t^n, & 1 \le n \le N, \\ \forall t \in [0, T], & \mu_t^N = \frac{1}{N} \sum_{n=1}^N \delta_{X_t^{n,N}}. \end{cases} \quad (1.24) $$

The McKean-Vlasov equation (1.1) is considered as a limit equation of such particle system (1.24) when $N \to +\infty$, often appeared in the propagation of chaos theory in the literature. The particle method (1.8) in this paper can be considered as a temporal discretization of (1.24). Furthermore, from a point of view of unsupervised learning, the hybrid particle-quantization scheme (1.22) is essentially adding a $K$-means clustering step of the particle method, which can be computed by e.g. Python package `sklearn.cluster.KMeans` in practice. This adding step could be itself heuristic for the modelling since $K$-means clustering is usually used for feature extraction in unsupervised learning. For example, we can think the $K$-means cluster centers model some representative particles of the system. In Section 4.2, we show by numerical experiments that the hybrid particle-quantization scheme reduces the data volume while keeping a kind of stability on the mean and variance of the test function value.
## 2 Particle method and its convergence rate

In this section, we study the convergence of the particle method and prove Theorem 1.2 in three steps.

* Step 1 (Subsection 2.1): We state the space of stochastic processes where live the strong solution $X = (X_t)_{t \in [0,T]}$ of the McKean-Vlasov equation (1.1) and the process $\bar{X} = (\bar{X}_t)_{t \in [0,T]}$ defined by the continuous Euler scheme (1.7). We also define the space of probability distributions and the space of marginal distributions of these two processes $X = (X_t)_{t \in [0,T]}$, $\bar{X} = (\bar{X}_t)_{t \in [0,T]}$.

* Step 2 (Subsection 2.2): The system $(\bar{X}_{t_m}^{1,N}, \dots, \bar{X}_{t_m}^{N,N})_{0 \le m \le M}$ defined by (1.8) can be considered as a particle system with interaction through $\bar{\mu}_{t_m}^N = \frac{1}{N} \sum_{n=1}^N \delta_{\bar{X}_{t_m}^{n,N}}$ at each time $t_m$, $0 \le m \le M$. We define another particle system $(Y_t^1, \dots, Y_t^N)_{t \in [0,T]}$ without interaction, that is, $Y^n$, $1 \le n \le N$ are i.i.d. Itô processes with the same distribution as $\bar{X} = (\bar{X}_t)_{t \in [0,T]}$ defined by (1.7).

* Step 3 (Subsection 2.3): We compare these two particle systems and prove Theorem 1.2.

### 2.1 Spaces of stochastic processes and their probability distributions

Recall that $(C([0, T], \mathbb{R}^d), \|\cdot\|_{\sup})$ denotes the space of $\mathbb{R}^d$-valued continuous applications defined on $[0, T]$ and valued in $\mathbb{R}^d$, equipped with the uniform norm $\|\alpha\|_{\sup} := \sup_{t \in [0,T]} |\alpha_t|$. Let

$$L_{C([0,T], \mathbb{R}^d)}^p (\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P}) := \left\{ \begin{array}{l} Y = (Y_t)_{t \in [0,T]} : (\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P}) \to C([0,T], \mathbb{R}^d) \text{ such that} \\ Y \text{ is } (\mathcal{F}_t)\text{-adapted and } \|Y\|_p := [\mathbb{E} \|Y\|_{\sup}^p]^{1/p} < +\infty \end{array} \right\}. \quad (2.1)$$

For any random variable $Y \in L_{C([0,T], \mathbb{R}^d)}^p (\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P})$, its probability distribution $P_Y$ naturally lies in $\mathcal{P}_p(C([0,T], \mathbb{R}^d))$ defined by (1.9). Moreover, the space $\mathcal{P}_p(C([0,T], \mathbb{R}^d))$ equipped with the Wasserstein distance $W_p$ defined by (1.10) is complete and separable since $(C([0,T], \mathbb{R}^d), \|\cdot\|_{\sup})$ is a Polish space (see [Bol08]).

**Lemma 2.1.** ([LP22, Proposition 2.1]) Assume that Assumption I holds for an index $p \in [2, +\infty)$. Then the unique solution $X = (X_t)_{t \in [0,T]}$ of (1.1) and the process $\bar{X}^M = (\bar{X}_t^M)_{t \in [0,T]}$ defined by the continuous Euler scheme (1.7) satisfy the following inequality

$$\forall M \ge 1, \quad \left\| \sup_{t \in [0,T]} |X_t| \right\|_p \vee \left\| \sup_{t \in [0,T]} |\bar{X}_t^M| \right\|_p \le C_{p,d,b,\sigma,T,L} (1 + \|X_0\|_p), \quad (2.2)$$

where $C_{p,d,b,\sigma,T,L}$ is a constant depending on $p, d, b, \sigma, T, L$. Hence, $X, \bar{X} \in L_{C([0,T], \mathbb{R}^d)}^p (\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P})$ and their respective probability distributions $\mu = \mathbb{P} \circ X^{-1}$ and $\bar{\mu} = \mathbb{P} \circ \bar{X}^{-1}$ lie in $\mathcal{P}_p(C([0,T], \mathbb{R}^d))$.

Let us consider now

$$C([0, T], \mathcal{P}_p(\mathbb{R}^d)) := \left\{ (\mu_t)_{t \in [0, T]} \text{ s.t. } t \mapsto \mu_t \text{ is a continuous application from } [0, T] \text{ to } (\mathcal{P}_p(\mathbb{R}^d), \mathcal{W}_p) \right\}$$

equipped with the distance

$$d_C((\mu_t)_{t \in [0,T]}, (\nu_t)_{t \in [0,T]}) := \sup_{t \in [0,T]} \mathcal{W}_p(\mu_t, \nu_t). \quad (2.3)$$

As $(\mathcal{P}_p(\mathbb{R}^d), \mathcal{W}_p)$ is a complete space (see [Bol08]), the space $C([0, T], \mathcal{P}_p(\mathbb{R}^d))$ equipped with the uniform distance $d_C$ is also complete. For any $t \in [0, T]$, we define $\pi_t : C([0, T], \mathbb{R}^d) \to \mathbb{R}^d$ by $\alpha \mapsto \pi_t(\alpha) = \alpha_t$ and define an application $\iota : \mathcal{P}_p(C([0, T], \mathbb{R}^d)) \to C([0, T], \mathcal{P}_p(\mathbb{R}^d))$ by

$$\mu \mapsto \iota(\mu) = (\mu \circ \pi_t^{-1})_{t \in [0,T]} = (\mu_t)_{t \in [0,T]} \quad (2.4)$$

The well-posedness of the application $\iota$ has been proved in [Liu19, Lemma 5.1.2], which also implies that the marginal distributions of processes $X$ and $\bar{X}$ lie in the space $C([0, T], \mathcal{P}_p(\mathbb{R}^d))$.
**Lemma 2.2.** Let $(\mu_t)_{t \in [0,T]}$ and $(\bar{\mu}_t)_{t \in [0,T]}$ be respective marginal distributions of the unique solution $X$ of (1.1) and the process $\bar{X}$ defined by the continuous Euler scheme (1.7). Then $(\mu_t)_{t \in [0,T]}$, $(\bar{\mu}_t)_{t \in [0,T]} \in C([0,T], \mathcal{P}_p(\mathbb{R}^d))$.

In [Lac18], the author defines an application $W_{p,t}$, called “truncated Wasserstein distance”, on $\mathcal{P}_p(C([0,T], \mathbb{R}^d)) \times \mathcal{P}_p(C([0,T], \mathbb{R}^d))$ by

$$ \forall \mu, \nu \in \mathcal{P}_p(C([0,T], \mathbb{R}^d)), W_{p,t}(\mu, \nu) := \left[ \inf_{\pi \in \Pi(\mu, \nu)} \int_{C([0,T], \mathbb{R}^d) \times C([0,T], \mathbb{R}^d)} \sup_{s \in [0,t]} |x_s - y_s|^p \pi(dx, dy) \right]^{\frac{1}{p}}, \quad (2.5) $$

where $\bar{\Pi}(\mu, \nu)$ is the same set as in (1.10). The following Lemma, whose proof is postponed to Appendix A, shows the relation between $W_{p,t}(\mu, \nu)$ and $\sup_{s \in [0,t]} W_p(\mu_s, \nu_s)$.

**Lemma 2.3.** For every $\mu, \nu \in \mathcal{P}_p(C([0,T], \mathbb{R}^d))$, we have, for every $t \in [0,T]$,

$$ \sup_{s \in [0,t]} W_p(\mu_s, \nu_s) \le W_{p,t}(\mu, \nu), $$

where $\mu_s = \mu \circ \pi_s^{-1}$. In particular, for every $\mu, \nu \in \mathcal{P}_p(C([0,T], \mathbb{R}^d))$, $d_C(\iota(\mu), \iota(\nu)) \le W_p(\mu, \nu)$. Hence the application $\iota$ defined in (2.4) is 1-Lipschitz continuous.

## 2.2 Definition of the particle systems with and without interaction

We define the continuous expansion $(X_t^{1,N}, ..., X_t^{N,N})_{t \in [0,T]}$ of the particle system (1.8) as follows: set $\bar{X}_0^{n,N}, 1 \le n \le N \stackrel{\text{i.i.d.}}{\sim} X_0$; for any $n \in \{1, ..., N\}$ and for any $t \in (t_m, t_{m+1}]$, set

$$ \bar{X}_t^{n,N} = \bar{X}_{t_m}^{n,N} + (t - t_m)b(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) + \sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N)(B_t^n - B_{t_m}^n) \text{ with } \bar{\mu}_{t_m}^N := \frac{1}{N} \sum_{n=1}^N \delta_{\bar{X}_{t_m}^{n,N}}. \quad (2.6) $$

For every $t \in [t_m, t_{m+1})$, we define $\underline{t} = t_m$. Then, the continuous expansion $(X_t^{1,N}, ..., X_t^{N,N})_{t \in [0,T]}$ are solutions to the following equations with initial random variable $\bar{X}_0^{n,N}, 1 \le n \le N$,

$$ \bar{X}_t^{n,N} = \bar{X}_0^{n,N} + \int_0^t b(\underline{u}, \bar{X}_{\underline{u}}^{n,N}, \bar{\mu}_{\underline{u}}^N) du + \int_0^t \sigma(\underline{u}, \bar{X}_{\underline{u}}^{n,N}, \bar{\mu}_{\underline{u}}^N) dB_t^n, \quad 1 \le n \le N, \text{ with } \bar{\mu}_{\underline{t}}^N := \frac{1}{N} \sum_{n=1}^N \delta_{\bar{X}_{\underline{t}}^{n,N}}. \quad (2.7) $$

**Lemma 2.4.** Suppose that Assumption I holds for some $p \in [2, +\infty)$.

(a) The coefficients $b$ and $\sigma$ have a linear growth in the sense that there exists a constant $C_{b,\sigma,L,T}$ depending on $b, \sigma, L$ and $T$ such that

$$ \forall t \in [0,T], \forall x \in \mathbb{R}^d, \forall \mu \in \mathcal{P}_p(\mathbb{R}^d), |b(t, x, \mu)| \vee \|\sigma(t, x, \mu)\| \le C_{b,\sigma,L,T}(1 + |x| + W_p(\mu, \delta_0)). \quad (2.8) $$

(b) Let $(\bar{X}_t^{1,N}, ..., \bar{X}_t^{N,N})_{t \in [0,T]}$ be processes defined by (2.6). Then for a fixed temporal discretization number $M \in \mathbb{N}^*$, we have

$$ \forall n \in \{1, ..., N\}, \quad \left\| \sup_{t \in [0,T]} |\bar{X}_t^{n,N}| \right\|_p < +\infty. \quad (2.9) $$

The proof of Lemma 2.4 is postponed to Appendix A. The system $(\bar{X}_t^{1,N}, ..., \bar{X}_t^{N,N})_{t \in [0,T]}$ defined by (2.6) can be considered as a particle system having interaction through $\bar{\mu}_{t_m}^N = \frac{1}{N} \sum_{n=1}^N \delta_{\bar{X}_{t_m}^{n,N}}, 0 \le m \le M$. Now we define another particle system $(Y_t^1, ..., Y_t^N)_{t \in [0,T]}$ without interaction, which are essentially i.i.d. processes of $\bar{X} = (\bar{X}_t)_{t \in [0,T]}$ defined by the continuous Euler scheme (1.7). Based on the same Brownian motions $B^n, n = 1, ..., N$ as in (2.6), these processes $Y^n = (Y_t^n)_{t \in [0,T]}$, $1 \le n \le N$, are defined by the following equations

$$ \forall n = 1, ..., N, \forall m = 0, ..., M-1, \forall t \in (t_m, t_{m+1}] $$
$$Y_t^n = Y_{t_m}^n + (t - t_m)b(t_m, Y_{t_m}^n, \bar{\mu}_{t_m}) + \sigma(t_m, Y_{t_m}^n, \bar{\mu}_{t_m})(B_t^n - B_{t_m}^n), \quad Y_0^n = \bar{X}_0^{n,N} \quad (2.10)$$

where for every $t \in [0, T]$, $\bar{\mu}_t$ in (2.10) is the marginal distribution at time $t$ of the process $\bar{X} = (\bar{X}_t)_{t \in [0,T]}$ defined by (1.7). By using the same notation $\underline{t}$ as in (2.7), for every $n = 1, \dots, N$, $(Y_t^n)_{t \in [0,T]}$ defined by (2.10) is the solution of

$$Y_t^n = Y_0^n + \int_0^t b(\underline{u}, Y_{\underline{u}}, \bar{\mu}_{\underline{u}}) du + \int_0^t \sigma(\underline{u}, Y_{\underline{u}}, \bar{\mu}_{\underline{u}}) dB_u^n \quad \text{with} \quad Y_0^n = \bar{X}_0^{n,N}.$$

Moreover, it is obvious that the processes $Y^n, 1 \le n \le N$, are i.i.d copies of $\bar{X}$ and then

$$\nu^{N,\omega} := \frac{1}{N} \sum_{n=1}^{N} \delta_{Y^n(\omega)}, \quad \omega \in \Omega \quad (2.11)$$

is the empirical measure of $\bar{\mu} = \mathbb{P} \circ \bar{X}^{-1}$. When there is no ambiguity, we will write $\nu^N$ instead of $\nu^{N,\omega}$. The random measure $\nu^{N,\omega}$ is valued in $\mathcal{P}_p(\mathcal{C}([0,T], \mathbb{R}^d))$. In fact, for every $\omega \in \Omega$, $Y^n(\omega)$ lies in $\mathcal{C}([0,T], \mathbb{R}^d)$ so that $\|Y^n(\omega)\|_{\text{sup}} < +\infty$. Hence, for every $\omega \in \Omega$,

$$\int_{\mathcal{C}([0,T], \mathbb{R}^d)} \|\xi\|_{\text{sup}}^p \nu^{N,\omega}(d\xi) = \frac{1}{N} \sum_{n=1}^{N} \|Y^n(\omega)\|_{\text{sup}}^p < +\infty.$$

As before, we write $\nu_t^{N,\omega} := \nu^{N,\omega} \circ \pi_t^{-1} = \frac{1}{N} \sum_{n=1}^N \delta_{Y_t^n(\omega)}$. Thus, $\iota(\nu^{N,\omega}) = (\nu_t^{N,\omega})_{t \in [0,T]} \in \mathcal{C}([0,T], \mathcal{P}_p(\mathbb{R}^d))$.

## 2.3 Proof of Theorem 1.2

The proof of Theorem 1.2 relies on a variant version of Gronwall's Lemma (see e.g. [Pag18, Lemma 7.3] for the proof) and the following Theorem 2.6 from [FG15], which shows a non-asymptotic upper bound of the convergence rate in the Wasserstein distance of the empirical measures.

**Lemma 2.5** ("À la Gronwall" Lemma). Let $f : [0, T] \to \mathbb{R}_+$ be a Borel, locally bounded, non-negative and non-decreasing function and let $\psi : [0, T] \to \mathbb{R}_+$ be a non-negative non-decreasing function satisfying

$$\forall t \in [0, T], \quad f(t) \le A \int_0^t f(s)ds + B \left( \int_0^t f^2(s)ds \right)^{\frac{1}{2}} + \psi(t),$$

where $A, B$ are two positive real constants. Then, for any $t \in [0, T]$, $f(t) \le 2e^{(2A+B^2)t}\psi(t)$.

**Theorem 2.6.** ([FG15, Theorem 1]) Let $p > 0$ and let $\mu \in \mathcal{P}_q(\mathbb{R}^d)$ for some $q > p$. Let $U^1(\omega), \dots, U^n(\omega), \dots$ be i.i.d random variables with distribution $\mu$. Let $\mu_n^\omega$ denote the empirical measure of $\mu$ defined by

$$\mu_n^\omega := \frac{1}{n} \sum_{i=1}^{n} \delta_{U^i(\omega)}.$$

Then, there exists a real constant $C$ only depending on $p, d, q$ such that, for all $n \ge 1$,

$$\mathbb{E}(W_p^p(\mu_n^\omega, \mu)) \le CM_q^{p/q}(\mu) \times \begin{cases} n^{-1/2} + n^{-(q-p)/q} & \text{if } p > d/2 \text{ and } q \neq 2p \\ n^{-1/2} \log(1+n) + n^{-(q-p)/q} & \text{if } p = d/2 \text{ and } q \neq 2p \\ n^{-p/d} + n^{-(q-p)/q} & \text{if } p \in (0, d/2) \text{ and } q \neq d/(d-p) \end{cases},$$

$$\text{where } M_q(\mu) := \int_{\mathbb{R}^d} |\xi|^q \mu(d\xi).$$

In particular, Theorem 2.6 implies that for $p \ge 2$,

$$\|W_p(\mu_n^\omega, \mu)\|_p \le CM_q^{1/q}(\mu) \times \begin{cases} n^{-1/2p} + n^{-(q-p)/qp} & \text{if } p > d/2 \text{ and } q \neq 2p \\ n^{-1/2p} (\log(1+n))^{1/p} + n^{-(q-p)/qp} & \text{if } p = d/2 \text{ and } q \neq 2p \\ n^{-1/d} + n^{-(q-p)/qp} & \text{if } p \in (0, d/2) \text{ and } q \neq d/(d-p) \end{cases}. \quad (2.12)$$
Moreover, we need the following Lemma, whose proof is postponed to Appendix A, for the proof of Theorem 1.2.

**Lemma 2.7.** Assume that Assumption I holds for an index $p \in [2, +\infty)$. Then the coefficients $b$ and $\sigma$ satisfy the following inequalities

$$ 
\begin{aligned}
& \forall (X, (\mu_t)_{t \in [0,T]}), (Y, (\nu_t)_{t \in [0,T]}) \in L_{C([0,T], \mathbb{R}^d)}^p (\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P}) \times \mathcal{C}([0,T], \mathcal{P}_p(\mathbb{R}^d)), \forall t \in [0,T] \\
& \left\| \sup_{s \in [0,t]} \left| \int_0^s [b(u, X_u, \mu_u) - b(u, Y_u, \nu_u)] du \right| \right\|_p \le L \int_0^t \left[ \|X_u - Y_u\|_p + \|\mathcal{W}_p(\mu_u, \nu_u)\|_p \right] du, \quad (2.13) \\
& \left\| \sup_{s \in [0,t]} \left| \int_0^s [\sigma(u, X_u, \mu_u) - \sigma(u, Y_u, \nu_u)] dB_u \right| \right\|_p \le C_{d,p,L} \left\{ \int_0^t \left[ \|X_u - Y_u\|_p^2 + \|\mathcal{W}_p(\mu_u, \nu_u)\|_p^2 \right] du \right\}^{\frac{1}{2}},
\end{aligned}
$$

where $C_{d,p,L}$ is a constant depending on $d, p, L$.

**Proof of Theorem 1.2.** (a) For every $n \in \{1, ..., N\}$, we have

$$ |Y_t^n - \bar{X}_t^{n,N}| = \left| \int_0^t \left[ b(\underline{u}, Y_{\underline{u}}^n, \bar{\mu}_{\underline{u}}) - b(\underline{u}, \bar{X}_{\underline{u}}^{n,N}, \bar{\mu}_{\underline{u}}^N) \right] du + \int_0^t \left[ \sigma(\underline{u}, Y_{\underline{u}}^n, \bar{\mu}_{\underline{u}}) - \sigma(\underline{u}, \bar{X}_{\underline{u}}^{n,N}, \bar{\mu}_{\underline{u}}^N) \right] dB_u \right| $$

Hence,

$$ 
\begin{aligned}
& \left\| \sup_{s \in [0,t]} |Y_s^n - \bar{X}_s^{n,N}| \right\|_p \\
& = \left\| \sup_{s \in [0,t]} \left| \int_0^s \left[ b(\underline{u}, Y_{\underline{u}}^n, \bar{\mu}_{\underline{u}}) - b(\underline{u}, \bar{X}_{\underline{u}}^{n,N}, \bar{\mu}_{\underline{u}}^N) \right] du + \int_0^s \left[ \sigma(\underline{u}, Y_{\underline{u}}^n, \bar{\mu}_{\underline{u}}) - \sigma(\underline{u}, \bar{X}_{\underline{u}}^{n,N}, \bar{\mu}_{\underline{u}}^N) \right] dB_u \right| \right\|_p \\
& \le \left\| \sup_{s \in [0,t]} \left| \int_0^s \left[ b(\underline{u}, Y_{\underline{u}}^n, \bar{\mu}_{\underline{u}}) - b(\underline{u}, \bar{X}_{\underline{u}}^{n,N}, \bar{\mu}_{\underline{u}}^N) \right] du + \sup_{s \in [0,t]} \left| \int_0^s \left[ \sigma(\underline{u}, Y_{\underline{u}}^n, \bar{\mu}_{\underline{u}}) - \sigma(\underline{u}, \bar{X}_{\underline{u}}^{n,N}, \bar{\mu}_{\underline{u}}^N) \right] dB_u \right| \right| \right\|_p \\
& \le \left\| \sup_{s \in [0,t]} \left| \int_0^s \left[ b(\underline{u}, Y_{\underline{u}}^n, \bar{\mu}_{\underline{u}}) - b(\underline{u}, \bar{X}_{\underline{u}}^{n,N}, \bar{\mu}_{\underline{u}}^N) \right] du \right| \right\|_p \\
& \qquad + \left\| \sup_{s \in [0,t]} \left| \int_0^s \left[ \sigma(\underline{u}, Y_{\underline{u}}^n, \bar{\mu}_{\underline{u}}) - \sigma(\underline{u}, \bar{X}_{\underline{u}}^{n,N}, \bar{\mu}_{\underline{u}}^N) \right] dB_u \right| \right\|_p \quad \text{(by the Minkowski inequality)} \\
& \le L \int_0^t \left[ \|Y_{\underline{u}}^n - \bar{X}_{\underline{u}}^{n,N}\|_p + \|\mathcal{W}_p(\bar{\mu}_{\underline{u}}, \bar{\mu}_{\underline{u}}^N)\|_p \right] du \\
& \qquad + C_{d,p,L} \left[ \int_0^t \left[ \|Y_{\underline{u}}^n - \bar{X}_{\underline{u}}^{n,N}\|_p^2 + \|\mathcal{W}_p(\bar{\mu}_{\underline{u}}, \bar{\mu}_{\underline{u}}^N)\|_p^2 \right] du \right]^{\frac{1}{2}} \quad \text{(by Lemma 2.7)} \\
& \le L \int_0^t \left\| \sup_{v \in [0,u]} |Y_v^n - \bar{X}_v^{n,N}| \right\|_p du + C_{d,p,L} \left[ \int_0^t \left\| \sup_{v \in [0,u]} |Y_v^n - \bar{X}_v^{n,N}| \right\|_p^2 du \right]^{\frac{1}{2}} + \psi(t),
\end{aligned}
$$

where

$$ \psi(t) = L \int_0^t \left\| \mathcal{W}_p(\bar{\mu}_{\underline{u}}, \bar{\mu}_{\underline{u}}^N) \right\|_p du + C_{d,p,L} \left[ \int_0^t \left\| \mathcal{W}_p(\bar{\mu}_{\underline{u}}, \bar{\mu}_{\underline{u}}^N) \right\|_p^2 du \right]^{\frac{1}{2}}, \quad (2.14) $$

owing to $\sqrt{a+b} \le \sqrt{a} + \sqrt{b}$ for any $a \ge 0, b \ge 0$. Then by Lemma 2.5, we have

$$ \left\| \sup_{s \in [0,t]} |Y_s^n - \bar{X}_s^{n,N}| \right\|_p \le 2 e^{(2L+C_{d,p,L}^2)t} \psi(t). $$

Here we can apply Lemma 2.5 since Lemma 2.1 and Lemma 2.4-(b) imply

$$ Y^n, \bar{X}^{n,N} \in L_{C([0,T], \mathbb{R}^d)}^p (\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P}), \quad n = 1, ..., N $$

as $Y^n, n = 1, ..., N \stackrel{\text{i.i.d.}}{\sim} \bar{X}$, thus the process $Y^n - \bar{X}^{n,N}$ is also in $L_{C([0,T], \mathbb{R}^d)}^p (\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P})$ for every
$n = 1, \dots, N$.

Moreover, the empirical measure $\frac{1}{N} \sum_{n=1}^{N} \delta_{(\bar{X}^{n,N}, Y^n)}$ is a coupling of the random measures $\bar{\mu}^N$ and $\nu^N$. Thus

$$
\begin{align*}
\mathbb{E} \mathbb{W}_{p,t}^p (\bar{\mu}^N, \nu^N) &= \mathbb{E} \left[ \inf_{\pi \in \Pi(\bar{\mu}^N, \nu^N)} \int_{C([0,T], \mathbb{R}^d) \times C([0,T], \mathbb{R}^d)} \sup_{s \in [0,t]} |x_s - y_s|^p \pi(dx, dy) \right] \\
&\le \mathbb{E} \left[ \int_{C([0,T], \mathbb{R}^d) \times C([0,T], \mathbb{R}^d)} \sup_{s \in [0,t]} |x_s - y_s|^p \frac{1}{N} \sum_{n=1}^N \delta_{(\bar{X}^{n,N}, Y^n)}(dx, dy) \right] \\
&= \mathbb{E} \left[ \frac{1}{N} \sum_{n=1}^N \sup_{s \in [0,t]} |\bar{X}_s^{n,N} - Y_s^n|^p \right] = \frac{1}{N} \sum_{n=1}^N \left\| \sup_{s \in [0,t]} |\bar{X}_s^{n,N} - Y_s^n| \right\|_p^p \\
&\le \left[ 2 e^{(2L+C_{d,p,L}^2)t} \psi(t) \right]^p \le \left[ 2 e^{(2L+C_{d,p,L}^2)^T} \psi(t) \right]^p.
\end{align*}
$$

Lemma 2.3 implies $\sup_{s \in [0,t]} \mathbb{W}_p^p(\bar{\mu}_s^N, \nu_s^N) \le \mathbb{W}_{p,t}^p(\bar{\mu}^N, \nu^N)$. Hence,

$$
\left\| \sup_{s \in [0,t]} \mathcal{W}_p(\bar{\mu}_s^N, \nu_s^N) \right\|_p \le C_{d,p,L,T} \psi(t) \quad (2.15)
$$

with $C_{d,p,L,T} = 2 e^{(2L+C_{d,p,L}^2)^T}$. It follows that,

$$
\begin{align*}
\left\| \sup_{s \in [0,t]} \mathcal{W}_p(\bar{\mu}_s^N, \bar{\mu}_s) \right\|_p &\le \left\| \sup_{s \in [0,t]} \left[ \mathcal{W}_p(\bar{\mu}_s^N, \nu_s^N) + \mathcal{W}_p(\nu_s^N, \bar{\mu}_s) \right] \right\|_p \\
&\le \left\| \sup_{s \in [0,t]} \mathcal{W}_p(\bar{\mu}_s^N, \nu_s^N) \right\|_p + \left\| \sup_{s \in [0,t]} \mathcal{W}_p(\nu_s^N, \bar{\mu}_s) \right\|_p \\
&\le C_{d,p,L,T} \psi(t) + \left\| \sup_{s \in [0,t]} \mathcal{W}_p(\nu_s^N, \bar{\mu}_s) \right\|_p \quad \text{(by applying (2.15))} \\
&\le \left\| \sup_{s \in [0,t]} \mathcal{W}_p(\nu_s^N, \bar{\mu}_s) \right\|_p + C_{d,p,L,T} \cdot L \int_0^t \left\| \mathcal{W}_p(\bar{\mu}_{\underline{u}}, \bar{\mu}_{\underline{u}}^N) \right\|_p du + C_{d,p,L,T} \cdot C_{d,p,L} \left[ \int_0^t \left\| \mathcal{W}_p(\bar{\mu}_{\underline{u}}, \bar{\mu}_{\underline{u}}^N) \right\|_p^2 du \right]^{\frac{1}{2}} \\
&\qquad \text{(by the definition of } \psi(t) \text{ in (2.14))} \\
&\le \left\| \sup_{s \in [0,t]} \mathcal{W}_p(\nu_s^N, \bar{\mu}_s) \right\|_p + C_{d,p,L,T} \cdot L \int_0^t \left\| \sup_{v \in [0,u]} \mathcal{W}_p(\bar{\mu}_v, \bar{\mu}_v^N) \right\|_p du \\
&\qquad + C_{d,p,L,T} \cdot C_{d,p,L} \left[ \int_0^t \left\| \sup_{v \in [0,u]} \mathcal{W}_p(\bar{\mu}_v, \bar{\mu}_v^N) \right\|_p^2 du \right]^{\frac{1}{2}}.
\end{align*}
$$

Then, by Lemma 2.5, we obtain

$$
\left\| \sup_{s \in [0,t]} \mathcal{W}_p(\bar{\mu}_s^N, \bar{\mu}_s) \right\|_p \le 2e^{(2A+B^2)T} \left\| \sup_{s \in [0,t]} \mathcal{W}_p(\bar{\mu}_s, \nu_s^N) \right\|_p, \quad (2.16)
$$

where $A = C_{d,p,L,T}L$ and $B = C_{d,p,L,T} \cdot C_{d,p,L}$. Finally,

$$
\begin{align*}
\left\| \sup_{0 \le m \le M} \mathcal{W}_p(\bar{\mu}_{t_m}^N, \bar{\mu}_m) \right\|_p &\le 2e^{(2A+B^2)T} \left\| \sup_{s \in [0,T]} \mathcal{W}_p(\bar{\mu}_s, \nu_s^N) \right\|_p \\
&\le 2e^{(2A+B^2)T} \left\| \mathbb{W}_p(\bar{\mu}, \nu^N) \right\|_p \to 0 \quad \text{as } N \to +\infty,
\end{align*}
$$

where the second inequality above follows from Lemma 2.3. For the convergence of $\|\mathbb{W}_p(\bar{\mu}, \nu^N)\|_p \to 0$ as $N \to +\infty$, we refer to [Lac18, Corollary 2.14] among many other references.

(b) If $\|X_0\|_{p+\varepsilon} < +\infty$ for some $\varepsilon > 0$, the Assumption I holds with index $p + \varepsilon$ since for every $\mu, \nu \in \mathcal{P}_p(\mathbb{R}^d)$, $\mathcal{W}_p(\mu, \nu) \le \mathcal{W}_{p+\varepsilon}(\mu, \nu)$. Then, the inequality (2.2) implies

$$
\| \bar{X} \|_{p+\varepsilon} = \left\| \sup_{u \in [0,T]} |\bar{X}_u| \right\|_{p+\varepsilon} \le C_{p,d,b,\sigma,T,L} (1 + \|X_0\|_{p+\varepsilon}) < +\infty.
$$

Thus $\bar{\mu} \in \mathcal{P}_{p+\varepsilon}(\mathcal{C}([0,T], \mathbb{R}^d))$, which implies that $\bar{\mu}_s \in \mathcal{P}_{p+\varepsilon}(\mathbb{R}^d)$ for any $s \in [0,T]$.
For any $s \in [0, T]$, $\nu_s^N$ is the empirical measure of $\bar{\mu}_s$. It follows from Theorem 2.6 that for any $s \in [0, T]$,

$$ \begin{aligned} \|W_p(\nu_s^N, \bar{\mu}_s)\|_p &\le CM_{p+\varepsilon}^{1/p+\varepsilon}(\bar{\mu}_s) \\ &\quad \times \begin{cases} N^{-1/2p} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p > d/2 \text{ and } \varepsilon \ne p \\ N^{-1/2p} (\log(1+N))^{1/p} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p = d/2 \text{ and } \varepsilon \ne p \\ N^{-1/d} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p \in (0, d/2) \text{ and } p + \varepsilon \ne \frac{d}{(d-p)} \end{cases} \end{aligned} \quad (2.17) $$

where $C$ is a constant depending on $p, \varepsilon, d$. Moreover, the inequality (2.2) implies that

$$ \begin{aligned} \sup_{s \in [0,T]} M_{p+\varepsilon}(\bar{\mu}_s) &= \sup_{s \in [0,T]} \mathbb{E}[|X_s|^{p+\varepsilon}] \le \left\| \sup_{s \in [0,T]} |X_s| \right\|_{p+\varepsilon}^{p+\varepsilon} \\ &\le \left[ C_{p,d,b,\sigma,L,T} (1 + \|X_0\|_{p+\varepsilon}) \right]^{p+\varepsilon} < +\infty. \end{aligned} $$

Hence,

$$ \sup_{s \in [0,T]} \|W_p(\nu_s^N, \bar{\mu}_s)\|_p \le C' \times \begin{cases} N^{-1/2p} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p > d/2 \text{ and } \varepsilon \ne p \\ N^{-1/2p} (\log(1+N))^{1/p} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p = d/2 \text{ and } \varepsilon \ne p \\ N^{-1/d} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p \in (0, d/2) \text{ and } p + \varepsilon \ne \frac{d}{(d-p)} \end{cases}, \quad (2.18) $$

where $C' = C \sup_{s \in [0,T]} M_{p+\varepsilon}^{\frac{1}{p+\varepsilon}}(\bar{\mu}_s)$ which is a constant depending on $p, \varepsilon, d, b, \sigma, L, T$ and $\|X_0\|_{p+\varepsilon}$.
Moreover, the inequality (2.15) implies that

$$ \sup_{s \in [0,t]} \|W_p(\bar{\mu}_s^N, \nu_s^N)\|_p \le \left\| \sup_{s \in [0,t]} W_p(\bar{\mu}_s^N, \nu_s^N) \right\|_p \le C_{d,p,L,T} \psi(t). \quad (2.19) $$

Then,

$$ \begin{aligned} \sup_{s \in [0,t]} \|W_p(\bar{\mu}_s^N, \bar{\mu}_s)\|_p &\le \sup_{s \in [0,t]} \|W_p(\bar{\mu}_s^N, \nu_s^N) + W_p(\nu_s^N, \bar{\mu}_s)\|_p \\ &\le \sup_{s \in [0,t]} \|W_p(\bar{\mu}_s^N, \nu_s^N)\|_p + \sup_{s \in [0,t]} \|W_p(\nu_s^N, \bar{\mu}_s)\|_p \\ &\le C_{d,p,L,T} \psi(t) + \sup_{s \in [0,t]} \|W_p(\nu_s^N, \bar{\mu}_s)\|_p \quad (\text{by applying (2.19)}) \\ &\le \sup_{s \in [0,t]} \|W_p(\nu_s^N, \bar{\mu}_s)\|_p + C_{d,p,L,T} \cdot L \int_0^t \|W_p(\bar{\mu}_{\underline{u}}, \bar{\mu}_{\underline{u}}^N)\|_p du \\ &\quad + C_{d,p,L,T} \cdot C_{d,p,L} \left[ \int_0^t \|W_p(\bar{\mu}_{\underline{u}}, \bar{\mu}_{\underline{u}}^N)\|_p^2 du \right]^{\frac{1}{2}} \quad (\text{by the definition of } \psi(t) \text{ in (2.14)}) \\ &\le \sup_{s \in [0,t]} \|W_p(\nu_s^N, \bar{\mu}_s)\|_p + C_{d,p,L,T} \cdot L \int_0^t \sup_{v \in [0,u]} \|W_p(\bar{\mu}_v, \bar{\mu}_v^N)\|_p du \\ &\quad + C_{d,p,L,T} \cdot C_{d,p,L} \left[ \int_0^t \sup_{v \in [0,u]} \|W_p(\bar{\mu}_v, \bar{\mu}_v^N)\|_p^2 du \right]^{\frac{1}{2}}, \end{aligned} $$

Then, by Lemma 2.5, we obtain

$$ \sup_{s \in [0,T]} \|W_p(\bar{\mu}_s^N, \bar{\mu}_s)\|_p \le 2e^{(2A+B^2)T} \sup_{s \in [0,T]} \|W_p(\bar{\mu}_s, \nu_s^N)\|_p, \quad (2.20) $$

where $A = C_{d,p,L,T}L$ and $B = C_{d,p,L,T} \cdot C_{d,p,L}$. Finally, it follows from (2.18) that

$$ \sup_{0 \le m \le M} \|W_p(\bar{\mu}_{t_m}^N, \bar{\mu}_m)\|_p \le \sup_{s \in [0,T]} \|W_p(\bar{\mu}_s^N, \bar{\mu}_s)\|_p $$
$$ \le \tilde{C} \times \begin{cases} N^{-\frac{1}{2p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p > d/2 \text{ and } \varepsilon \ne p \\ N^{-\frac{1}{2p}} \left[ \log(1+N) \right]^{\frac{1}{p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p = d/2 \text{ and } \varepsilon \ne p \\ N^{-\frac{1}{d}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p \in (0, d/2) \text{ and } p+\varepsilon \ne \frac{d}{(d-p)} \end{cases}, \quad (2.21) $$

where $\tilde{C}$ is a constant depending on $p, \varepsilon, d, b, \sigma, L, T$ and $\|X_0\|_{p+\varepsilon}$.

By the definition of $\psi(t)$ in (2.14), we have

$$ \begin{aligned} \psi(t) &= L \int_0^t \|W_p(\bar{\mu}_u, \bar{\mu}_u^N)\|_p \, du + C_{d,p,L} \left[ \int_0^t \|W_p(\bar{\mu}_u, \bar{\mu}_u^N)\|_p^2 \, du \right]^{\frac{1}{2}}, \\ &\le L T \sup_{u \in [0,T]} \|W_p(\bar{\mu}_u, \bar{\mu}_u^N)\|_p + C_{d,p,L} \sqrt{T} \cdot \sup_{u \in [0,T]} \|W_p(\bar{\mu}_u, \bar{\mu}_u^N)\|_p \\ &\le C_{d,p,L,T} \sup_{u \in [0,T]} \|W_p(\bar{\mu}_u, \bar{\mu}_u^N)\|_p \end{aligned} $$

so that (2.15) and (2.21) imply that

$$ \begin{aligned} \left\| \sup_{s \in [0,T]} W_p(\bar{\mu}_s^N, \nu_s^N) \right\|_p &\le C_{d,p,L,T} \sup_{u \in [0,T]} \left\| W_p(\bar{\mu}_u, \bar{\mu}_u^N) \right\|_p \\ &\le C \times \begin{cases} N^{-\frac{1}{2p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p > d/2 \text{ and } \varepsilon \ne p \\ N^{-\frac{1}{2p}} \left[ \log(1+N) \right]^{\frac{1}{p}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p = d/2 \text{ and } \varepsilon \ne p \\ N^{-\frac{1}{d}} + N^{-\frac{\varepsilon}{p(p+\varepsilon)}} & \text{if } p \in (0, d/2) \text{ and } p+\varepsilon \ne \frac{d}{(d-p)} \end{cases} \end{aligned} $$

where $C$ is a constant depending on $p, \varepsilon, d, b, \sigma, L, T$ and $\|X_0\|_{p+\varepsilon}$. $\square$

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

# 3 Quantization-based schemes and their error analyses

This section is devoted to the error analyses of the quantization based schemes (1.16), (1.18) and (1.22). We start with a review of optimal quantization, as well as its connection with the K-means clustering. We also prove Theorem 1.3, the $L^2$-error analysis of the quantization based scheme (1.16) in Subsection 3.1. Next in Subsection 3.2, we give computational details of transition probabilities (1.18) for the recursive quantization scheme and show how to optimally compute the quantizer $x^{(m)}$ by integrating Lloyd's algorithm to reduce the quantization error. Finally, in Subsection 3.3, we prove Proposition 1.4, the error analysis of the hybrid particle-quantization scheme.

## 3.1 Preliminaries on the optimal quantization

In this section, $X: (\Omega, \mathcal{F}, \mathbb{P}) \to (\mathbb{R}^d, \mathcal{B}(\mathbb{R}^d))$ is a random variable having $p$-th finite moment, $p \ge 1$. Its probability distribution is denoted by $\mu$. Let $K \in \mathbb{N}^*$ be the fixed quantization level. We call $x = (x_1, ..., x_K) \in (\mathbb{R}^d)^K$ a quantizer of level $K$. The quantization error function of $\mu$ at level $K$ and order $p$, denoted by $e_{K,p}(\mu, \cdot)$, is defined by

$$ x = (x_1, ..., x_K) \in (\mathbb{R}^d)^K \mapsto e_{K,p}(\mu, x) := \left[ \int_{\mathbb{R}^d} \min_{1 \le k \le K} |\xi - x_k|^p \mu(d\xi) \right]^{\frac{1}{p}} \in \mathbb{R}_+. \quad (3.1) $$

Sometimes we write $e_{K,p}(X, \cdot)$ for the quantization error function since $P_X = \mu$, but the value of this function depends only on the probability distribution $\mu$. For a quantizer $x = (x_1, ..., x_K)$ such that $x_i \ne x_j, i \ne j$ and for $(V_k(x))_{1 \le k \le K}$ a Voronoï partition generated by $x$ defined in (1.12), the quantization error of $x$ satisfying

$$ e_{K,p}^p(\mu, x) = \sum_{k=1}^{K} \int_{V_k(x)} |\xi - x_k|^p \mu(d\xi). \quad (3.2) $$
It is obvious that the error value does not depend on the choice of Voronoï partition. Moreover, still considering this Voronoï partition $(V_k(x))_{1 \le k \le K}$, we have (see e.g. [GL00, Lemma 3.4] <fnref n="3" />)

$$e_{K,p}(\mu, x) = \|X - \hat{X}^x\|_p = W_p(\mu, \hat{\mu}^x) \quad (3.3)$$

where $\hat{X}^x := \text{Proj}_x(X)$ and $\hat{\mu}^x := \mu \circ \text{Proj}_x^{-1} = \mathcal{L}(\hat{X}^x)$ with the projection function $\text{Proj}_x$ defined by (1.13).

The optimal quantizer at level K and order p for $\mu$, denoted by $x^* = (x_1^*, ..., x_K^*)$, is defined by

$$x^* = (x_1^*, ..., x_K^*) \in \underset{y \in (\mathbb{R}^d)^K}{\operatorname{argmin}} e_{K,p}(\mu, y). \quad (3.4)$$

The existence of such an optimal quantizer has been proved in [Pol82]. There exists few results in the literature for the uniqueness of the optimal quantizer (see e.g. [Kie83]) but in practice, we only need one optimal quantizer for a further simulation as all optimal quantizers have the same quantization error by (3.4). We call

$$e_{K,p}^*(X) = e_{K,p}^*(\mu) := \inf_{x \in (\mathbb{R}^d)^K} e_{K,p}(\mu, x)$$

the optimal (quantization) error at level K and order p for $\mu$ (or for X).

Now we recall several basic properties of the optimal quantization.

**Proposition 3.1.** Let $K \in \mathbb{N}^*$ be the quantization level and let $p \ge 1$. We consider an $\mathbb{R}^d$-valued random variable X having probability distribution $\mu \in \mathcal{P}_p(\mathbb{R}^d)$. Let $x^* = (x_1^*, ..., x_K^*)$ be an optimal quantizer of X at level K and order p.

(1) (see [GL00, Theorem 4.1 and Theorem 4.12]) Assume that the support of $\mu$ satisfies $\text{card}(\text{supp}(\mu)) \ge K$, where $\text{card}(\text{supp}(\mu))$ denotes the cardinality of the support of $\mu$. Then any K-level optimal quantizer $x^*$ contains K different points, i.e. for every $i, j \in \{1, ..., K\}$ such that $i \ne j$, we have $x_i^* \ne x_j^*$.

(2) (Zador's theorem, see [LP08] and [Pag18, Theorem 5.2]) If $\mu \in \mathcal{P}_{p+\epsilon}(\mathbb{R}^d)$ for some $\epsilon > 0$, the optimal error has the following non-asymptotic upper bound for every quantization level $K \ge 1$,

$$e_{p,K}^*(X) = e_{p,K}^*(\mu) \le C_{d,p,\epsilon} \mathcal{M}_{p+\epsilon}(\mu) K^{-\frac{1}{d}}, \quad (3.5)$$

where $C_{d,p,\epsilon} \ge 0$ is a constant depending only on d, p and $\epsilon$, and where

$$\mathcal{M}_p(\mu) = \mathcal{M}_p(X) := \inf_{a \in \mathbb{R}^d} \|X - a\|_p. \quad (3.6)$$

(3) (Optimal discrete representation, see [GL00, Lemma 3.4] and [Pag18, Proposition 5.1-(b)]) Let $\hat{X}^x = \text{Proj}_{x^*}(X)$ and let $\hat{\mu}^x := \mu \circ \text{Proj}_{x^*}^{-1} = \mathcal{L}(\hat{X}^x)$. Then

$$\|X - \hat{X}^{x^*}\|_p = \min \left\{ \|X - Y\|_p \mid Y : (\Omega, \mathcal{F}, \mathcal{P}) \to (\mathbb{R}^d, \mathcal{B}(\mathbb{R}^d)) \text{ random variable such that} \right. \\ \text{the support of } Y \text{ has at most } K \text{ points} \right\}$$

$$\text{and } W_p(\mu, \hat{\mu}^{x^*}) = \min \left\{ W_p(\mu, \nu) \mid \nu \in \mathcal{P}_p(\mathbb{R}^d) \text{ such that } \text{card}(\text{supp}(\nu)) \le K \right\}.$$

(4) (Stationary property when $p=2$, see e.g. [Pag18, Proposition 5]) For the quadratic optimal quantization (i.e. $p=2$), any optimal quantizer $x^*$ is stationary in the sense that

$$\mathbb{E}[X | \hat{X}^{x^*}] = \hat{X}^{x^*}. \quad (3.7)$$

(5) (Consistency of the optimal quantization, see [Pol82, Theorem 9] and [LP20, Theorem 4 and Appendix A]) Consider a probability distribution sequence $\mu_n$, $\mu \in \mathcal{P}_2(\mathbb{R}^d)$, $n \ge 1$, such that $W_2(\mu_n, \mu) \to$

<footnote n="3">³The statement in [GL00, Lemma 3.4] is established only for an optimal quantizer. However, its proof is also valid for an arbitrary quantizer.</footnote>
0 when $n \to 0$. For every $n \ge 1$, let $x^{(n)}$ be a quadratic optimal quantizer of $\mu_n$. Then any limiting point of $(x^{(n)})_{n \ge 1}$ is a quadratic optimal quantizer of $\mu$. Moreover, we have

$$ (e_{K,2}^*(\mu_n))^2 - (e_{K,2}^*(\mu))^2 \le 4e_{K,2}^*(\mu) W_2(\mu_n, \mu) + 4W_2^2(\mu_n, \mu). \quad (3.8) $$

From a numerical point of view, the main idea of the optimal quantization is to use $\hat{X}^{x^*} = \text{Proj}_{x^*}(X)$ (respectively, $\hat{\mu}^{x^*} = \mu \circ \text{Proj}_{x^*}^{-1} = \mathcal{L}(\hat{X}^{x^*})$) as an approximation of the target random variable $X$ (resp. the target probability distribution $\mu$). Proposition 3.1-(3) shows that $\hat{X}^{x^*}$ (respectively $\hat{\mu}^{x^*}$) is the closest discrete representation of $X$ (resp. of $\mu$) with respect to the $L^p$-norm (resp. the Wasserstein distance $W_p$), among all random variables (resp. probability distributions) having a support of at most $K$ points. Furthermore, the Proposition 3.1-(5) shows the connection between the optimal quantization and the $K$-means clustering in unsupervised learning (see e.g. [Pol82] and [LP20]). Consider a sample $\{\eta_1, ..., \eta_n\} \subset \mathbb{R}^d$, the $K$-means clustering is essentially to find an optimal quantizer with respect to the empirical measure $\bar{\mu}_n := \frac{1}{n} \sum_{i=1}^n \delta_{\eta_i}$ over the sample. In the unsupervised learning context, such an optimal quantizer is called cluster centers and $\eta_1, ..., \eta_n$ are often considered as i.i.d. samples having probability distribution $\mu$ so that $W_p(\bar{\mu}_n, \mu) \to 0$ a.s. if $\mu \in \mathcal{P}_p(\mathbb{R}^d)$ (see e.g. [Pol82, Thereom 7]).

Finally, we recall that there exist several numerical methods to find an optimal quantizer in the quadratic setting ($p=2$) such as the (stochastic) gradient descent algorithm (see e.g. [Pag15, Section 3.2]) and Lloyd's fixed point algorithm (see e.g. [Llo82] and [PY16] for its convergence). In this paper, we implement Lloyd's algorithm to obtain the quadratic optimal quantizer in the simulations but numerical results can be obtained by other methods as well. Lloyd's algorithm, whose objective is to find a stationary quantizer in the sense of (3.7), is described as follows.

<caption>**Algorithm 1:** Lloyd's (fixed point) algorithm for a probability measure $\mu \in \mathcal{P}_2(\mathbb{R}^d)$</caption>

Set $K \in \mathbb{N}^*$ for the quantization level.

**Input:** Initial quantizer $x^{[0]} = (x_1^{[0]}, ..., x_K^{[0]})$ such that $x_k^{[0]} \in \text{supp}(\mu)$, $1 \le k \le K$.

**Repeat**

$$ \text{Compute } x^{[l+1]} = \mathbb{E}[X | \hat{X}^{x^{[l]}}](\Omega), \text{ i.e. for every } k = 1, ..., K, \\ x_k^{[l+1]} = \begin{cases} x_k^{[l]}, & \text{if } \mu(V_k(x^{[l]})) = 0, \\ \mathbb{E}[X | \hat{X}^{x^{[l]}} = x_k^{[l]}] = \frac{\int_{V_k(x^{[l]})} \xi \mu(d\xi)}{\mu(V_k(x^{[l]}))}, & \text{otherwise.} \end{cases} \quad (3.9) $$

**until** Some stopping criterion occurs (for example, $|x^{[l+1]} - x^{[l]}| < \varepsilon$ for a fixed threshold $\varepsilon > 0$);

**Return** $x^{[l]} = (x_1^{[l]}, ..., x_K^{[l]})$

Now we give the proof of Theorem 1.3.

*Proof of Theorem 1.3.* (1) To simplify the notation, we will denote by

$$ \bar{b}_m := b(t_m, \bar{X}_{t_m}, \bar{\mu}_{t_m}), \quad \bar{\sigma}_m := \sigma(t_m, \bar{X}_{t_m}, \bar{\mu}_{t_m}), \quad \hat{b}_m := b(t_m, \hat{X}_{t_m}, \hat{\mu}_{t_m}), \quad \hat{\sigma}_m := \sigma(t_m, \hat{X}_{t_m}, \hat{\mu}_{t_m}). $$

The definitions of $\bar{X}_{t_m}$, $\tilde{X}_{t_m}$ and in (1.6) and (1.16) directly imply that

$$ |\bar{X}_{t_{m+1}} - \tilde{X}_{t_{m+1}}| = |(\bar{X}_{t_m} - \hat{X}_{t_m}) + h(\bar{b}_m - \hat{b}_m) + \sqrt{h}(\bar{\sigma}_m - \hat{\sigma}_m)Z_{m+1}|. $$

Hence,

$$ \begin{aligned} & \mathbb{E}[|\bar{X}_{t_{m+1}} - \tilde{X}_{t_{m+1}}|^2] \\ &= \mathbb{E}[|(\bar{X}_{t_m} - \hat{X}_{t_m}) + h(\bar{b}_m - \hat{b}_m)|^2] + h\mathbb{E}[|(\bar{\sigma}_m - \hat{\sigma}_m)Z_{m+1}|^2] \\ & \quad + 2\sqrt{h}\mathbb{E}[\langle (\bar{X}_{t_m} - \hat{X}_{t_m}) + h(\bar{b}_m - \hat{b}_m) | (\bar{\sigma}_m - \hat{\sigma}_m)Z_{m+1} \rangle] \end{aligned} $$
$$ = \mathbb{E} \left[ |\bar{X}_{t_m} - \hat{X}_{t_m}|^2 \right] + h^2 \mathbb{E} \left[ |\bar{b}_m - \hat{b}_m|^2 \right] + 2h \mathbb{E} \left[ \langle \bar{X}_{t_m} - \hat{X}_{t_m} | \bar{b}_m - \hat{b}_m \rangle \right] \\ + h \mathbb{E} \left[ |(\bar{\sigma}_m - \hat{\sigma}_m) Z_{m+1}|^2 \right] + 2\sqrt{h} \mathbb{E} \left[ \langle (\bar{X}_{t_m} - \hat{X}_{t_m}) + h(\bar{b}_m - \hat{b}_m) | (\bar{\sigma}_m - \hat{\sigma}_m) Z_{m+1} \rangle \right]. \quad (3.10) $$

Let $\mathcal{F}_0$ denote the $\sigma$-algebra generated by $X_0$. For every $m \in \{1, ..., M\}$, we define $\mathcal{F}_m$ the $\sigma$-algebra generated by $X_0, Z_1, ..., Z_m$. Then, as $Z_{m+1}$ is independent of $\mathcal{F}_m$ and $\bar{X}_{t_m}, \hat{X}_{t_m}, \bar{b}_m, \hat{b}_m, \bar{\sigma}_m, \hat{\sigma}_m$ are $\mathcal{F}_m$-measurable, we have

$$ \begin{aligned} & \mathbb{E} \left[ \left\langle (\bar{X}_{t_m} - \hat{X}_{t_m}) + h(\bar{b}_m - \hat{b}_m) \mid (\bar{\sigma}_m - \hat{\sigma}_m) Z_{m+1} \right\rangle \right] \\ &= \mathbb{E} \left[ \left[ (\bar{X}_{t_m} - \hat{X}_{t_m}) + h(\bar{b}_m - \hat{b}_m) \right]^\top (\bar{\sigma}_m - \hat{\sigma}_m) Z_{m+1} \right] \\ &= \mathbb{E} \left\{ \mathbb{E} \left[ \left[ (\bar{X}_{t_m} - \hat{X}_{t_m}) + h(\bar{b}_m - \hat{b}_m) \right]^\top (\bar{\sigma}_m - \hat{\sigma}_m) Z_{m+1} \mid \mathcal{F}_m \right] \right\} \\ &= \mathbb{E} \left\{ \left[ \left[ (\bar{X}_{t_m} - \hat{X}_{t_m}) + h(\bar{b}_m - \hat{b}_m) \right]^\top (\bar{\sigma}_m - \hat{\sigma}_m) \right] \mathbb{E} \left[ Z_{m+1} \right] \right\} = 0. \end{aligned} $$

Moreover, Assumption I implies that

$$ \mathbb{E} \left[ |\bar{b}_m - \hat{b}_m|^2 \right] \le 2L^2 \left[ \mathbb{E} \left[ |\bar{X}_{t_m} - \hat{X}_{t_m}|^2 \right] + W_2^2 (\bar{\mu}_m, \hat{\mu}_m) \right] \le 4L^2 \mathbb{E} \left[ |\bar{X}_{t_m} - \hat{X}_{t_m}|^2 \right] $$

and

$$ \mathbb{E} \left[ \left\langle \bar{X}_{t_m} - \hat{X}_{t_m} \mid \bar{b}_m - \hat{b}_m \right\rangle \right] \le \left\| \bar{X}_{t_m} - \hat{X}_{t_m} \right\|_2 \left\| \bar{b}_m - \hat{b}_m \right\|_2 \le 2L \mathbb{E} \left[ \left\| \bar{X}_{t_m} - \hat{X}_{t_m} \right\|_2 \right], $$

as well as

$$ \begin{aligned} \mathbb{E} \left[ \left| (\bar{\sigma}_m - \hat{\sigma}_m) Z_{m+1} \right|^2 \right] &\le \mathbb{E} \left[ \left\| \bar{\sigma}_m - \hat{\sigma}_m \right\|^2 | Z_{m+1} \right]^2 \\ &\le \mathbb{E} \left[ \mathbb{E} \left[ \left\| \bar{\sigma}_m - \hat{\sigma}_m \right\|^2 | Z_{m+1} \right]^2 \mid \mathcal{F}_m \right] \\ &= \mathbb{E} \left[ \left\| \bar{\sigma}_m - \hat{\sigma}_m \right\|^2 \mathbb{E} \left[ | Z_{m+1} \right]^2 \right] \\ &\le 4L^2 q \mathbb{E} \left[ \left\| \bar{X}_{t_m} - \hat{X}_{t_m} \right\|_2 \right], \end{aligned} $$

where we recall that $\|\cdot\|$ denotes the operator norm defined by $\|A\| := \sup_{|z| \le 1} |Az|$. Consequently,

$$ \mathbb{E} \left[ \left| \bar{X}_{t_{m+1}} - \tilde{X}_{t_{m+1}} \right|^2 \right] \le \left[ 1 + 4Lh(1 + Lh + Lq) \right] \mathbb{E} \left[ \left| \bar{X}_{t_m} - \hat{X}_{t_m} \right|^2 \right] $$

so that

$$ \begin{aligned} \left\| \bar{X}_{t_{m+1}} - \tilde{X}_{t_{m+1}} \right\|_2 &\le \sqrt{1 + 4Lh(1 + Lh + Lq)} \left\| \bar{X}_{t_m} - \hat{X}_{t_m} \right\|_2 \\ &\le \left[ 1 + 2Lh(1 + Lh + Lq) \right] \left\| \bar{X}_{t_m} - \hat{X}_{t_m} \right\|_2 \end{aligned} $$

where in the last inequality, we use the fact that $\sqrt{1+2x} \le 1+x$ when $x \ge 0$. Hence,

$$ \begin{aligned} \left\| \bar{X}_{t_{m+1}} - \hat{X}_{t_{m+1}} \right\|_2 &\le \left\| \bar{X}_{t_{m+1}} - \tilde{X}_{t_{m+1}} \right\|_2 + \left\| \tilde{X}_{t_{m+1}} - \hat{X}_{t_{m+1}} \right\|_2 \\ &\le \left[ 1 + 2Lh(1 + Lh + Lq) \right] \left\| \bar{X}_{t_m} - \hat{X}_{t_m} \right\|_2 + \Xi_{m+1}, \end{aligned} $$

where $\Xi_m = \|\tilde{X}_{t_m} - \hat{X}_{t_m}\|_2$ denotes the quadratic quantization error at time $t_m$ (see (3.3)). This directly implies

$$ \left\| \bar{X}_{t_m} - \hat{X}_{t_m} \right\|_2 \le \sum_{j=0}^{m} \left[ 1 + 2Lh(1 + Lh + Lq) \right]^j \Xi_{m-j}. $$

(2) STEP 1. We first prove that there exists a constant $\tilde{C}_{\max} > 0$ such that

$$ \forall m \in \{0, ..., M\}, \quad \|\tilde{X}_{t_m}\|_{2+\epsilon} \le \tilde{C}_{\max}. \quad (3.11) $$

As for every $m \in \{0, ..., M\}$, the quantizer $x^{(m)} = (x_1^{(m)}, ..., x_K^{(m)})$ is quadratic optimal for $\tilde{X}_{t_m}$, Propo-
sition 3.1-(4) implies that $\mathbb{E}[\tilde{X}_{t_m} | \hat{X}_{t_m}] = \hat{X}_{t_m}$. Hence, for every $p \ge 1$,

$$ \mathbb{E}[|\tilde{X}_{t_m}|^p] = \mathbb{E}\left[\left|\mathbb{E}[\tilde{X}_{t_m} | \hat{X}_{t_m}]\right|^p\right] \le \mathbb{E}\left[\mathbb{E}[|\tilde{X}_{t_m}|^p | \hat{X}_{t_m}]\right] = \mathbb{E}[|\tilde{X}_{t_m}|^p]. \quad (3.12) $$

On the other hand, we have

$$ \begin{aligned} \|\tilde{X}_{t_{m+1}}\|_{p} &\le \|\hat{X}_{t_m}\|_{p} + h\|b(t_m, \hat{X}_{t_m}, \hat{\mu}_{t_m})\|_{p} + \sqrt{h}\left\|\left\|\sigma(t_m, \hat{X}_{t_m}, \hat{\mu}_{t_m})\right\| Z_{m+1}\right\|_{p} \\ &\le \|\hat{X}_{t_m}\|_{p} + h\|b(t_m, \hat{X}_{t_m}, \hat{\mu}_{t_m})\|_{p} + \sqrt{h}\left\|\left\|\sigma(t_m, \hat{X}_{t_m}, \hat{\mu}_{t_m})\right\| Z_{m+1}\right\|_{p} \\ &\le (1 + 2hC_{b,\sigma,L,T} + 2\sqrt{h}\|Z_{m+1}\|_{p} C_{b,\sigma,L,T}) \|\hat{X}_{t_m}\|_{p} + hC_{b,\sigma,L,T} + \sqrt{h}\|Z_{m+1}\|_{p} C_{b,\sigma,L,T} \end{aligned} $$

where the first inequality follows from the definition of $\tilde{X}_{t_m}$ in (1.16) and Minkowski's inequality, the second inequality follows from the fact that $\sigma(t_m, \hat{X}_{t_m}, \hat{\mu}_{t_m}) \perp Z_{m+1}$ and the third inequality follows from Lemma 2.4-(a) and the inequality $W_p(\hat{\mu}_{t_m}, \delta_0) \le \|\tilde{X}_{t_m}\|_p$. Thus, the inequality (3.12) implies that

$$ \|\tilde{X}_{t_{m+1}}\|_{p} \le C_{1} \|\tilde{X}_{t_m}\|_{p} + C_{2} $$

with $C_1 := 1 + 2hC_{b,\sigma,L,T} + 2\sqrt{h}\|Z_{m+1}\|_p C_{b,\sigma,L,T}$ and $C_2 := hC_{b,\sigma,L,T} + \sqrt{h}\|Z_{m+1}\|_p C_{b,\sigma,L,T}$. As $\tilde{X}_0 = X_0$, we have

$$ \forall m \in \{1, \dots, M\}, \quad \|\tilde{X}_{t_m}\|_p \le C_1^m \|X_0\|_p + \sum_{j=0}^{m-1} C_1^j C_2 \quad (3.13) $$

Thus (3.11) is obtained by considering $p = 2 + \varepsilon$ and by letting $\tilde{C}_{\max} := C_1^M \|X_0\|_{2+\varepsilon} + \sum_{j=0}^{M-1} C_1^j C_2$.

STEP 2. As for every $m \in \{0, \dots, M\}$, the quantizer $x^{(m)} = (x_1^{(m)}, \dots, x_K^{(m)})$ is quadratic optimal for $\tilde{X}_{t_m}$, Zador's theorem (Proposition 3.1-(2)) implies that

$$ \forall m \in \{0, \dots, M\}, \quad \Xi_m \le C_{d,\varepsilon} \mathcal{M}_{2+\varepsilon}(\tilde{X}_{t_m}) K^{-\frac{1}{d}}. $$

Then we can conclude by remarking $\mathcal{M}_{2+\varepsilon}(\tilde{X}_{t_m}) \le \|\tilde{X}_{t_m}\|_{2+\varepsilon} \le \tilde{C}_{\max}$. $\square$

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

## 3.2 Recursive quantization-based scheme for the Vlasov equation

In this section, we present the recursive quantization scheme in the Vlasov setting and show how to integrate Lloyd's algorithm into this scheme to reduce the quantization error at each time step.

Recall that in the Vlasov setting, there exist $\beta : [0, T] \times \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d$ and $a : [0, T] \times \mathbb{R}^d \times \mathbb{R}^d \to \mathcal{M}_{d,q}(\mathbb{R})$ such that

$$ b(t, x, \mu) = \int_{\mathbb{R}^d} \beta(t, x, u)\mu(du) \quad \text{and} \quad \sigma(t, x, \mu) = \int_{\mathbb{R}^d} a(t, x, u)\mu(du). $$

Consider a quantizer sequence $x^{(m)} = (x_1^{(m)}, \dots, x_K^{(m)}) \in (\mathbb{R}^d)^K$, $0 \le m \le M$ such that for every $m = \{0, \dots, M\}$, we have $x_i^{(m)} \ne x_j^{(m)}$ if $i \ne j$. Write $(V_k(x^{(m)}))_{1 \le k \le K}$ for a Voronoï partition generated by $x^{(m)}$. Then the transition step of theoretical quantization-based scheme (1.16) can be written as

$$ \tilde{X}_{t_{m+1}} = \hat{X}_{t_m} + h \sum_{k=1}^{K} p_k^{(m)} \beta(t_m, \hat{X}_{t_m}, x_k^{(m)}) + \sqrt{h} \left[ \sum_{k=1}^{K} p_k^{(m)} a(t_m, \hat{X}_{t_m}, x_k^{(m)}) \right] Z_{m+1}, $$

where $p_k^{(m)} = \mathbb{P}(\hat{X}_{t_m} = x_k^{(m)}) = \mathbb{P}(\tilde{X}_{t_m} \in V_k(x^{(m)}))$ is the weight of $x_k^{(m)}$. Hence, given $\hat{X}_{t_m}$ and $p^{(m)}$, $\tilde{X}_{t_{m+1}}$ has a normal distribution

$$ \tilde{X}_{t_{m+1}} \sim \mathcal{N}(\hat{X}_{t_m} + h \sum_{k=1}^{K} p_k^{(m)} \beta(t_m, \hat{X}_{t_m}, x_k^{(m)}), h \left[ \sum_{k=1}^{K} p_k^{(m)} a(t_m, \hat{X}_{t_m}, x_k^{(m)}) \right]^\top \left[ \sum_{k=1}^{K} p_k^{(m)} a(t_m, \hat{X}_{t_m}, x_k^{(m)}) \right] ) $$
since $Z_{m+1} \sim \mathcal{N}(0, \mathbf{I}_q)$. It follows that

$$
\begin{align}
\mathbb{P}(\hat{X}_{t_{m+1}} = x_j^{(m+1)} | \hat{X}_{t_m} = x_i^{(m)}, p^{(m)}) \\
&= \mathbb{P}(\tilde{X}_{t_{m+1}} \in V_j(x^{(m+1)}) | \hat{X}_{t_m} = x_i^{(m)}, p^{(m)}) \\
&= \mathbb{P}\left[(x_i^{(m)} + h \sum_{k=1}^K p_k^{(m)} \beta(t_m, x_i^{(m)}, x_k^{(m)}) + \sqrt{h} \sum_{k=1}^K p_k^{(m)} a(t_m, x_i^{(m)}, x_k^{(m)}) Z_{m+1}) \in V_j(x^{(m+1)})\right] \tag{3.14}
\end{align}
$$

and obviously, by letting

$$
\mathcal{E}_i(x^{(m)}, p^{(m)}, Z_{m+1}) := x_i^{(m)} + h \sum_{k=1}^{K} p_k^{(m)} \beta(t_m, x_i^{(m)}, x_k^{(m)}) + \sqrt{h} \sum_{k=1}^{K} p_k^{(m)} a(t_m, x_i^{(m)}, x_k^{(m)}) Z_{m+1},
$$

we get

$$
\begin{align}
\mathbb{P}(\hat{X}_{t_{m+1}} = x_j^{(m+1)} | p^{(m)}) &= \mathbb{P}(\tilde{X}_{t_{m+1}} \in V_j(x^{(m+1)}) | p^{(m)}) \nonumber \\
&= \sum_{i=1}^{K} \mathbb{P}(\hat{X}_{t_{m+1}} = x_j^{(m+1)} | \hat{X}_{t_m} = x_i^{(m)}, p^{(m)}) \cdot \mathbb{P}(\hat{X}_{t_m} = x_i^{(m)}) \nonumber \\
&= \sum_{i=1}^{K} \mathbb{P}(\mathcal{E}_i(x^{(m)}, p^{(m)}, Z_{m+1}) \in V_j(x^{(m+1)})) \cdot p_i^{(m)} =: g(x^{(m)}, p^{(m)}, Z_{m+1}) \tag{3.15}
\end{align}
$$

Thus one can compute the weights $p_j^{(m+1)} = \mathbb{P}(\hat{X}_{t_{m+1}} = x_j^{(m+1)})$, $1 \le j \le K$, starting from the weight vector $p^{(m)}$ at time $t_m$ by applying (3.15)

The formula (3.15) is valid for any quantizer sequence $x^{(m)}$ with distinct components. Nonetheless,
the quantization errors $\Xi_{t_m} = \|\hat{X}_{t_m} - \hat{X}_{t_m}\|_2$ of each time step $t_m$ impact on the simulation result,
which is also indicated in the error analysis (1.20). One way to reduce the quantization error $\Xi_{t_m}$ is
to integrate Lloyd's algorithm (see Algorithm 1) into the recursive quantization scheme. Given $x^{(m)}$
and $p^{(m)}$, the Lloyd's iteration (3.9) sending $x^{(m+1, [l])} = (x_1^{(m+1, [l])}, \dots, x_K^{(m+1, [l])})$ to $x^{(m+1, [l+1])} =
(x_1^{(m+1, [l+1])}, \dots, x_K^{(m+1, [l+1])})$ for time step $t_{m+1}$ is

$$
x_k^{(m+1, [l+1])} = \frac{\int_{V_k(x^{(m+1, [l])})} \xi \tilde{\mu}_{t_{m+1}}(d\xi)}{\tilde{\mu}_{t_{m+1}}(V_k(x^{(m+1, [l])}))}, \quad k = 1, \dots, K. \quad (3.16)
$$

The denominator of (3.16) can be directly computed by (3.15) while the numerator of (3.16) is essentially
to compute the following value

$$
\begin{align}
& \mathbb{E}[\tilde{X}_{t_{m+1}} \mathbb{1}_{V_j(x^{(m+1)})}(\tilde{X}_{t_{m+1}}) | p^{(m)}] \nonumber \\
&= \mathbb{E}[\mathbb{E}[\tilde{X}_{t_{m+1}} \mathbb{1}_{V_j(x^{(m+1)})}(\tilde{X}_{t_{m+1}}) | \hat{X}_{t_m} = x_i^{(m)}, p^{(m)}] | p^{(m)}] \nonumber \\
&= \sum_{i=1}^{K} \mathbb{E}[\tilde{X}_{t_{m+1}} \mathbb{1}_{V_j(x^{(m+1)})}(\tilde{X}_{t_{m+1}}) | \hat{X}_{t_m} = x_i^{(m)}, p^{(m)}] \cdot \mathbb{P}(\hat{X}_{t_m} = x_i^{(m)}) \nonumber \\
&= \sum_{i=1}^{K} \mathbb{E}[\tilde{X}_{t_{m+1}} \mathbb{1}_{V_j(x^{(m+1)})}(\tilde{X}_{t_{m+1}}) | \hat{X}_{t_m} = x_i^{(m)}, p^{(m)}] \cdot p_i^{(m)}, \tag{3.17}
\end{align}
$$

where

$$
\mathbb{E}[\tilde{X}_{t_{m+1}} \mathbb{1}_{V_j(x^{(m+1)})(\tilde{X}_{t_{m+1}})} | \hat{X}_{t_m} = x_i^{(m)}, p^{(m)}] = \mathbb{E}[Y \mathbb{1}_{V_j(x^{(m+1)})(Y)}] \quad (3.18)
$$

with

$$
Y \sim \mathcal{N}(x_i^{(m)} + h \sum_{k=1}^{K} p_k^{(m)} \beta(t_m, x_i^{(m)}, x_k^{(m)}), h \left[ \sum_{k=1}^{K} p_k^{(m)} a(t_m, x_i^{(m)}, x_k^{(m)}) \right]^\top \left[ \sum_{k=1}^{K} p_k^{(m)} a(t_m, x_i^{(m)}, x_k^{(m)}) \right]). \quad (3.19)
$$
Remark 3.2. In dimension 1, Lloyd’s iteration above is numerically low cost, since the Voronoï cells in dimension 1 are in fact intervals of $\mathbb{R}$. For example, let $x = (x_1, ..., x_K) \in \mathbb{R}^K$ be a quantizer such that $x_i < x_{i+1}$, $i = 1, ..., K-1$, one can choose a Voronoï partition as follows:

$$ 
\begin{aligned}
V_1(x) &= \left(-\infty, \frac{x_1 + x_2}{2}\right) \quad \text{and} \quad V_K(x) = \left[\frac{x_{K-1} + x_K}{2}, +\infty\right), \\
V_k(x) &= \left[\frac{x_{k-1} + x_k}{2}, \frac{x_k + x_{k+1}}{2}\right), \quad k = 2, ..., K-1.
\end{aligned}
$$

Let $x^{(m)} = (x_1^{(m)}, ..., x_K^{(m)})$ be the quantizer of the $m$-th time step. The transition probability $\pi_{ij}^{(m)} = P(\hat{X}_{t_{m+1}} = x_j^{(m+1)} | \hat{X}_{t_m} = x_i^{(m)}, p^{(m)})$ in (3.14) reads

$$ F_{m,\sigma^2}\left(\frac{x_{j+1}^{(m)} + x_j^{(m)}}{2}\right) - F_{m,\sigma^2}\left(\frac{x_{j-1}^{(m)} + x_j^{(m)}}{2}\right), $$

where $F_{m,\sigma^2}$ denotes the cumulative distribution function of $N(m, \sigma^2)$ with

$$ m = x_i^{(m)} + h \sum_{k=1}^{K} p_k^{(m)} \beta(x_i^{(m)}, x_k^{(m)}) \quad \text{and} \quad \sigma = \sqrt{h} \left[ \sum_{k=1}^{K} p_k^{(m)} a(x_i^{(m)}, x_k^{(m)}) \right]. \quad (3.20) $$

Moreover, Lloyd’s iteration (3.18) depends on the value

$$ \int_{(x_{j-1}^{(m)} + x_j^{(m)})/2}^{(x_{j+1}^{(m)} + x_j^{(m)})/2} \xi \cdot f_{m,\sigma^2}(\xi) d\xi \quad (3.21) $$

where $f_{m,\sigma^2}(\xi)$ is the density function of $N(m, \sigma^2)$ with $m$ and $\sigma$ defined in (3.20). In fact, to avoid computing the integral, (3.21) can be alternatively calculated by the following formula: for every $a, b \in \mathbb{R}$,

$$ \int_a^b \xi \cdot f_{m,\sigma^2}(\xi) d\xi = \int_a^b \frac{1}{\sqrt{2\pi}\sigma^2} \xi e^{-\frac{(\xi-m)^2}{2\sigma^2}} d\xi = -\frac{\sigma}{\sqrt{2\pi}} \left[ e^{-\frac{(\xi-m)^2}{2\sigma^2}} \right]_a^b + m[F_{m,\sigma^2}(b) - F_{m,\sigma^2}(a)]. \quad (3.22) $$

However, in higher dimension, the main difficulty of this method is computing the integral of the density function of a multi-dimensional normal distribution over a Voronoï cell (see (3.15), (3.18)). There exists several numerical solutions such as `pyhull` package in `Python` (see also the website `www.qhull.com`) for the cubature formulas of the numerical integration over a convex set in low and medium dimensions ($d = 2, 3, 4$).

## 3.3 $L^2$-error analysis of the hybrid particle-quantization scheme

We prove Proposition 1.4 in this section.

*Proof of Proposition 1.4.* To simplify the notation, we will denote by

$$ b_m^{Q,n} := b(t_m, \tilde{X}_{t_m}^{n,N}, \hat{\mu}_{t_m}^K), \quad b_m^{P,n} := b(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N), \quad \sigma_m^{Q,n} := \sigma(t_m, \tilde{X}_{t_m}^{n,N}, \hat{\mu}_{t_m}^K), \quad \sigma_m^{P,n} := \sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N), $$

where the superscript “P” indicates the *Particle method* and the superscript “Q” indicates the *hybrid particle-Quantization scheme*. Moreover, let $\mathcal{F}_m$ be the $\sigma$-algebra generated by $X_0, \bar{X}_0^{1,N}, ..., \bar{X}_0^{N,N}$, $Z_j^n, n = 1, ..., N, j = 1, ..., m$. Then $b_m^{P,n}, b_m^{Q,n}, \sigma_m^{P,n}, \sigma_m^{Q,n}$ are $\mathcal{F}_m$-measurable and $Z_{m+1}^n, 1 \le n \le N$ are independent of $\mathcal{F}_m$. It follows that

$$ \tilde{X}_{t_{m+1}}^{n,N} - \bar{X}_{t_{m+1}}^{n,N} = \tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N} + h(b_m^{Q,n} - b_m^{P,n}) + \sqrt{h}(\sigma_m^{Q,n} - \sigma_m^{P,n})Z_{m+1}^n, \quad (3.23) $$

then

$$ \begin{aligned} \mathbb{E} \left[ \left| \tilde{X}_{t_{m+1}}^{n,N} - \bar{X}_{t_{m+1}}^{n,N} \right|^2 \right] = \mathbb{E} \left[ \left| (\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}) + h(b_m^{Q,n} - b_m^{P,n}) \right|^2 \right] + \mathbb{E} \left[ \left| \sqrt{h}(\sigma_m^{Q,n} - \sigma_m^{P,n}) Z_{m+1}^n \right|^2 \right] \right. \\ \left. + 2 \mathbb{E} \left[ \left| (\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}) + h(b_m^{Q,n} - b_m^{P,n}) \right| \sqrt{h}(\sigma_m^{Q,n} - \sigma_m^{P,n}) Z_{m+1}^n \right] \right]. \quad (3.24) \end{aligned} $$
The third term in (3.24) equals to 0 since $Z_{m+1}^n, n \in \{1, ..., N\}$, are independent of $\mathcal{F}_m$. Moreover, Assumption I implies,

$$ \mathbb{E}[|b_m^{Q,n} - b_m^{P,n}|^2] \le 2L^2 \{\mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2] + \mathbb{E}[W_2^2(\hat{\mu}_{t_m}^K, \bar{\mu}_{t_m}^N)]\}, $$

and

$$ \begin{aligned} \mathbb{E}[|\sqrt{h}(\sigma_m^{Q,n} - \sigma_m^{P,n})Z_{m+1}^n|^2] &\le h\mathbb{E}[\mathbb{E}[\|\sigma_m^{Q,n} - \sigma_m^{P,n}\|^2 | Z_{m+1}^n|^2 | \mathcal{F}_m]] \\ &= h\mathbb{E}[\|\sigma_m^{Q,n} - \sigma_m^{P,n}\|^2 \mathbb{E}[|Z_{m+1}^n|^2]] = hq\mathbb{E}[\|\sigma_m^{Q,n} - \sigma_m^{P,n}\|^2] \\ &\le 2L^2hq[\mathbb{E}|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2 + \mathbb{E}W_2^2(\hat{\mu}_{t_m}^K, \bar{\mu}_{t_m}^N)]. \end{aligned} $$

Hence, (3.24) becomes

$$ \begin{aligned} & \mathbb{E}[|\tilde{X}_{t_{m+1}}^{n,N} - \bar{X}_{t_{m+1}}^{n,N}|^2] \\ &= \mathbb{E}[| (\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}) + h(b_m^{Q,n} - b_m^{P,n})|^2] + \mathbb{E}[|\sqrt{h}(\sigma_m^{Q,n} - \sigma_m^{P,n})Z_{m+1}^n|^2] \\ &= \mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2] + \mathbb{E}[h^2|b_m^{Q,n} - b_m^{P,n}|^2] + 2h\mathbb{E}[\langle \tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N} | b_m^{Q,n} - b_m^{P,n} \rangle] \\ &\quad + \mathbb{E}[|\sqrt{h}(\sigma_m^{Q,n} - \sigma_m^{P,n})Z_{m+1}^n|^2] \\ &\le \mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2] + 2L^2(h^2 + hq)[\mathbb{E}|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2 + \mathbb{E}W_2^2(\hat{\mu}_{t_m}^K, \bar{\mu}_{t_m}^N)] \\ &\quad + h\mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2 + |b_m^{Q,n} - b_m^{P,n}|^2] \\ &\le \mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2] + 2L^2(h^2 + hq)[\mathbb{E}|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2 + \mathbb{E}W_2^2(\hat{\mu}_{t_m}^K, \bar{\mu}_{t_m}^N)] \\ &\quad + h\mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2] + 2L^2h[\mathbb{E}|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2 + \mathbb{E}W_2^2(\hat{\mu}_{t_m}^K, \bar{\mu}_{t_m}^N)], \\ &\le (1 + 2L^2(h^2 + hq) + h + 2L^2h)\mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2] + (2L^2(h^2 + hq) + 2L^2h)\mathbb{E}W_2^2(\hat{\mu}_{t_m}^K, \bar{\mu}_{t_m}^N). \end{aligned} $$

Remark that for any $m \in \{1, ..., M\}$, the measure $\frac{1}{N} \sum_{n=1}^N \delta_{(\tilde{X}_{t_m}^{n,N}, \bar{X}_{t_m}^{n,N})}$ is a random coupling of $\frac{1}{N} \sum_{n=1}^N \delta_{\tilde{X}_{t_m}^{n,N}}$ and $\bar{\mu}_{t_m}^N = \frac{1}{N} \sum_{n=1}^N \delta_{\tilde{X}_{t_m}^{n,N}}$. Thus, for any $m \in \{1, ..., M\}$,

$$ \begin{aligned} \mathbb{E}[W_2^2(\frac{1}{N} \sum_{n=1}^N \delta_{\tilde{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N})] &\le \mathbb{E}[\int_{\mathbb{R}^d \times \mathbb{R}^d} |x-y|^2 \frac{1}{N} \sum_{n=1}^N \delta_{(\tilde{X}_{t_m}^{n,N}, \bar{X}_{t_m}^{n,N})}(dx, dy)] \\ &= \mathbb{E}[\frac{1}{N} \sum_{n=1}^N |\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2]. \end{aligned} \quad (3.25) $$

Hence, by letting $C_1 := 1 + 2L^2(h^2 + hq) + h + 2L^2h$ and $C_2 := 2L^2(h^2 + hq) + 2L^2h$ and by denoting the quantization error of the time $t_m$ by $\hat{\Xi}_m = W_2(\hat{\mu}_{t_m}^K, \frac{1}{N} \sum_{n=1}^N \delta_{\tilde{X}_{t_m}^{n,N}})$, we get

$$ \begin{aligned} & \sqrt{\frac{1}{N} \sum_{n=1}^{N} \mathbb{E}[|\tilde{X}_{t_{m+1}}^{n,N} - \bar{X}_{t_{m+1}}^{n,N}|^2]} \\ &= \sqrt{\frac{C_1}{N} \sum_{n=1}^{N} \mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2] + C_2 \mathbb{E}[W_2^2(\hat{\mu}_{t_m}^K, \bar{\mu}_{t_m}^N)]} \\ &\le \sqrt{\frac{C_1}{N} \sum_{n=1}^{N} \mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2] + C_2 \{2\mathbb{E}[W_2^2(\hat{\mu}_{t_m}^K, \frac{1}{N} \sum_{n=1}^N \delta_{\tilde{X}_{t_m}^{n,N}})] + 2\mathbb{E}[W_2^2(\frac{1}{N} \sum_{n=1}^N \delta_{\tilde{X}_{t_m}^{n,N}}, \bar{\mu}_{t_m}^N)]\}} \end{aligned} $$
$$ \begin{aligned} & \le \sqrt{(C_1 + 2C_2) \cdot \frac{1}{N} \sum_{n=1}^{N} \mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2] + 2C_2 \mathbb{E} \hat{\Xi}_m^2} \quad (\text{by } (3.25)) \\ & \le \sqrt{C_1 + 2C_2} \sqrt{\frac{1}{N} \sum_{n=1}^{N} \mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2] + \sqrt{2C_2} \sqrt{\mathbb{E} \hat{\Xi}_m^2}} \end{aligned} \quad (3.26) $$

Let $\bar{C}_1 := \sqrt{C_1 + 2C_2}$ and $\bar{C}_2 = \sqrt{2C_2}$. The inequality (3.26) implies

$$ \sqrt{\frac{1}{N} \sum_{n=1}^{N} \mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2]} \le \bar{C}_2 \sum_{j=0}^{m-1} \bar{C}_1^j \sqrt{\mathbb{E} \hat{\Xi}_{m-1-j}^2} $$

Hence, it follows from (3.25) that

$$ \begin{aligned} \mathbb{E} \left[ W_2 \left( \frac{1}{N} \sum_{n=1}^{N} \delta_{\tilde{X}_{t_m}^{n,N}}, \tilde{\mu}_{t_m}^{N} \right) \right] & \le \sqrt{\mathbb{E} W_2^2 \left( \frac{1}{N} \sum_{n=1}^{N} \delta_{\tilde{X}_{t_m}^{n,N}}, \tilde{\mu}_{t_m}^{N} \right)} \\ & \le \sqrt{\frac{1}{N} \sum_{n=1}^{N} \mathbb{E}[|\tilde{X}_{t_m}^{n,N} - \bar{X}_{t_m}^{n,N}|^2]} \\ & \le \bar{C}_2 \sum_{j=0}^{m-1} \bar{C}_1^j \sqrt{\mathbb{E} \hat{\Xi}_{m-1-j}^2} \end{aligned} $$

Consequently,

$$ \begin{aligned} \mathbb{E} \left[ W_2(\tilde{\mu}_{t_m}^K, \tilde{\mu}_{t_m}^N) \right] & \le \mathbb{E} \left[ W_2 \left( \frac{1}{N} \sum_{n=1}^{N} \delta_{\tilde{X}_{t_m}^{n,N}}, \tilde{\mu}_{t_m}^N \right) \right] + \mathbb{E} \left[ W_2 \left( \frac{1}{N} \sum_{n=1}^{N} \delta_{\tilde{X}_{t_m}^{n,N}}, \hat{\mu}_{t_m}^K \right) \right] \\ & \le \bar{C}_2 \sum_{j=0}^{m-1} \bar{C}_1^j \sqrt{\mathbb{E} \hat{\Xi}_{m-1-j}^2} + \mathbb{E} \hat{\Xi}_m. \end{aligned} \quad \square $$

Remark 3.3. For every $m = 1, \dots, M$, it follows from (3.3) that

$$ \mathbb{E} \hat{\Xi}_m = \mathbb{E} \left[ W_2 \left( \hat{\mu}_{t_m}^K, \frac{1}{N} \sum_{n=1}^{N} \delta_{\tilde{X}_{t_m}^{n,N}} \right) \right] = \mathbb{E} \left[ e^{K \left( \frac{1}{N} \sum_{n=1}^{N} \delta_{\tilde{X}_{t_m}^{n,N}} (x^{(m)}) \right)} \right]. $$

Thus one can implement Lloyd's algorithm at each time step in order to reduce the error bound on the right-hand side of (1.23). Moreover, for a fixed $\omega \in \Omega$, finding optimal quantizer for the empirical measure $\frac{1}{N} \sum_{n=1}^{N} \delta_{\tilde{X}_{t_m}^{n,N}}(\omega)$ is equivalent to compute the K-means cluster centers of the sampling $(\tilde{X}_{t_m}^{1,N}(\omega), \dots, \tilde{X}_{t_m}^{N,N}(\omega))$ for which we can use e.g. `sklearn.cluster.KMeans` package in Python.

# 4 Simulation examples

In this section, we illustrate two simulation examples. The first one is the Burgers equation introduced and already considered for numerical tests in [BT97] and [GPPP05]. This is a one-dimensional example with an explicit solution so we use this example to compare the accuracy and computational time of the different numerical methods under consideration. The second example, the Network of FitzHugh-Nagumo neurons, already numerically investigated in [BFFT12], [BFT15] (see also [dRES22]), is a 3-dimensional example. All examples are implemented in Python 3.7.

## 4.1 Simulation of the Burgers equation, comparison of three algorithms

In [BT97], the authors analyse the solution and study the particle method of the Burgers equation

$$ \begin{cases} dX_t = \int_{\mathbb{R}} H(X_t - y)\mu_t(dy)dt + \sigma dB_t, \\ \forall t \in [0, T], \mu_t = P_{X_t} \text{ and } X_0 := (\Omega, \mathcal{F}, \mathbb{P}) \to (\mathbb{R}, \mathcal{B}(\mathbb{R})), \end{cases} \quad (4.1) $$
where $H$ is the Heaviside function $H(z) = 1_{\{z \ge 0\}}$ and $\sigma$ is a real constant. The cumulative distribution function $V(t, x)$ of $\mu_t$ satisfies

$$ \frac{\partial V}{\partial t} = \frac{1}{2}\sigma^2\frac{\partial^2 V}{\partial x^2} - V\frac{\partial V}{\partial x}, \quad V(0, x) = V_0(x). $$

Moreover, if the initial cumulative distribution function $V_0$ satisfies $\int_0^x V_0(y)dy = \mathcal{O}(x)$, then the function $V$ has a closed form given by (see [Hop50])

$$ V(t, x) = \frac{\int_{\mathbb{R}} V_0(y)\exp\left(-\frac{1}{\sigma^2}\left[\frac{(x-y)^2}{2t} + \int_0^y V_0(z)dz\right]\right)dy}{\int_{\mathbb{R}} \exp\left(-\frac{1}{\sigma^2}\left[\frac{(x-y)^2}{2t} + \int_0^y V_0(z)dz\right]\right)dy}, \quad (t, x) \in [0, T] \times \mathbb{R}. $$

Hence, if we consider $X_0 = 0$, then the cumulative distribution function at time $T = 1$ reads

$$ F_{T=1}(x) = \frac{\int_{\mathbb{R}_+} \exp\left(-\frac{1}{\sigma^2}\left[\frac{(x-y)^2}{2} + y\right]\right)dy}{\int_{\mathbb{R}} \exp\left(-\frac{1}{\sigma^2}\left[\frac{(x-y)^2}{2} + y1_{\{y \ge 0\}}\right]\right)dy}. \quad (4.2) $$

So we can compare the accuracy of the different numerical methods proposed in the former sections by computing

$$ \|F_{\text{simu}} - F_{\text{true}}\|_{\text{sup}}, \quad (4.3) $$

where $F_{\text{simu}}$ represents the simulated cumulative distribution function by different numerical methods and $F_{\text{true}}$ is the true cumulative distribution function (4.2). For two probability distributions $\mu, \nu \in \mathcal{P}_p(\mathbb{R}^d)$ with respective cumulative distribution function $F$ and $G$, the Wasserstein distance $W_p(\mu, \nu)$ can be computed by

$$ W_p^p(\mu, \nu) = \int_0^1 |F^{-1}(u) - G^{-1}(u)|^p du, \quad p \ge 1. \quad (4.4) $$

However, it is computationally extremely costly to directly compute the inverse function of the cumulative distribution function (4.2) and if we compute (4.4) by using Monte-Carlo simulation, it will induce its own statistical error which may perturb our comparisons. Thus, instead of considering (4.4), we choose to directly compute (4.3) by

$$ \|F_{\text{simu}}(x) - F_{\text{true}}(x)\|_{\text{sup}} \approx \sup_{x \in \text{Unifset}} |F_{\text{simu}}(x) - F_{\text{true}}(x)|, $$

where Unifset is a uniformly spaced point set in $[-2.5, 3.5]$.

In the following simulation, we choose $\sigma^2 = 0.2$ and $M = 50$ so that we have the same time step $h = T/M = 0.02$ for each method. We first give preliminary illustrations of the simulated cumulative distribution function by the particle method (1.8), the recursive quantization scheme (1.16)-(1.18)-(1.19) (with and without Lloyd's algorithm) and the hybrid particle-quantization scheme (1.22) in Figure 3, 4, 5, 6 and 7. The true cumulative distribution functions are marked in red while the simulated ones are marked in blue.

<figure><img src="image_6.png" /><figcaption>Figure 3: Simulation by the particle method (1.8)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: Simulation by the recursive quantization scheme without Lloyd's iteration (1.16)-(1.19)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5: Simulation by the recursive quantization scheme with 5 Lloyd's iterations at each Euler step (1.16)-(1.19) and (3.15)-(3.18)</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6: Simulation by the hybrid particle-quantization scheme without Lloyd iteration (1.22)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Simulation by the hybrid particle-quantization scheme with 5 Lloyd iterations at each Euler step (1.22)</figcaption></figure>

A detailed comparison of the different methods is displayed in Table 1. Remark that the particle method (1.8) and the hybrid particle-quantization scheme (1.22) are random algorithms. Hence, their accuracy are computed by taking an average error computed over 50 independent simulations.

<table><caption>Table 1: Comparison of three algorithms</caption><thead><tr><th></th><th>Particle size and quantizer size</th><th>Computing time for each Euler step</th><th>Error ||F<sub>simu</sub> - F<sub>true</sub>||<sub>sup</sub></th></tr></thead><tbody><tr><td>Particle method (1.8)</td><td>Particle size N = 10000</td><td>0.00320s</td><td>0.01021</td></tr><tr><td>Recursive quantization scheme without Lloyd iterations (1.16)-(1.19)</td><td>Quantizer size K = 500</td><td>0.00205s</td><td>0.01054</td></tr><tr><td>Recursive quantization scheme with 5 Lloyd iterations at each Euler step (1.16)-(1.19) and (3.15)-(3.18)</td><td>Quantizer size K = 500</td><td>8.21598s</td><td>0.01029</td></tr><tr><td>Hybrid particle-quantization scheme (1.22) without Lloyd iterations</td><td>Particle size N = 10000<br>Quantizer size K = 500</td><td>6.09480s</td><td>0.01626</td></tr><tr><td>Hybrid particle-quantization scheme (1.22) with 5 Lloyd iterations at each Euler step</td><td>Particle size N = 10000<br>Quantizer size K = 500</td><td>9.37229s</td><td>0.01013</td></tr></tbody></table>

In a second phase, we show the converging rate of the simulation error (4.3) of the particle method (1.8) and the recursive quantization scheme without Lloyd quantizer optimization (1.16)-(1.18)-(1.19) respectively according to N and K. Table 2 shows the simulation error (4.3) of the particle method (1.8) with respect to the numbers of particle N = 2<sup>8</sup>, ..., 2<sup>13</sup> for a fixed M = 50. As the particle method (1.8) is a random algorithm, the error shows in Table 2 is computed by reruning independently 500 times for each value of N.

<table><caption>Table 2: Error of the particle method (1.8) with respect to the particle size N</caption><thead><tr><th>N</th><th>2<sup>8</sup></th><th>2<sup>9</sup></th><th>2<sup>10</sup></th><th>2<sup>11</sup></th><th>2<sup>12</sup></th><th>2<sup>13</sup></th></tr></thead><tbody><tr><td>Error ||F<sub>simu</sub> - F<sub>true</sub>||<sub>sup</sub></td><td>0.04691</td><td>0.03409</td><td>0.02438</td><td>0.01785</td><td>0.01407</td><td>0.01131</td></tr><tr><td>Standard deviation</td><td>0.01207</td><td>0.00939</td><td>0.00687</td><td>0.00469</td><td>0.00408</td><td>0.00294</td></tr></tbody></table>

Figure 8, 9 and 10 show the curve of the particle method errors with respect to N, the respective log-log
error and the standard deviation of the error.

<figure><img src="image_11.png" /><figcaption>Figure 8: Error of the particle method (1.8) with respect to the particle size $N$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: Log-error of the particle method (1.8) with respect to $\log_2(N)$. The slope is approximately equal to $-0.285$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10: The standard deviation of the error of the particle method (1.8) with respect to $N$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11: Error of the recursive quantization scheme (1.16)-(1.18)-(1.19) with respect to the quantizer size $K$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 12: Log-error of the recursive quantization scheme (1.16)-(1.18)-(1.19) with respect to $\log_2(K)$. The slope is approximately equal to $-0.436$.</figcaption></figure>

Table 3 shows the convergence rate of the error of the recursive quantization scheme (1.16)-(1.18) with respect to the quantizer size $K = 2^5, 2^6, 2^7, 2^8, 2^9, 2^{10}$. Here we use a fixed quantizer sequence which is a uniformly spaced point set in [-2.5, 3.5] without Lloyd's algorithm for each time step $t_m$.

<table><caption>Table 3: Error of the recursive quantization scheme (1.16)-(1.18) with respect to the quantizer size $K$. Figure 11 and 12 above show the curve of the error with respect to $K$ and the respective log-log error.</caption><thead><tr><th>K</th><th>2<sup>5</sup></th><th>2<sup>6</sup></th><th>2<sup>7</sup></th><th>2<sup>8</sup></th><th>2<sup>9</sup></th><th>2<sup>10</sup></th></tr></thead><tbody><tr><td>Error ||F<sub>simu</sub> - F<sub>true</sub>||<sub>sup</sub></td><td>0.07347</td><td>0.04176</td><td>0.02360</td><td>0.01471</td><td>0.01043</td><td>0.00829</td></tr></tbody></table>

Now we give some comments on the numerical performance of different methods mainly through two aspects: the accuracy and the computing time.

### Comparison of the computing time.

The particle method (1.8) and the recursive quantization scheme (1.16)-(1.18)-(1.19) without Lloyd iteration are the two fastest methods. In fact, these two methods are essentially computing a Markov chain in $\mathbb{R}^N$ and $\mathbb{R}^K$ respectively. The application of Lloyd's iteration in the recursive quantization scheme is slightly faster than in the hybrid particle-quantization scheme since we used the closed formulas (3.22). We can also remark that, when the quantization level $K$ is large, Lloyd's iteration is numerically costly even in dimension 1 and the application of Lloyd's algorithm only slightly reduces the quantization error.
* Comparison of the accuracy computed by $||F_{simu}(x) - F_{true}(x)||_{sup}$.

- The particle method (1.8) and the hybrid particle-quantization scheme (1.22) are random algorithms, whose simulation results, including the error $||F_{simu}(x) - F_{true}(x)||_{sup}$, depend on $\omega$ in $(\Omega, \mathcal{F}, \mathbb{P})$. Therefore, when considering the accuracy of these two algorithms, one need to compute in the same time the standard deviation of errors. However, the recursive quantization scheme (1.16)-(1.18)-(1.19) is deterministic.

- Comparing with the particle size $N$ in Table 2 and the quantizer size $K$ in Table 3, one can remark that to achieve the same accuracy, we need fewer points in the quantizer than in the particle. So if we need a discrete representation of the cumulative distribution function $F$ (or equivalently, a discrete representation of the probability distribution $\mu$) to compute a further functional of $\mu$, such as an integral with respect to $\mu$, the recursive quantization based scheme provides a smaller data set ($K$-size quantizer and $K$-size weight vector) than the particle method.

- The error of the recursive quantization scheme (1.16)-(1.18)-(1.19), especially when we implement without Lloyd's quantizer optimization, strongly depends on the choice of quantizer. Generally, a practical way to choose the initial quantizer of a probability distribution $\mu$ is to use self-quantization technique for which we refer to [DGLP06], [GL00][Section 7.1 and Section 14], [PP03] and [PPP04]. Another efficient trick is to rely on a so-called “splitting method”, that is, using the trained quantizer of the Euler step $t_m$ as an initial quantizer of the Euler step $t_{m+1}$.

In this one dimensional example, we did not remark any obvious advantage of the hybrid particle quantization scheme (1.22) comparing with other methods. However, in the next section, we will show that the hybrid method provides a fair balance between the accuracy and the data volume.

## 4.2 Simulation of the network of FitzHugh-Nagumo neurons in dimension 3

In this section, we consider the network of FitzHugh-Nagumo neurons introduced in [BFFT12] and [BFT15]:

$$dX_t = b(X_t, \mu_t)dt + \sigma(X_t, \mu_t)dB_t \quad (4.5)$$

where $b: \mathbb{R}^3 \times \mathcal{P}(\mathbb{R}^3) \to \mathbb{R}^3$ and $\sigma: \mathbb{R}^3 \times \mathcal{P}(\mathbb{R}^3) \to \mathbb{M}_{3 \times 3}$ are defined by

$$b(x, \mu) := \begin{pmatrix} x_1 - (x_1)^3/3 - x_2 + I - \int_{\mathbb{R}^3} J(x_1 - V_{rev})z_3 \mu(dz) \\ \frac{c(x_1 + a - bx_2)}{a_r \frac{T_{\max}(1-x_3)}{1+\exp(-\lambda(x_1-V_T))} - a_d x_3} \end{pmatrix}, \quad (4.6)$$

$$\sigma(x, \mu) := \begin{pmatrix} \sigma_{ext} & 0 & -\int_{\mathbb{R}^3} \sigma_J(x_1 - V_{rev})z_3 \mu(dz) \\ 0 & 0 & 0 \\ 0 & \sigma_{32}(x) & 0 \end{pmatrix},$$

with

$$\sigma_{32} := \mathbb{1}_{x_3 \in (0,1)} \sqrt{a_r \frac{T_{\max}(1-x_3)}{1 + \exp(-\lambda(x_1 - V_T))} + a_d x_3 \Gamma \exp(-\frac{\Lambda}{1 - (2x_3 - 1)^2})},$$

and where

$$X_0 \sim \mathcal{N} \left( \begin{pmatrix} V_0 \\ \omega_0 \\ y_0 \end{pmatrix}, \begin{pmatrix} \sigma_{V_0} & 0 & 0 \\ 0 & \sigma_{\omega_0} & 0 \\ 0 & 0 & \sigma_{y_0} \end{pmatrix} \right)$$

with the following parameter values

<table><tr><td>V<sub>0</sub> = 0</td><td>σ<sub>V<sub>0</sub></sub> = 0.4</td><td>a = 0.7</td><td>b = 0.8</td><td>c = 0.08</td><td>I = 0.5</td><td>σ<sub>ext</sub> = 0.5</td></tr><tr><td>ω<sub>0</sub> = 0.5</td><td>σ<sub>ω<sub>0</sub></sub> = 0.4</td><td>V<sub>rev</sub> = 1</td><td>a<sub>r</sub> = 1</td><td>a<sub>d</sub> = 1</td><td>T<sub>max</sub> = 1</td><td>λ = 0.2</td></tr><tr><td>y<sub>0</sub> = 0.3</td><td>σ<sub>y<sub>0</sub></sub> = 0.05</td><td>J = 1</td><td>σ<sub>J</sub> = 0.2</td><td>V<sub>T</sub> = 2</td><td>Γ = 0.1</td><td>Λ = 0.5.</td></tr></table>
This model was first studies in [BFFT12] and then rigorously investigated and analyzed in [BFT15]. We are aware that the equation (4.5) does not fulfill Assumption I since the drift b defined by (4.6) is only locally Lipschitz in $(x, \mu)$. However, the drift b satisfies

$$ \langle x | b(x, \mu) \rangle \leq c_1 (|x|^2 + 1) $$

for some $c_1 > 0$ and the coefficient $\sigma$ is bounded. This ensures the existence of a strong solution living on the whole $\mathbb{R}_+$. In fact the presence of the term $-\frac{(x_1)^3}{3}$ also induces a mean-reversion property that makes possible to control the long-run behaviour of the equation. On the other hand such drift with non linear growth (in norm) usually induces an instability of the Euler scheme as emphasized e.g. in [Lem07], at least when using equidistant time steps. We nevertheless chose this model for a 3D numerical illustration due to its challenging feature, but with a refined time step ($M = 5000$, see below) to ensure its stability. Other numerical schemes, typically semi-implicit could be more stable but out of the scope of this paper.

In this section, we compare the performance of the particle method (1.8) and the hybrid particle-quantization scheme (1.22) in two aspects. First, we intuitively compare these two methods by simulating the density function of $(x_1, x_2)$ for $T = 1.5$, as in the original paper [BFFT12][Page 31, Figure 4, the third one in the right]. In this step, we choose $M = 5000$ for the temporal discretization to reduce (as much as possible) the error induced by the Euler scheme. The images of the density function simulated respectively by the particle method and the hybrid method are displayed in Figure 13, 14 and Figure 15, 16, 17 respectively.

Next, as the particle method and the hybrid method are both random methods, we take

$$ \varphi(\mu_T^{\text{simu}}) := \int_{\mathbb{R}^3} |\xi|^2 \mu_T^{\text{simu}}(d\xi) = \mathbb{E}[|X_T^{\text{simu}}|^2] $$

as a test function for the simulated distribution $\mu_T^{\text{simu}}$ at time $T$, rerun 200 times for each method and compare the mean and the standard deviation of $\varphi(\mu_T^{\text{simu}})$.

The particle method (1.8):

<figure><img src="image_16.png" /><figcaption>Figure 13: The first and second coordinates of 5000 particles at time $T = 1.5$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 14: The simulated density function smoothened by the Gaussian kernel method (bandwith = 0.241).</figcaption></figure>
The hybrid particle-quantization scheme (1.22):

<figure><img src="image_18.png" /><figcaption>Figure 15: The quantizer of $(x_0, x_1)$, simulated with particle number $N = 5000$, quantizer size $K = 300$ and 10 Lloyd iterations at each Euler step.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 16: The Voronoï cells of the above quantizer. The colour of each Voronoï cell represents the weight of this cell (the darker the heavier).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 17: The density function simulated by the hybrid particle-quantization method (1.22). The vertical axis is the weight divided by the area of the corresponding Voronoï cell.</figcaption></figure>

The obtained density functions have a similar form by these two methods but the data volume obtained by the particle method is

5000 (the number of particle) × 3 (dimension),

while the data volume obtained by the hybrid method is

300 (the quantizer size) × 4 (dimension + weight for each quantizer).

For a more precise comparison, we set now the time discretization number $M = 150$ and we consider the following test function for the simulated distribution $\mu_T^{\text{simu}}$ at $T = 1.5$

$$ \varphi(\mu_T^{\text{simu}}) := \int_{\mathbb{R}^3} |\xi|^2 \mu_T^{\text{simu}}(d\xi) = \mathbb{E}[|X_T^{\text{simu}}|^2] $$

and rerun 200 times for each method. Table 4 shows the simulation results.
<table><caption>Table 4: Comparison of the simulation result $\varphi(\mu_T^{\text{simu}})$.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Particle method</th><th>Hybrid method</th></tr><tr><th>N=300</th><th>N=5000</th><th>N = 5000 and K = 300</th></tr></thead><tbody><tr><td>Data volume for μ<sub>T</sub><sup>simu</sup></td><td>300×3</td><td>5000×3</td><td>300×4</td></tr><tr><td>Average computing time for each Euler step</td><td>0.034s</td><td>0.608s</td><td>2.670s</td></tr><tr><td rowspan="2">Test function ϕ(μ<sub>T</sub><sup>simu</sup>)</td><td>Mean</td><td>1.194</td><td>1.205</td><td>1.192</td></tr><tr><td>Standard deviation</td><td>0.058</td><td>0.015</td><td>0.015</td></tr></tbody></table>

Intuitively, the hybrid method can be considered as adding a “feature extraction” step on the particle method. Comparing the third and fourth columns of the above Table 4, one can notice that this added step needs more computing time but reduces the size of the output data for the further computing of the test function $\varphi(\mu_T^{\text{simu}})$ without increasing the standard deviation. However, the second column of the above table shows that if we implement the particle method with a similar data size than the quantization level, the test value of $\varphi(\mu_T^{\text{simu}})$ provides a much larger standard deviation.

## A Proofs of Lemma 2.3, Lemma 2.4 and Lemma 2.7

*Proof of Lemma 2.3.* We consider the canonical space $\Omega = C([0, T], \mathbb{R}^d) \times C([0, T], \mathbb{R}^d)$ equipped with the $\sigma$-algebra $\mathcal{F}$ generated by the distance

$$d((\omega^1, \omega^2), (\alpha^1, \alpha^2)) := \|\omega^1 - \alpha^1\|_{\sup} \lor \|\omega^2 - \alpha^2\|_{\sup}$$

and $\mathbb{P} \in \bar{\Pi}(\mu, \nu)$ where $\bar{\Pi}(\mu, \nu)$ is the set of probability measures with marginals $\mu$ and $\nu$. For any $\omega = (\omega^1, \omega^2) \in \Omega$, we define the canonical projections $X : \Omega \to C([0, T], \mathbb{R}^d)$ and $Y : \Omega \to C([0, T], \mathbb{R}^d)$ by

$$\forall \omega = (\omega^1, \omega^2), \quad \forall t \in [0, T], \quad X_t(\omega) = \omega_t^1 \text{ and } Y_t(\omega) = \omega_t^2.$$

The couple $(X, Y)$ makes up the canonical process on $\Omega$. Since $\mathbb{P} \in \bar{\Pi}(\mu, \nu)$, then $X$ has probability distribution $\mu$ and $Y$ has probability distribution $\nu$. Moreover, we have

$$\sup_{s \in [0,t]} \mathcal{W}_p^p(\mu_s, \nu_s) \le \sup_{s \in [0,t]} \mathbb{E} |X_s - Y_s|^p \le \mathbb{E} \sup_{s \in [0,t]} |X_s - Y_s|^p.$$

Then we can choose by the usual arguments $\mathbb{P} \in \bar{\Pi}(\mu, \nu)$ such that

$$\mathbb{E} \sup_{s \in [0,t]} |X_s - Y_s|^p = (\mathcal{W}_{p,t}(\mu_s, \nu_s))^p$$

to conclude the proof.

□

*Proof of Lemma 2.4.* (a) For any $x \in \mathbb{R}^d$ and for any $\mu \in \mathcal{P}_p(\mathbb{R}^d)$, Assumption I implies that

$$\forall t \in [0, T], \quad |b(t, x, \mu)| - |b(t, 0, \delta_0)| \le |b(t, x, \mu) - b(t, 0, \delta_0)| \le L(|x| + \mathcal{W}_p(\mu, \delta_0)).$$

Hence,

$$|b(t, x, \mu)| \le |b(t, 0, \delta_0)| + L(|x| + \mathcal{W}_p(\mu, \delta_0)) \le (|b(t, 0, \delta_0)| \lor L)(1 + |x| + \mathcal{W}_p(\mu, \delta_0))$$

Similarly, we have $\|\sigma(t, x, \mu)\| \le (\|\sigma(t, 0, \delta_0)\| \lor L)(1 + |x| + \mathcal{W}_p(\mu, \delta_0))$, so we can take $C_{b,\sigma,L,T} = \sup_{t \in [0,T]} |b(t, 0, \delta_0)| \lor \sup_{t \in [0,T]} \|\sigma(t, 0, \delta_0)\| \lor L$ to complete the proof.
(b) STEP 1. We prove $\sup_{0 \le m \le M} \sup_{1 \le n \le N} \|\bar{X}_{t_m}^{n,N}\|_p < +\infty$ by induction in this step. First, we know that for every $n = 1, \dots, N$, $\|\bar{X}_0^{n,N}\|_p = \|X_0\|_p < +\infty$ by the definition of $\bar{X}_0^{n,N}$ so that $\sup_{1 \le n \le N} \|\bar{X}_0^{n,N}\|_p < +\infty$. Now assume that $\sup_{1 \le n \le N} \|\bar{X}_{t_m}^{n,N}\|_p < +\infty$. For every $n = 1, \dots, N$, we have

$$ \begin{aligned} \|\bar{X}_{t_{m+1}}^{n,N}\|_p &\le \|\bar{X}_{t_m}^{n,N}\|_p + h \cdot \|b(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N)\|_p + \sqrt{h} \cdot \|\sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) Z_{m+1}^n\|_p \\ &\le \|\bar{X}_{t_m}^{n,N}\|_p + h \cdot \|b(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N)\|_p + \sqrt{h} \cdot \left\|\|\sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N)\|\right\|_p \|Z_{m+1}^n\|_p \end{aligned} $$

where the first inequality above follows from the Minkowski inequality and the second inequality follows from the fact that $Z_{m+1}^n$ is independent of $\mathcal{F}_{t_m}$ and $\sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N)$ is $\mathcal{F}_{t_m}$ measurable. By applying the inequality (2.8) in the first part, we obtain

$$ \|b(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N)\|_p \vee \left\| \|\sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N)\| \right\|_p \le C_{b,\sigma,L,T} \left( 1 + \|\bar{X}_{t_m}^{n,N}\|_p + \|\mathcal{W}_p(\bar{\mu}_{t_m}^N, \delta_0)\|_p \right) \quad (A.1) $$

and by the definition of $\bar{\mu}_{t_m}^N$, we have

$$ \| \mathcal{W}_p^p(\bar{\mu}_{t_m}^N, \delta_0) \|_p^p = \mathbb{E} \left[ \mathcal{W}_p^p \left( \frac{1}{N} \sum_{n=1}^N \delta_{\bar{X}_{t_m}^{n,N}}, \delta_0 \right) \right] \le \frac{1}{N} \sum_{n=1}^N \mathbb{E} \left[ |\bar{X}_{t_m}^{n,N}|^p \right] < +\infty. \quad (A.2) $$

Hence, for every $n = 1, \dots, N$, $\|\bar{X}_{t_{m+1}}^{n,N}\|_p < +\infty$. Then we can conclude $\sup_{0 \le m \le M} \sup_{1 \le n \le N} \|\bar{X}_{t_m}^{n,N}\|_p < +\infty$ by induction.

STEP 2. We prove (2.9) in this step. The first step implies that

$$ \left\| \sup_{0 \le m \le M} \sup_{1 \le n \le N} |\bar{X}_{t_m}^{n,N}| \right\|_p^p \le \sum_{m=0}^{M} \sum_{n=1}^{N} \mathbb{E} \left[ |\bar{X}_{t_m}^{n,N}|^p \right] \le (M+1)N \cdot \sup_{0 \le m \le M} \sup_{1 \le n \le N} \|\bar{X}_{t_m}^{n,N}\|_p^p < +\infty. \quad (A.3) $$

Then,

$$ \begin{aligned} & \left\| \sup_{t \in [0,T]} |\bar{X}_t^{n,N}| \right\|_p = \left\| \sup_{0 \le m \le M} \sup_{t \in [t_m, t_{m+1}]} |\bar{X}_t^{n,N}| \right\|_p \\ &\le \left\| \sup_{0 \le m \le M} \sup_{t \in [t_m, t_{m+1}]} |\bar{X}_{t_m}^{n,N} + (t-t_m)b(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) + \sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N)(B_t^n - B_{t_m}^n)| \right\|_p \\ &\le \left\| \sup_{0 \le m \le M} \sup_{t \in [t_m, t_{m+1}]} |\bar{X}_{t_m}^{n,N}| + |b(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N)| \cdot |t-t_m| + \left\| \sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) \right\| \|B_t^n - B_{t_m}^n\| \right\|_p \\ &\le \left\| \sup_{0 \le m \le M} |\bar{X}_{t_m}^{n,N}| + h \left\| b(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) \right\| \right\|_p + \left\| \sup_{0 \le m \le M} \left\| \sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) \right\| \sup_{t \in [t_m, t_{m+1}]} |B_t^n - B_{t_m}^n| \right\|_p \\ &\quad \text{(by applying the Minkowski inequality)} \end{aligned} $$

$$ \begin{aligned} & \le \left\| \sup_{0 \le m \le M} |\bar{X}_{t_m}^{n,N}| + h \left\| b(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) \right\| \right\|_p + \sum_{m=0}^{M} \left\| \left\| \sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) \right\| \sup_{t \in [t_m, t_{m+1}]} |B_t^n - B_{t_m}^n| \right\|_p \\ &\quad \text{(as for every $a, b \in \mathbb{R}_+$, $p \ge 2$, $\{\mathbb{E}[a^p \vee b^p]\}^{1/p} \le \{\mathbb{E}[a^p] + \mathbb{E}[b^p]\}^{1/p} \le \{\mathbb{E}[a^p]\}^{1/p} + \{\mathbb{E}[b^p]\}^{1/p})} \\ & \le \left\| \sup_{0 \le m \le M} |\bar{X}_{t_m}^{n,N}| + h \left\| b(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) \right\| \right\|_p + \sum_{m=0}^{M} \left\| \left\| \sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) \right\| \sup_{t \in [t_m, t_{m+1}]} |B_t^n - B_{t_m}^n| \right\|_p \\ &\quad \text{(as } \left\| \sigma(t_m, \bar{X}_{t_m}^{n,N}, \bar{\mu}_{t_m}^N) \right\| \text{ is } \mathcal{F}_{t_m} \text{ mesurable and } \sup_{t \in [t_m, t_{m+1}]} |B_t^n - B_{t_m}^n| \text{ is independent of } \mathcal{F}_{t_m}) \\ &< +\infty \end{aligned} $$

by applying (2.8), (A.1) (A.2) and (A.3)

□
For the proof of Lemma 2.7, we need the following two inequalities and we refer to [Pag18][Section 7.8] among other references for more details.

**Lemma A.1** (The Generalized Minkowski Inequality). For any (bi-measurable) process $X = (X_t)_{t \ge 0}$, for every $p \in [1, \infty)$ and for every $T \in [0, +\infty]$,

$$ \left\| \int_{0}^{T} X_t dt \right\|_{p} \le \int_{0}^{T} \| X_t \|_{p} dt. $$

**Lemma A.2** (Burkölder-Davis-Gundy Inequality (continuous time)). For every $p \in (0, +\infty)$, there exists two real constants $c_p^{BDG} > 0$ and $C_p^{BDG} > 0$ such that, for every continuous local martingale $(X_t)_{t \in [0,T]}$ null at 0,

$$ c_p^{BDG} \left\| \sqrt{\langle X \rangle_T} \right\|_p \le \left\| \sup_{t \in [0,T]} |X_t| \right\|_p \le C_p^{BDG} \left\| \sqrt{\langle X \rangle_T} \right\|_p. $$

In particular, if $(B_t)$ is an $(\mathcal{F}_t)$-standard Brownian motion and $(H_t)_{t \ge 0}$ is an $(\mathcal{F}_t)$-progressively measurable process having values in $M_{d,q}(\mathbb{R})$ such that $\int_0^T \|H_t\|^2 dt < +\infty \quad \mathbb{P}-\text{a.s.}$, then the $d$-dimensional local martingale $\int_0^T H_s dB_s$ satisfies

$$ \left\| \sup_{t \in [0,T]} \left| \int_0^t H_s dB_s \right| \right\|_p \le C_{d,p}^{BDG} \left\| \sqrt{\int_0^T \|H_t\|^2 dt} \right\|_p. \quad (A.4) $$

*Proof of Lemma 2.7.* For every

$$ (X, (\mu_t)_{t \in [0,T]}), (Y, (\nu_t)_{t \in [0,T]}) \in L_{C([0,T], \mathbb{R}^d)}^p (\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \ge 0}, \mathbb{P}) \times \mathcal{C}([0,T], \mathcal{P}_p(\mathbb{R}^d)) $$

for every $t \in [0, T]$, we have

$$ \begin{aligned} & \left\| \sup_{s \in [0,t]} \left| \int_0^t [b(u, X_u, \mu_u) - b(u, Y_u, \nu_u)] du \right| \right\|_p \\ & \le \left\| \sup_{s \in [0,t]} \int_0^s |b(u, X_u, \mu_u) - b(u, Y_u, \nu_u)| du \right\|_p = \left\| \int_0^t |b(u, X_u, \mu_u) - b(u, Y_u, \nu_u)| du \right\|_p \\ & \le \int_0^t \|b(u, X_u, \mu_u) - b(u, Y_u, \nu_u)\|_p du \quad (\text{by Lemma A.1}) \\ & \le \int_0^t \|L[|X_u - Y_u| + W_p(\mu_u, \nu_u)]\|_p du \le L \int_0^t [\|X_u - Y_u\|_p + \|W_p(\mu_u, \nu_u)\|_p] du, \end{aligned} $$

and

$$ \begin{aligned} & \left\| \sup_{s \in [0,t]} \left| \int_0^s [\sigma(u, X_u, \mu_u) - \sigma(u, Y_u, \nu_u)] dB_u \right| \right\|_p \\ & \le C_{d,p}^{BDG} \left\| \sqrt{\int_0^t \|\sigma(u, X_u, \mu_u) - \sigma(u, Y_u, \nu_u)\|^2 du} \right\|_p \quad (\text{by Lemma A.2}) \\ & \le C_{d,p}^{BDG} \left\| \int_0^t \|\sigma(u, X_u, \mu_u) - \sigma(u, Y_u, \nu_u)\|^2 du \right\|_{\frac{p}{2}}^{\frac{1}{2}} \quad (\text{since } \left\| \sqrt{U} \right\|_p = [\mathbb{E}U^{\frac{p}{2}}]^{\frac{2}{p} \times \frac{1}{2}} = \|U\|_{\frac{p}{2}}^{\frac{1}{2}}, \text{ when } U \ge 0) \\ & \le C_{d,p}^{BDG} \left[ \int_0^t \left\| \|\sigma(u, X_u, \mu_u) - \sigma(u, Y_u, \nu_u)\|^2 \right\|_{\frac{p}{2}} du \right]^{\frac{1}{2}} \quad (\text{by Minkowski's inequality, since } p \in [2, +\infty)) \\ & \le C_{d,p}^{BDG} \left[ \int_0^t \left\| \|\sigma(u, X_u, \mu_u) - \sigma(u, Y_u, \nu_u)\| \right\|_p^2 du \right]^{\frac{1}{2}} \quad (\text{since } \left\| U \right\|_{\frac{p}{2}}^2 = \left[ (\mathbb{E}|U|^p)^{\frac{1}{p}} \right]^2 = \|U\|_p^2) \\ & \le C_{d,p}^{BDG} \left[ \int_0^t \|L[|X_u - Y_u| + W_p(\mu_u, \nu_u)]\|_p^2 du \right]^{\frac{1}{2}} \quad (\text{by Assumption I}) \end{aligned} $$
$$
\begin{aligned}
&\le C_{d,p}^{BDG} L \left[ \int_0^t \left[ \|X_u - Y_u\|_p + \|W_p(\mu_u, \nu_u)\|_p \right]^2 du \right]^{\frac{1}{2}} \\
&\le \sqrt{2} C_{d,p}^{BDG} L \left[ \int_0^t \left[ \|X_u - Y_u\|_p^2 + \|W_p(\mu_u, \nu_u)\|_p^2 \right] du \right]^{\frac{1}{2}}.
\end{aligned}
$$

Then we can conclude the proof by setting $C_{d,p,L} = \sqrt{2} C_{d,p}^{BDG} L$. $\square$

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

B The organisation of different numerical schemes

<figure><img src="image_22.png" /><figcaption>Figure 18: Numerical methods of this paper</figcaption></figure>
References

[AKH02] Fabio Antonelli and Arturo Kohatsu-Higa. Rate of convergence of a particle method to the solution of the McKean-Vlasov equation. *The Annals of Applied Probability*, 12(2):423-476, 2002.

[BFFT12] Javier Baladron, Diego Fasoli, Olivier Faugeras, and Jonathan Touboul. Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons. *Journal of Mathematical Neuroscience*, 2:Art. 10, 50, 2012.

[BFT15] Mireille Bossy, Olivier Faugeras, and Denis Talay. Clarification and complement to “Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and FitzHugh-Nagumo neurons”. *The Journal of Mathematical Neuroscience*, 5(1):1-23, 2015.

[BJR22] Mireille Bossy, Jean-François Jabir, and Kerlyns Martinez Rodriguez. Instantaneous turbulent kinetic energy modelling based on lagrangian stochastic approach in CFD and application to wind energy. *Journal of Computational Physics*, 464:110929, 2022.

[Bol08] François Bolley. Separability and completeness for the Wasserstein distance. In *Séminaire de probabilités XLI*, pages 371-377. Springer, 2008.

[BT97] Mireille Bossy and Denis Talay. A stochastic particle method for the McKean-Vlasov and the Burgers equation. *Mathematics of Computation*, 66(217):157-192, 1997.

[CD18] René Carmona and François Delarue. *Probabilistic theory of mean field games with applications*. I-II, volume 83 of *Probability Theory and Stochastic Modelling*. Springer, Cham, 2018. Mean field FBSDEs, control, and games.

[CL18] Pierre Cardaliaguet and Charles-Albert Lehalle. Mean field game of controls and an application to trade crowding. *Mathematics and Financial Economics*, 12(3):335-363, 2018.

[CST22] Jean-François Chassagneux, Lukasz Szpruch, and Alvin Tse. Weak quantitative propagation of chaos via differential calculus on the space of measures. *The Annals of Applied Probability*, 32(3):1929-1969, 2022.

[DGLP06] Sylvain Delattre, Siegfried Graf, Harald Luschgy, and Gilles Pagès. Quantization of probability distributions under norm-based distortion measures ii: Self-similar distributions. *Journal of mathematical analysis and applications*, 318(2):507-516, 2006.

[dRES22] Gonçalo dos Reis, Stefan Engelhardt, and Greig Smith. Simulation of McKean-Vlasov SDEs with super-linear growth. *IMA Journal of Numerical Analysis*, 42(1):874-922, 2022.

[Edw11] David A Edwards. On the Kantorovich-Rubinstein theorem. *Expositiones Mathematicae*, 29(4):387-398, 2011.

[FG15] Nicolas Fournier and Arnaud Guillin. On the rate of convergence in Wasserstein distance of the empirical measure. *Probability Theory and Related Fields*, 162(3-4):707-738, 2015.

[G88] Jürgen Gärtner. On the McKean-Vlasov limit for interacting diffusions. *Mathematische Nachrichten*, 137:197-248, 1988.

[GL00] Siegfried Graf and Harald Luschgy. *Foundations of quantization for probability distributions*, volume 1730 of *Lecture Notes in Mathematics*. Springer-Verlag, Berlin, 2000.

[GPPP05] Emmanuel Gobet, Gilles Pagès, Huyên Pham, and Jacques Printems. *Discretization and simulation for a class of SPDEs with applications to Zakai and McKean-Vlasov equations*. Preprint LPMA, 2005.

[HK02] Rainer Hegselmann and Ulrich Krause. Opinion dynamics and bounded confidence models, analysis, and simulation. *Journal of artificial societies and social simulation*, 5(3), 2002.
[HL22] Marc Hoffmann and Yating Liu. A statistical approach for simulating the density solution of a McKean-Vlasov equation. *In preparation*, 2022.

[Hop50] Eberhard Hopf. The partial differential equation $u_t + uu_x = \mu u_{xx}$. *Communications on Pure and Applied Mathematics*, 3:201–230, 1950.

[Kie82] John C. Kieffer. Exponential rate of convergence for Lloyd’s method. I. *IEEE Transactions on Information Theory*, 28(2):205–210, 1982.

[Kie83] John C. Kieffer. Uniqueness of locally optimal quantizer for log-concave density and convex error weighting function. *IEEE Transactions on Information Theory*, 29(1):42–47, 1983.

[Lac18] Daniel Lacker. Mean field games and interacting particle systems. *Preprint*, 2018.

[Lem07] Vincent Lemaire. An adaptive scheme for the approximation of dissipative systems. *Stochastic Processes and their Applications*, 117(10):1491–1518, 2007.

[Liu19] Yating Liu. Optimal Quantization: Limit Theorems, Clustering and Simulation of the McKean-Vlasov Equation. Ph.D. thesis, Sorbonne Université, UPMC ; Laboratoire de Probabilités, Statistique et Modelisation (LPSM), 2019.

[LL18] Jean-Michel Lasry and Pierre-Louis Lions. Mean-field games with a major player. *Comptes Rendus Mathematique*, 356(8):886–890, 2018.

[Llo82] Stuart P. Lloyd. Least squares quantization in PCM. *IEEE Transactions on Information Theory*, 28(2):129–137, 1982.

[LP08] Harald Luschgy and Gilles Pagès. Functional quantization rate and mean regularity of processes with an application to Lévy processes. *The Annals of Applied Probability*, 18(2):427–469, 2008.

[LP20] Yating Liu and Gilles Pagès. Convergence rate of optimal quantization and application to the clustering performance of the empirical measure. *Journal of Machine Learning Research*, 21(86):1–36, 2020.

[LP22] Yating Liu and Gilles Pagès. Functional convex order for the scaled McKean-Vlasov processes. *The Annals of Applied Probability*, to appear, 2022.

[MA01] Nicolas Martzel and Claude Aslangul. Mean-field treatment of the many-body Fokker-Planck equation. *Journal of Physics A: Mathematical and General*, 34(50):11225, 2001.

[McK67] 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.

[Pag15] Gilles Pagès. Introduction to vector quantization and its applications for numerics. *CEM-RACS 2013—modelling and simulation of complex systems: stochastic and deterministic approaches*, 48:29–79, 2015.

[Pag18] Gilles Pagès. *Numerical probability*. Universitext. Springer, Cham, 2018. An introduction with applications to finance.

[Pol82] David Pollard. Quantization and the method of K-means. *IEEE Transactions on Information Theory*, 28(2):199–205, 1982.

[PP03] Gilles Pagès and Jacques Printems. Optimal quadratic quantization for numerics: the Gaussian case. *Monte Carlo Methods and Applications*, 9(2):135–165, 2003.

[PPP04] Gilles Pages, Huyên Pham, and Jacques Printems. An optimal markovian quantization algorithm for multi-dimensional stochastic control problems. *Stochastics and dynamics*, 4(04):501–545, 2004.
[PS15] Gilles Pagès and Abass Sagna. Recursive marginal quantization of the Euler scheme of a diffusion process. *Applied Mathematical Finance*, 22(5):463–498, 2015.

[PY16] Gilles Pagès and Jun Yu. Pointwise convergence of the Lloyd I algorithm in higher dimension. *SIAM Journal on Control and Optimization*, 54(5):2354–2382, 2016.

[Szn91] Alain-Sol Sznitman. Topics in propagation of chaos. In *École d'Été de Probabilités de Saint-Flour XIX—1989*, volume 1464 of *Lecture Notes in Math.*, pages 165–251. Springer, Berlin, 1991.

[Vil09] Cédric Villani. *Optimal transport, Old and new*, volume 338 of *Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]*. Springer-Verlag, Berlin, 2009.