An adaptive Euler–Maruyama scheme for McKean–Vlasov SDEs with super-linear growth and application to the mean-field FitzHugh–Nagumo model

Abstract

In this paper, we introduce fully implementable, adaptive Euler–Maruyama schemes for McKean–Vlasov stochastic differential equations (SDEs) assuming only a standard monotonicity condition on the drift and diffusion coefficients but no global Lipschitz continuity in the state variable for either, while global Lipschitz continuity is required for the measure component. We prove moment stability of the discretised processes and a strong convergence rate of $1/2$. Several numerical examples, centred around a mean-field model for FitzHugh–Nagumo neurons, illustrate that the standard uniform scheme fails and that the adaptive approach shows in most cases superior performance to tamed approximation schemes. In addition, we introduce and analyse an adaptive Milstein scheme for a certain sub-class of McKean–Vlasov SDEs with linear measure-dependence of the drift.

Introduction

A McKean–Vlasov equation (introduced by McKean [1]) for a $d$-dimensional process $X={\left(X_t\right)}_{t\in [0,T]}$, with a given $T>0$, is an SDE where the underlying coefficients depend on the current state $X_t$ and, additionally, on the law of $X_t$, i.e., they have the form

$$\mathrm{d}{X}_t=b(t,X_t,{\mathcal{L}}_{X_t})\mathrm{d}t+\sigma (t,X_t,{\mathcal{L}}_{X_t})\mathrm{d}{W}_t,X_0\in L_0^m\left({\mathbb{R}}^d\right),$$

where $W={\left(W_t\right)}_{t\in [0,T]}$ is a $k$-dimensional standard Brownian motion and ${\mathcal{L}}_{X_t}$ denotes the marginal law of the process $X$ at time $t$. Here, for $m\ge 2$, $L_0^m\left({\mathbb{R}}^d\right)$ denotes the space of ${\mathbb{R}}^d$-valued, ${\mathcal{F}}_0$-measurable random variables $Y$ satisfying $\mathbb{E}\left[{\left|Y\right|}^m\right]<\infty$.

The existence and uniqueness theory for strong solutions of McKean–Vlasov SDEs with coefficients of linear growth and Lipschitz conditions with respect to the state and the measure is well-established; see, e.g., [2] for the classical setting. For further specific existence and uniqueness results on weak and strong solutions of McKean–Vlasov SDEs we refer to [3], [4], [5] and references cited therein. Most relevant to this work, in the case of super-linear drift and diffusion, it is known from [6], [7] that McKean–Vlasov SDEs admit a unique strong solution under a monotonicity condition and Lipschitz continuity of the coefficients in the measure component.

Our motivation to study McKean–Vlasov SDEs with non-globally Lipschitz drift (and diffusion) is that several important models in practice involve mean-field terms and do not exhibit the classical global Lipschitz conditions on the coefficients of the SDE. For example, some models in neuroscience, such as the Hodgkin–Huxley and FitzHugh–Nagumo models [8], [9], or mean-field equations describing the behaviour of a (large) network of interacting spiking neurons [10], do not satisfy the classical assumptions (e.g., linear growth) on the drift coefficient.

The simulation of McKean–Vlasov SDEs of the form (1) typically involves two steps: a particle approximation and a time stepping-scheme. For a seminal work on numerical schemes for McKean–Vlasov SDEs we refer the reader to [11].

In the first step, at each time $t$, the true measure ${\mathcal{L}}_{X_t}$ is approximated by the empirical measure

$${\mu }_t^N\left(\mathrm{d}x\right)≔\frac{1}{N}\sum _{j=1}^N{\delta }_{X_t^{j,N}}\left(\mathrm{d}x\right),$$

where ${(X_t^{1,N},\dots ,X_t^{N,N})}_{t\in [0,T]}^T$ (so-called interacting particle system) is the solution to the ${\left({\mathbb{R}}^d\right)}^N$-dimensional SDE with components

$$\mathrm{d}{X}_t^{i,N}=b(t,X_t^{i,N},{\mu }_t^N)\mathrm{d}t+\sigma (t,X_t^{i,N},{\mu }_t^N)\mathrm{d}{W}_t^i,X_0^{i,N}=X_0^i,$$

and $(W^i,X_0^i)$, $i\in \{1,\dots ,N\}$, are independent Brownian motions (independent of $W$) and i. i. d. random initial values with $\text{Law}\left(X_0^i\right)={\mu }_0$, respectively. Furthermore, ${\delta }_x$ denotes the Dirac measure at $x$.

Key to the convergence as $N\to \infty$ is the concept of propagation of chaos. Consider $N$ (independent) processes $X^i$ (with initial value $X_0^i$) driven by independent Brownian motions $W^i$, $i\in \{1,\dots ,N\}$

$$\mathrm{d}{X}_t^i=b(t,X_t^i,{\mathcal{L}}_{X_t^i})\mathrm{d}t+\sigma (t,X_t^i,{\mathcal{L}}_{X_t^i})\mathrm{d}{W}_t^i,$$

which are standard McKean–Vlasov SDEs (and by uniqueness ${\mathcal{L}}_{X_t^i}={\mathcal{L}}_{X_t}$). Then pathwise propagation of chaos refers to the property

$$\underset{N\to \infty }{lim}\underset{i\in \{1,\dots ,N\}}{max}\mathbb{E}\left[\underset{t\in [0,T]}{sup}{|X_t^{i,N}-X_t^i|}^2\right]=0.$$

This type of result is classical under global Lipschitz conditions (see, e. g., [12]), and has been proven in [13, Proposition 3.1] under super-linear growth conditions of the drift. While the rates for the propagation of chaos shown in [13] suffer from the curse of dimensionality (i.e., the rate depends on $d$ and diminishes as $d\to \infty$)  [14, Lemma 5.1] gives for $p\le 4$ and any $i\in \{1,\dots ,N\}$ the bound

$$\mathbb{E}\left[\underset{t\in [0,T]}{sup}{|X_t^i-X_t^{i,N}|}^p\right]\le C{N}^{-p/2},$$

where the constant $C>0$ is dimension-independent. However, in [14, Lemma 5.1], the drift coefficient of the underlying McKean–Vlasov SDE is assumed to satisfy strong regularity assumptions (differentiability, along with Lipschitz and boundedness conditions on the state and measure components) and the diffusion coefficient is assumed to be constant, i.e., the framework is much more restrictive than the setting of the present article (in particular, the super-linear growth case is not covered in [14]).

As McKean–Vlasov equations typically arise from a mean-field approximation to a high- but finite-dimensional particle system, a particle approximation to the McKean–Vlasov equation simply reverses this step. To be more concrete, the FitzHugh–Nagumo model for $N$ interacting neurons from $P$ populations with different characteristics is a system of $N$ three-dimensional SDEs of the form

$$\mathrm{d}{X}_t^{i,N}=f_{\alpha }(t,X_t^{i,N})\mathrm{d}t+g_{\alpha }(t,X_t^{i,N})\left[\begin{array}{c}\hfill \mathrm{d}{W}_t^i\hfill \\ \hfill \mathrm{d}{W}_t^{i,y}\hfill \end{array}\right]+\sum _{\gamma =1}^P\frac{1}{N_{\gamma }}\sum _{j,p\left(j\right)=\gamma }\left(b_{\alpha \gamma }(X_t^{i,N},X_t^{j,N})\mathrm{d}t+{\beta }_{\alpha \gamma }(X_t^{i,N},X_t^{j,N})\mathrm{d}{W}_t^{i,\gamma }\right),$$

where $N={\sum }_{\gamma =1}^P{N}_{\gamma }$ and $f_{\alpha }$ satisfies a monotone growth condition but is not globally Lipschitz continuous (see Example 4 in Section 4.4 for a detailed description). It is shown in [8], [9] that the mean-field limit $N\to \infty$ is described by the three-dimensional McKean–Vlasov SDE

$$\mathrm{d}{\bar{X}}_t^{\alpha }=f_{\alpha }(t,{\bar{X}}_t^{\alpha })\mathrm{d}t+\sum _{\gamma =1}^P{\mathbb{E}}_{\bar{Z}}\left[b_{\alpha \gamma }({\bar{X}}_t^{\alpha },{\bar{Z}}_t^{\gamma })\right]\mathrm{d}t+g_{\alpha }(t,{\bar{X}}_t^{\alpha })\mathrm{d}{W}_t^{\alpha }+\sum _{\gamma =1}^P{\mathbb{E}}_{\bar{Z}}\left[{\beta }_{\alpha \gamma }({\bar{X}}_t^{\alpha },{\bar{Z}}_t^{\gamma })\right]\mathrm{d}{W}_t^{\alpha \gamma },$$

for $\alpha \in \{1,\dots ,P\}$, where $\bar{Z}$ is a vector-valued process independent of $\bar{X}$ with the same law. Conversely, (5) can be interpreted as a particle approximation to (6).

The focus of our paper is a time-stepping scheme for the particle system in (3) over some finite time horizon $[0,T]$. We have in mind systems which arise from particle approximations of McKean–Vlasov equations, but equally the scheme is of interest for interacting particle systems, potentially high-dimensional ones, in their own right. Note that at each time-step, the computational cost of simulating the particle system is usually $\mathcal{O}\left(N^2\right)$, as interaction terms due to the dependence of the coefficients on the empirical measure have to be computed for each particle. Strategies to reduce this cost include the projected particle method of [15].

In the context of classical SDEs, it is well-documented that the explicit Euler–Maruyama scheme (see [16]), for which strong convergence results of order $1/2$ are known under Lipschitz-type conditions, is generally not appropriate if we work with drift terms of super-linear growth. Specifically, it was shown in [17] for a certain SDE that the uniform time-step Euler–Maruyama scheme is not stable, i.e., the moments of the discretised process explode as the mesh-size tends to zero, even though a unique solution of the original SDE with bounded moments exists. A similar phenomenon is observed in the context of particle systems in [13], namely that during the simulation some of the particles strongly oscillate and eventually diverge (as also illustrated by Fig. 5(a) in our paper for the FitzHugh–Nagumo model).

To overcome this problem in the setting of classical SDEs, several stable time-discretisations have been introduced, including tamed explicit Euler and Milstein schemes [18], [19], [20], balanced schemes [21], [22], an explicit adaptive Euler–Maruyama method [23], a truncated Euler method [24] and an implicit Euler scheme [25]. For these methods, stability of the discretised process and strong convergence results have been proven. In [26], an explicit and a semi-implicit adaptive scheme for classical SDEs with non-globally Lipschitz continuous coefficients is studied, assuming a monotonicity property resulting from a Khasminskii-type condition [27], and additionally differentiability of the coefficients. A backstop scheme with a strict upper and lower bounds on the time-step immediately gives attainability of the final time $T$, and convergence arbitrarily close to order 1/2 is shown for a choice of the time-step which depends solely on the growth of the drift term. An adaptive time-stepping Milstein method for SDEs with one-sided Lipschitz drift was developed in [28], which also relies on such a backstop scheme.

In [13], the tamed Euler scheme and an implicit scheme were introduced for McKean–Vlasov SDEs with super-linearly growing drift and a diffusion coefficient which is globally Lipschitz in the state variable. Both coefficients are assumed to be Lipschitz continuous in the measure. In this setting, besides pathwise propagation of chaos, the paper establishes strong convergence of the time-discretised particle system to the SDE (3). The recent papers [29], [30], [31] introduce and analyse tamed Milstein schemes for (delay) McKean–Vlasov SDEs with a one-sided Lipschitz condition on the drift and a globally Lipschitz continuous diffusion coefficient (in both, the state and measure component), while [7] extends this analysis to the framework of McKean–Vlasov SDEs with common noise, where the diffusion coefficient is also allowed to grow super-linearly. We refer to [7] for examples of equations exhibiting such super-linear diffusions.

The first contribution of the current article is to introduce implementable adaptive time-stepping methods for McKean–Vlasov SDEs with super-linear growth in the drift and diffusion, assuming only a monotonicity condition (see (H.3(2)) in Section 3.1). The adaptive scheme has to cope with the fact that different time meshes will generally be used for different particles and the mean-field term needs to be approximated efficiently on all such time meshes, a difficulty not encountered for standard SDEs.

First, we will present an explicit adaptive Euler–Maruyama scheme for the setting where the diffusion coefficient is globally Lipschitz in both components. Here, we need only a mild assumption on the choice of the adaptive time-step (in the spirit of [23]) in order to achieve stability of the scheme. Then, we remove the Lipschitz assumption on the diffusion coefficient and impose only a certain Khasminskii-type monotonicity condition. However, we will need stronger assumptions on the time-step function, in order to also control the growth in the diffusion coefficient. The original analysis is not applicable in this more general setting and a new time-step control with a suitable stopping time argument is needed, which ultimately gives attainability of the final time and stability of the scheme. This extension is achieved without the need for a (semi-implicit) backstop scheme if the proposed adaptive time-step becomes too small and equally forms a contribution to the existing literature on adaptive time-stepping schemes for classical (non-measure dependent) SDEs.

Moreover, we prove moment stability and strong convergence of order $1/2$ in the time-step. Several numerical examples in Section 4 demonstrate that our adaptive scheme often outperforms the tamed Euler scheme introduced in [13], i.e., in many cases it gives more accurate numerical approximations or even achieves superior strong convergence rates.

Additionally, we present an adaptive Milstein scheme for a sub-class of McKean–Vlasov SDEs where the drift depends on the state and linearly on its law but the diffusion coefficient is only state dependent, see Eq. (16) for details. We consider this special class firstly because many relevant examples arising in practical applications have this form (see, e.g., [32] and the references cited therein), and secondly to avoid introducing a notion of derivatives for functions depending on probability measures (e.g., the Lions derivative; see [30], [31], [33] for details). We do not foresee any fundamental difficulties in extending the analysis to a more general setting, and leave it to future work.

A new difficulty for the computation and the analysis of adaptive schemes for particle systems, as mentioned above, is the approximation of the empirical measure at all adaptive time points. If the number of particles is large, as for particle approximations of McKean–Vlasov equations, this may result in a high computational cost if the empirical measure is re-evaluated at the smallest time-steps dictated by any particle. To avoid this, we present a second scheme where the empirical measure is kept constant over larger, fixed time-steps. This procedure avoids this extra cost but still allows us to obtain the same stability and convergence results as for the more expensive scheme. We prove stability for both these Euler–Maruyama schemes in the super-linear drift case, and convergence for the second scheme, while the convergence proof for the first scheme is found in the arXiv version of this article. Similarly, we only give the stability proofs for super-linear diffusion and the Milstein scheme for the first variant; while the stability results for the second variant follow along similar lines as the results presented in this paper, the proofs of strong convergence are subject to future research.

The remainder of this article is organised as follows: In Section 2, we describe the particle method and our adaptive schemes, and give the main convergence result (concerned with a globally Lipschitz continuous diffusion); the proof thereof is deferred to Section 5. Section 3 introduces an adaptive Euler–Maruyama scheme for a more general class where we allow both coefficients to grow super-linearly, and an adaptive Milstein scheme for a more specific class of McKean–Vlasov equations with linear measure-dependence; the proofs are again given in Section 5. The numerical results for several standard test cases taken from the literature, notably the FitzHugh–Nagumo model, are presented in Section 4.

We end this section by introducing some basic notations and give background results needed throughout this article.

Preliminaries

Throughout this paper, let $({\mathbb{R}}^d,〈\cdot ,\cdot 〉,\left|\cdot \right|)$ represent the $d$-dimensional Euclidean space, and $(\Omega ,\mathcal{F},\mathbb{F}={\left({\mathcal{F}}_t\right)}_{t\in [0,T]},\mathbb{P})$, for a given $T>0$, will be a complete filtered probability space, where $\mathbb{F}$ is the natural filtration of the $k$-dimensional Brownian motion ${\left(W_t\right)}_{t\in [0,T]}$ augmented with an independent $\sigma$-algebra ${\mathcal{F}}_0$.

In addition, we use $\mathcal{P}\left({\mathbb{R}}^d\right)$ to denote the family of all probability measures on $({\mathbb{R}}^d,\mathcal{B}\left({\mathbb{R}}^d\right))$, where $\mathcal{B}\left({\mathbb{R}}^d\right)$ denotes the Borel $\sigma$-field over ${\mathbb{R}}^d$, and define the subset of probability measures with finite second moment by

$${\mathcal{P}}_2\left({\mathbb{R}}^d\right)≔\left\{\mu \in \mathcal{P}\left({\mathbb{R}}^d\right)|{\int }_{{\mathbb{R}}^d}{\left|x\right|}^2\mu \left(\mathrm{d}x\right)<\infty \right\}.$$

As metric on the space ${\mathcal{P}}_2\left({\mathbb{R}}^d\right)$, we use the Wasserstein distance. For $\mu ,\nu \in {\mathcal{P}}_2\left({\mathbb{R}}^d\right)$, the Wasserstein distance between $\mu$ and $\nu$ is defined as

$${\mathcal{W}}_2(\mu ,\nu )≔\underset{\pi \in \mathcal{C}(\mu ,\nu )}{inf}{\left({\int }_{{\mathbb{R}}^d\times {\mathbb{R}}^d}{|x-y|}^2\pi (\mathrm{d}x,\mathrm{d}y)\right)}^{1/2},$$

where $\mathcal{C}(\mu ,\nu )$ means all the couplings of $\mu$ and $\nu$, i.e., $\pi \in \mathcal{C}(\mu ,\nu )$ if and only if $\pi (\cdot ,{\mathbb{R}}^d)=\mu \left(\cdot \right)$ and $\pi ({\mathbb{R}}^d,\cdot )=\nu \left(\cdot \right)$. Further, for $p\ge 2$, ${\mathcal{S}}^p\left([0,T]\right)$ refers to the space of ${\mathbb{R}}^d$-valued continuous, $\mathbb{F}$-adapted processes, defined on the interval $[0,T]$, with finite $p$-th moments, i.e., processes ${\left(X_t\right)}_{t\in [0,T]}$ satisfying $\mathbb{E}\left[{sup}_{t\in [0,T]}{\left|X_t\right|}^p\right]<\infty$.

For all linear operators (e.g., matrices $A\in {\mathbb{R}}^{d\times k}$) appearing in this article, we will use the standard Hilbert–Schmidt norm denoted by $\Vert \cdot \Vert$.