An explicit Milstein-type scheme for interacting particle systems and McKean–Vlasov SDEs with common noise and non-differentiable drift coefficients

For a given $T>0$, consider the following stochastic differential equation (SDE) of McKean–Vlasov type and with common noise,

almost surely for all $t\in[0,T]$, where $W:=\{W_{t}\}_{t\in[0,T]}$ and $W^{0}:=\{W_{t}^{0}\}_{t\in[0,T]}$ are, respectively, $m_{1}$ and $m_{0}$ dimensional independent Wiener processes and $\{\mathcal{L}^{1}(X_{t})\}_{t\in[0,T]}$ denotes the stochastic flow of conditional marginal laws of $X:=\{X_{t}\}_{t\in[0,T]}$ given $W^{0}$. The initial value $X_{0}$ is an $\mathscr{F}^{0}$- measurable random variable, independent of $W$ and $W^{0}$. The McKean–Vlasov SDE (1) can be viewed as an infinite-dimensional system of particles with $W$ representing the randomness inherent in the individual particle and $W^{0}$ the randomness common to all the particles. When $\sigma_{0}\equiv 0$, the particles are governed by only one source of randomness, $W$, and the stochastic flow $\{\mathcal{L}^{1}(X_{t})\}_{t\in[0,T]}$ becomes a deterministic one. Notice that McKean–Vlasov SDEs are different from standard SDEs due to the dependence of the coefficients on the (conditional) marginal law $\mathcal{L}^{1}(X_{t})$ of $X_{t}$ given $W^{0}$, which brings additional difficulties.

Due to their wide applications in areas such as Finance, mathematical neuroscience and biology, machine learning and physics — animal swarming, cell movement induced by chemotaxis, opinion dynamics, particle movement in porous media and electrical battery modelling, self-assembly of particles and dynamical density functional theory (see for example [31, 16, 8, 25, 2, 27, 33, 9, 13, 14, 28, 32])— McKean–Vlasov equations and associated interacting particle systems, with or without common noise, addressed via stochastic systems or associated Fokker Plank equations ([38, 22, 21]) have gained immense popularity.

As in the case of SDEs, explicit solutions of McKean–Vlasov SDEs are typically not available, which necessitates the development of numerical schemes to approximate them. The numerical approximation of McKean–Vlasov SDEs can be carried out in two steps, as explained below.

• •

  As a first step one builds the so-called interacting particle system, $\{X^{i,N}\}_{i\in\{1,\ldots,N\}}$, where one replaces the (conditional) marginal law appearing in the coefficients of (1) by the empirical law obtained from the particles. Concretely, taking $N$ i.i.d. copies $\{W^{i}\}_{i\in\{1,\ldots,N\}}$ of $W$ and $\{X^{i}_{0}\}_{i\in\{1,\ldots,N\}}$ of $X_{0}$,

  one defines the interacting particle system associated with the above McKean–Vlasov SDE by

  	$\displaystyle X_{t}^{i,N}=X_{0}^{i}$	$\displaystyle+\int_{0}^{t}b(s,X_{s}^{i,N},\mu_{s}^{X,N})\,ds+\sum_{\ell=1}^{m_{1}}\int_{0}^{t}\sigma_{1}^{\ell}(s,X_{s}^{i,N},\mu_{s}^{X,N})\,dW_{s}^{i,\ell}$
  		$\displaystyle+\sum_{\ell=1}^{m_{0}}\int_{0}^{t}\sigma_{0}^{\ell}(s,X^{i,N}_{s},\mu_{s}^{X,N})\,dW_{s}^{0,\ell},$		(2)

  almost surely for any $t\in[0,T]$ and $i\in\{1,\ldots,N\}$, where

  $\displaystyle\mu_{t}^{X,N}(\cdot):=\frac{1}{N}\sum_{i=1}^{N}\delta_{X_{t}^{i,N}}(\cdot)$

  is the empirical measure of $N$ particles. Subsequently, one needs to show, roughly put, that $\textrm{Law}(X^{i,N}_{\cdot})$ for some $i$ (fixed) convergences to $\textrm{Law}(X_{\cdot})$ of (1) as $N\to\infty$; see the seminal work by Sznitman [53] (and Proposition 2.5 below).

• •

  In the second step of approximation, the temporal discretization of the interacting particle system is performed to obtain fully implementable numerical schemes for the McKean–Vlasov SDEs such as Euler-type schemes and Milstein-type schemes. The main difficultly is to show that all estimates are independent of the number of particles $N$ (see Theorem 3.3 and Corollary 3.4 below). The direct application of results from classic SDE theory do not deliver this independence.

It is critical to note that the results of this manuscript double for either the numerical approximation of McKean–Vlasov SDEs if one’s starting point is (1), or for stand-alone systems of interacting particle SDE systems if one’s starting point is (2).

Results on the strong well-posedness and propagation of chaos for McKean–Vlasov SDEs are being extensively researched and we cannot possibly do justice to that growing body of literature. Nonetheless, we mention some of the milestones and more recent results fitting thematically with our manuscript: for PoC results, one starts from Sznitman’s seminal work [53] to the monographs [13, 14] and the recent developments by [23, 39, 40] to mention a few – overall, there are still gaps in the existing PoC results, especially across dimensions. Concretely, in [20, Section 3.5] the PoC rates across the dimension $d$ is estimated numerically for a diffusion $\sigma$ of polynomial growth and a drift $b$ of polynomial growth (a setting outside the scope of this work) and the rates estimated are better than the rates found in any theoretical results at present. For a focus on well-posedness, one starts from the McKean’s seminal work [43], to again [13, 14] and the recent developments [6, 17, 37, 44, 45, 52, 20] (and references therein).

The second step of numerical approximation of McKean–Vlasov SDEs saw relatively little development after the first Euler-type scheme was proposed and analysed (in a weak sense) in [10], but experience rapid advances in a string of recent papers, [3, 4, 5, 18, 19, 36, 37, 41, 24, 49]. In particular, [24] proposed an Euler-type numerical scheme for the interacting particle system associated with the McKean–Vlasov SDE, and its strong convergence is investigated when the drift coefficient grows super-linearly in the state variable. More precisely, the drift is assumed to be one-sided Lipschitz continuous in the state variable and Lipschitz continuous in the measure variable, while the diffusion coefficient is assumed to be Lipschitz continuous in both the state and measure variables. The rate of strong convergence of the scheme is shown to be equal to $1/2$. In [4, 36], a Milstein-type scheme for the interacting particle system associated with the McKean–Vlasov SDEs is proposed using the notion of Lions derivatives, introduced by P.-L. Lions in his lectures at the Collège de France and presented in [12], and its strong convergence is shown with rate $1$. The drift coefficient is assumed to satisfy a one-sided Lipschitz condition in the state variable and a polynomial Lipschitz condition, and the diffusion coefficient to satisfy Lipschitz condition; both are assumed to be Lipschitz continuous in the measure variable. Furthermore, only the diffusion coefficient is required to be once differentiable (in both state and measure variables). In [4], the authors additionally require second order differentiability of the coefficients. In [37], an Euler-type scheme and a Milstein-type scheme are developed for the interacting particle system connected with McKean–Vlasov SDEs with common noise, where all the coefficients are allowed to grow super-linearly in the state variable.

In this article, we develop a Milstein-type scheme for the interacting particle system corresponding to the McKean–Vlasov SDE without assuming the differentiability of the drift coefficient (in space or measure component), which is therefore more relaxed than the corresponding results in [4, 36, 37]. It is out of the scope of this work to lift the differentiability conditions on the diffusion coefficient and hence those assumptions match those already existing in the most recent literature. The relaxation of the regularity requirement of the drift coefficient is achieved by a certain randomisation strategy that needs to be applied to both state and measure components: the technical developments necessary to deal with this difficulty are the second contribution of this manuscript. In the case of SDEs (when the coefficients do not depend on the law of the solution process), the technique of randomisation has been studied in [35] (also [7, 34] and and more recently in [46, 48]) to construct a Milstein-type scheme without assuming the first order differentiability of the drift coefficient. As the coefficients in our settings depend on the law of the solution process as well, we require a two-fold randomisation – one with respect to the state variable and the other with respect to the measure variable. For this, we use a uniform random variable to generate a random point in each sub-interval of the time mesh and the Euler scheme is used to obtain values of the particles’ states at these random points, which are then used in the drift coefficient of the Milstein scheme, both for the state and empirical measure. The precise details of this randomisation can be found in Section 3. Critically, the technique developed in [35] for the analysis cannot be used directly in our settings and a novel approach is required. We observe the following.

• •

  The interacting particle system associated with the McKean–Vlasov SDE can be treated as an $(\mathbb{R}^{d})^{N}$-dimensional SDE and thus the results of [35] could be applied directly. However, all estimates would depend on $N$ and hence “explode” as $N$ tends to infinity. This implies that to establish our results a new tool must be developed in order to show the independence on $N$.

• •

  We propose a new notion of bistability and consistency of the numerical scheme that is appropriate for the context of high-dimensional interacting particle systems. Inspired by [35] we propose suitable stochastic Spijker norms capable of dealing with the interaction component of the particles. Further details can be found in Section 5.

• •

  A discussion on the practicalities of implementing our scheme is given in Section 3.1.1 which includes a critical view on complexity and the consequence of having common-noise.

• •

  In the simplest version possible of (1), three drift functions are well within the scope of our work are (linear interaction, convolution-functionals and linear interaction kernels)

  	$\displaystyle b(s,X_{s},\mathcal{L}^{1}(X_{s}))$	$\displaystyle=f_{1}(X_{s})+f_{2}\big{(}\tilde{\mathbb{E}}^{1}(g(X_{s})\big{)},$
  	$\displaystyle\widehat{b}(s,X_{s},\mathcal{L}^{1}(X_{s}))$	$\displaystyle=f_{1}(X_{s})+\int_{\mathbb{R^{d}}}g(X_{s}-y)\mathcal{L}^{1}(X_{s})(dy),$
  	$\displaystyle\tilde{b}(s,X_{s},\mathcal{L}^{1}(X_{s}))$	$\displaystyle=f_{1}(X_{s})+\int_{\mathbb{R^{d}}}\tilde{g}(X_{s},y)\mathcal{L}^{1}(X_{s})(dy),$

  for any real-valued functions $f_{1},f_{2},g,\tilde{g}$ that satisfy a standard Lipschitz condition in space (but are not differentiable), e.g., $g(x)=-|x|$ and more complex examples for $g,\tilde{g}$ can be found in [16, 13, 14, 30, 32] and $\tilde{\mathbb{E}}^{1}$ represents the conditional expectation given the common noise $W^{0}$. The 2nd drift example, $\hat{b}$, corresponds to the usual convolution operator fairly common in modelling with McKean–Vlasov SDE and associated interacting particle systems (with or without common noise) [15, 30, 32, 1]. Lastly, we point the reader to Example 3.15 in [1] that intuitively highlights why $\hat{b}$ is Lipschitz in the Wasserstein metric but is not Lions differentiable.

Let $T>0$ be a fixed constant. Consider probability spaces $(\Omega^{1},\mathscr{F}^{1},\mathbb{P}^{1})$ and $(\Omega^{0},\mathscr{F}^{0},\mathbb{P}^{0})$ equipped with filtrations $\mathbb{F}^{1}:=\{\mathscr{F}_{t}^{1}\}_{t\in[0,T]}$ and $\mathbb{F}^{0}:=\{\mathscr{F}_{t}^{0}\}_{t\in[0,T]}$, respectively. The filtrations $\mathbb{F}^{1}$ and $\mathbb{F}^{0}$ satisfy the usual conditions, i.e., they are complete and right continuous. Assume that $W:=\{W_{t}{\in\mathbb{R}^{m_{1}}}\}_{t\in[0,T]}$ and $W^{0}:=\{W_{t}^{0}{\in\mathbb{R}^{m_{0}}}\}_{t\in[0,T]}$ are independent Brownian motions defined on $(\Omega^{1},\mathscr{F}^{1},\mathbb{F}^{1},\mathbb{P}^{1})$ and $(\Omega^{0},\mathscr{F}^{0},\mathbb{F}^{0},\mathbb{P}^{0})$, respectively. In what follows, the interacting particles are governed by i.i.d. copies of $W$ and $W^{0}$ represents the noise common to all the particles. Let us define a product probability space $(\tilde{\Omega},\tilde{\mathscr{F}},\tilde{\mathbb{P}})$, where $\tilde{\Omega}:=\Omega^{1}\times\Omega^{0}$, $(\tilde{\mathscr{F}},\tilde{\mathbb{P}})$ is the completion of $(\mathscr{F}^{1}\otimes\mathscr{F}^{0},\mathbb{P}^{1}\otimes\mathbb{P}^{0})$ and $\tilde{\mathbb{F}}:=\{\tilde{\mathscr{F}}_{t}\}_{t\in[0,T]}$ is the completion and right-continuous augmentation of $\{\mathscr{F}_{t}^{1}\otimes\mathscr{F}_{t}^{0}\}_{t\in[0,T]}$. The expectation with respect to $\tilde{\mathbb{P}}$ is denoted by $\tilde{\mathbb{E}}$.

Let $b:[0,T]\times\mathbb{R}^{d}\times\mathscr{P}_{2}(\mathbb{R}^{d})\mapsto\mathbb{R}^{d}$, $\sigma_{1}:[0,T]\times\mathbb{R}^{d}\times\mathscr{P}_{2}(\mathbb{R}^{d})\mapsto\mathbb{R}^{d\times m_{1}}$ and $\sigma_{0}:[0,T]\times\mathbb{R}^{d}\times\mathscr{P}_{2}(\mathbb{R}^{d})\mapsto\mathbb{R}^{d\times m_{0}}$ be $\mathscr{B}([0,T])\otimes\mathscr{B}(\mathbb{R}^{d})\otimes\mathscr{B}(\mathscr{P}_{2}(\mathbb{R}^{d}))$-measurable functions.

In this article, we consider the following $\mathbb{R}^{d}$-valued McKean–Vlasov stochastic differential equations (SDEs) with common noise defined on $(\tilde{\Omega},\tilde{\mathscr{F}},\tilde{\mathbb{F}},\tilde{\mathbb{P}})$,

almost surely for all $t\in[0,T]$, where $\{\mathcal{L}^{1}(X_{t})\}_{t\in[0,T]}$ is the flow of conditional laws of $X_{t}$ given $W^{0}$ and $X_{0}$. A priori, it is not certain that the flow of conditional marginals $\{\mathcal{L}^{1}(X_{t})\}_{t\in[0,T]}$ is an $\mathbb{F}^{0}$-adapted continuous process. However, due to Lemma 2.5 in [14], when the McKean–Vlasov SDE (2) has a unique $\tilde{\mathbb{F}}$-adapted continuous solution with uniformly bounded second moment, then $\{\mathcal{L}^{1}(X_{t})\}_{t\in[0,T]}$ is an $\mathbb{F}^{0}$-adapted continuous process.

We make the following assumptions.

The proof of the following proposition can be found in [14, 37] and Appendix A.1.

To introduce the interacting particle system connected with the McKean–Vlasov SDEs (2), let us consider $N\in\mathbb{N}$ i.i.d. copies of $W$ and $X_{0}$, denoted by $W^{i}$ and $X_{0}^{i}$ for $i\in\{1,\ldots,N\}$, respectively. Define the following system of equations,

almost surely for any $t\in[0,T]$ and $i\in\{1,\ldots,N\}$. Notice that due to Proposition 2.11 in [14], $\mathbb{P}^{0}\big{(}\mathcal{L}^{1}(X_{t}^{1})=\mathcal{L}^{1}(X_{t}^{i})\mbox{ for all }t\in[0,T]\big{)}=1$. We remark that the proof of Proposition 2.11 in [14], which establishes the result in a more restrictive setting, uses (only) the well-posedness of the system (2) and Theorem 1.33 from [14] (Yamada Watanabe Theorem) and thus particles have the same law under our settings also. The system (2) is popularly known as the conditional non-interacting particle system. On approximating $\mathcal{L}^{1}(X_{t}^{1})$ by the empirical measure of the states $\{X_{t}^{i}\}_{i\in\{1,\ldots,N\}}$ of $N$ particles at time $t$, and denoted $\mu_{t}^{X,N}$, one obtains the following system of interacting particles,

almost surely for any $t\in[0,T]$ and $i\in\{1,\ldots,N\}$, where

is the empirical measure of $N$ particles.

The following lemma gives the time-regularity of the interacting particle system (2) and its proof is given in Appendix A.2.

The convergence of the interacting particle system (2) to the non-interacting particle system (2) is popularly known in the literature as the propagation of chaos and is stated in the following proposition (see Theorem 2.12 in [14] for details or Appendix A.3). For this, let us define the empirical measure of the non-interacting particle system (2) as

almost surely for any $t\in[0,T]$ and $N\in\mathbb{N}$. Further, use Theorem 5.8 in [13] and Proposition 2.2 to obtain the following estimate,

for any $t\in[0,T]$ and $N\in\mathbb{N}$, where $C_{3}:=C_{3}(d,\bar{p},\|X_{0}\|_{\mathscr{L}^{\bar{p}}(\tilde{\Omega})}^{\bar{p}})$ is a positive constant.

The proof of this result can be found in Appendix A.3.

All in all, the PoC rates here are in line with the PoC results found in the numerics for McKean Vlasov SDE literature [3, 4, 5, 18, 19, 52, 36, 37, 41, 24, 49]. There are improved results obtaining rate of $1/N$ and $1/N^{2}$ (e.g., [23, 39]) that hold under stricter conditions that do not fit the scope of the work presented in the manuscript. We have added text to this effect in the main body of the paper.

The drift-randomised Milstein scheme was originally proposed and analysed for standard SDEs with non-differentiable drift coefficient in [35]. In the case of McKean–Vlasov SDEs, for which the coefficients depend on the law of the solution process as well, one needs additional randomisation of the drift coefficient with respect to the measure component and a way to deal with the implications from interacting particles. The associated difficulties are tackled in this paper.