Study of micro–macro acceleration schemes for linear slow-fast stochastic differential equations with additive noise

Abstract

Computational multi-scale methods capitalize on a large time-scale separation to efficiently simulate slow dynamics over long time intervals. For stochastic systems, one often aims at resolving the statistics of the slowest dynamics. This paper looks at the efficiency of a micro–macro acceleration method that couples short bursts of stochastic path simulation with extrapolation of spatial averages forward in time. To have explicit derivations, we elicit an amenable linear test equation containing multiple time scales. We make derivations and perform numerical experiments in the Gaussian setting, where only the evolution of mean and variance matters. The analysis shows that, for this test model, the stability threshold on the extrapolation step is largely independent of the time-scale separation. In consequence, the micro–macro acceleration method increases the admissible time steps far beyond those for which a direct time discretization becomes unstable.

Appendices

A Monte-Carlo simulation

Let us discuss in more detail the discretisation in probability of the micro–macro acceleration algorithm presented in Sects. 2.1 and 4.2. We consider an initial distribution given as a random empirical measure

$$\begin{aligned} P^{\varDelta t,J}_{n}=\sum _{j=1}^{J}w^j_n\delta _{X^j_n}, \end{aligned}$$

where $X^j_n$, $j=1,\dotsc ,J$, are i.i.d. replicas with associated weights $w^j_n$.

Stage 1: Propagation of microscopic laws In this stage, we freeze the weights and propagate each replica over the K steps of inner integrator $\mathscr {S}^{\delta t}$ [compare with (4.4)]

$$\begin{aligned} X^{\delta t, j}_{n,k+1} = \mathscr {S}^{\delta t}(a, b, X^{\delta t, j}_{n,k},\xi ^{\delta t, j}),\quad X^{\delta t,j}_{n,0}=X^j_n, \end{aligned}$$

where $\xi ^{\delta t, j}$ are J i.i.d. centered normal variables with variance $\delta t$ (replicas of $\delta W_{k+1}$). In particular, when applying the Euler–Maruyama method to (2.1) we get

$$\begin{aligned} X^{\delta t, j}_{n,k+1} = X^{\delta t, j}_{n,k} + a\big (X^{\delta t, j}_{n,k}\big )\delta t + b\big (X^{\delta t, j}_{n,k}\big )\xi ^{\delta t, j}. \end{aligned}$$

The replicas $X^{\delta t, j}_{n,k}$ are associated to the following empirical measures

$$\begin{aligned} P^{\delta t,J}_{n,k} = \sum _{j=1}^{J}w^j_n\delta _{X^{\delta t,j}_{n,k}},\quad k=1,\dotsc ,K. \end{aligned}$$

Stage 2: Restriction to a finite number of observables In this stage, we evaluate the restriction operator $\mathscr {R}$ on the empirical measures to obtain K empirical observables as $\mathbf {m}^{J}_{n,k} = \mathscr {R}(P^{\delta t,J}_{n,k})$. For example, when restricting to the slow mean only as in (4.5), we find K vectors

$$\begin{aligned} \mathbf {m}^{s,J}_{n,k} = \sum _{j=1}^{J}w^j_nY^{\delta t,j}_{n,k}, \end{aligned}$$

where $Y^{\delta t,j}_{n,k}=\varPi ^sX^{\delta t,j}_{n,k}$ and $\varPi ^s$ is the projection onto the slow variables.

Stage 3: Extrapolation of macroscopic states We apply the extrapolation operator $\mathscr {E}^{\varDelta t - K\delta t}$ to the empirical macroscopic states $\mathbf {m}^{J}_{n,0},\dotsc ,\mathbf {m}^{J}_{n,K}$ and obtain $\mathbf {m}^{s,J}_{n+1}$ as in formula (2.4). In the special case of linear extrapolation one proceeds according to (4.7), with $\mathbf {m}^{s,J}_{n,0}$ and $\mathbf {m}^{s,J}_{n,K}$ plugged in the right-hand side.

Stage 4: Matching Finally, the matching amounts to re-weighting the replicas $X^{\delta t, j}_{n,K}$ using the Lagrange multipliers associated to procedure $(\hbox {ME})^{\mathrm{s}}$ or $(\hbox {MEV})^{\mathrm{s}}$ on page 12. Concentrating on $(\hbox {ME})^{\mathrm{s}}$ only, we approximate the Lagrange multipliers by applying the Newton–Raphson iteration to equation [compare with (2.7)]

$$\begin{aligned} \nabla _{\!\lambda }A^s\big (\lambda ,P^{\delta t,J}_{n,K}\big )=\mathbf {m}^{s,J}_{n+1},\quad \text {where}\quad A^s\big (\lambda ,P^{\delta t,J}_{n,K}\big ) = \ln \sum _{j=1}^{J}w^j_n\exp \big (\lambda \mathbin {\mathbf {\cdot }}\varPi ^sX^{\delta t,j}_{n,K}\big ). \end{aligned}$$

(A.1)

In this fashion, we obtain, up to a given tolerance, the vector $\overline{\lambda }^{s,J}_{n+1}$ of Lagrange multipliers with which, following (2.6), we evaluate weights

$$\begin{aligned} w^j_{n+1} = w^j_n\exp \Big (\overline{\lambda }^{s,J}_{n+1}\mathbin {\mathbf {\cdot }}\varPi ^sX^{\delta t,j}_{n,K} - A^s\big (\overline{\lambda }^{s,J}_{n+1},P^{\delta t,J}_{n,K}\big )\Big ). \end{aligned}$$

The matched empirical distribution reads

$$\begin{aligned} P^{\varDelta t,J}_{n+1} = \sum _{j=1}^{J}w^j_{n+1}\delta _{X^{\delta t,j}_{n,K}}. \end{aligned}$$

B Gap in the drift spectrum and the time scales of (3.2)

To elucidate how Assumption 3.1 influences the time scales present in the stochastic dynamics, recall first that Eq. (3.2) is related to the Ornstein–Uhlenbeck operator

$$\begin{aligned} \mathscr {L}= Ax\mathbin {\mathbf {\cdot }}\nabla _{\!x} + \frac{1}{2}\mathsf {Tr}(B\nabla ^2_{\!x}). \end{aligned}$$

(B.1)

SDE (3.2) and operator (B.1) are connected through the Markov semigroup $(e^{t\mathscr {L}})_{t\ge 0}$, generated by $\mathscr {L}$, that satisfy

$$\begin{aligned} \big (e^{t\mathscr {L}}f\big )(x) = \mathbb {E}_{}[f(X_t)|X_t=x], \end{aligned}$$

for every $t\ge 0$ and $f\in \mathscr {C}_{\mathrm {b}}(\mathbb {R}^d)$, the space of all continuous and bounded functions on $\mathbb {R}^d$. The assumptions that $\mathsf {Sp}(A)\subset \mathbb {C}_-$ and B is positive definite ensure the existence of a unique Gaussian invariant measure $\mathscr {N}_{0,V_\infty }$ for $(e^{t\mathscr {L}})_{t\ge 0}$, where $V_\infty$ is given by (3.4). To be more precise, the condition for invariance reads

$$\begin{aligned} \mathbb {E}_{\mathscr {N}_{0,V_\infty }}\big [e^{t\mathscr {L}}f\big ] = \mathbb {E}_{\mathscr {N}_{0,V_\infty }}[f], \end{aligned}$$

(B.2)

for all $f\in \mathscr {C}_{\mathrm {b}}(\mathbb {R}^d)$. Moreover, the semigroup $(e^{t\mathscr {L}})_{t\ge 0}$ extends to a strongly continuous semigroup of positive contractions in the complex Hilbert space $\mathscr {L}^2_{\mathbb {C}}(\mathbb {R}^d, d\mathscr {N}_{0,V_\infty })$ [34].

The time scales induced by the semigroup $(e^{t\mathscr {L}})_{t\ge 0}$ in the space $\mathscr {L}^2_{\mathbb {C}}(\mathbb {R}^d, d\mathscr {N}_{0,V_\infty })$ are determined by the eigenvalue problem $\mathscr {L}f = \gamma f$ [46, p. 371]. Every eigenpair $(\gamma ,\phi )$, with $\mathfrak {R}(\gamma )<0$ and $\Vert \phi \Vert _2=1$, is related to a decay of $e^{t\mathscr {L}}\phi$ towards the equilibrium on time scales of order $|2\gamma |^{-1}$, to wit

$$\begin{aligned} \mathbb {V}_{\!\mathscr {N}_{0,V_\infty }}\big [e^{t\mathscr {L}}\phi \big ] = \big \Vert e^{t\mathscr {L}}f - \mathbb {E}_{\mathscr {N}_{0,V_\infty }}[f]\big \Vert _2 = e^{2\mathfrak {R}(\gamma ) t}, \end{aligned}$$

where we used condition (B.2) and $\Vert \cdot \Vert _2$ denotes the associated $\mathscr {L}^2_{\mathbb {C}}$-norm. Having a complete orthonormal system $\{\phi _p\}_{p=1,2\dotsc }$ of eigenfunctions in $\mathscr {L}^2_{\mathbb {C}}(\mathbb {R}^d, d\mathscr {N}_{0,V_\infty })$, the Fourier expansion

$$\begin{aligned} f = \sum _{p=1}^{\infty }\langle f,\phi _p\rangle _2\,\phi _p \end{aligned}$$

(B.3)

decomposes the trend of $e^{t\mathscr {L}}f$ towards the equilibrium into separate modes, supported by the invariant subspaces generated by each $\phi _p$ and exponentially decaying with rates given by $\mathfrak {R}(\gamma _p)$.

The spectrum of the Ornstein–Uhlenbeck operator (B.1) in $\mathscr {L}^2_{\mathbb {C}}(\mathbb {R}^d, d\mathscr {N}_{0,V_\infty })$ consists of all complex numbers $\gamma = \sum _{i=1}^{d}n_i\kappa _i$ with $n_i\in \mathbb {N}_0, \kappa _i\in \mathsf {Sp}(A)$, and all the (generalised) eigenfunctions are polynomials and form a complete system [35]. Moreover, the invariant subspaces related to each $\gamma$ consist of homogeneous polynomials in the variables induced by the spectral decomposition of A and with degrees ranging throughout all $n_i>0$ that appear in the sums generating $\gamma$ [35, Sect. 4].

There are three main implications of these facts in our setting. First, the spectrum of $\mathscr {L}$ is independent of the diffusion matrix B. This demonstrates that B has no effect on the time scales of the dynamics and justifies the omission of any assumptions on its spectrum. Second, the eigenvalues of the drift matrix are embedded inside the spectrum of $\mathscr {L}$ and induce the most prominent time scales. Indeed, every $\gamma =\kappa _i$ is uniquely determined by $n_i=1$ and $n_j=0$ for $j\ne i$. The associated eigenfunction is a homogeneous polynomial of degree 1 in a variable associated to the invariant subspace of $\kappa _i$. Therefore, in the Fourier expansion (B.3), all eigenvalues $\kappa _i$ constitute the first approximation of f. Finally, the dynamics of $e^{t\mathscr {L}}$ has infinite number of different time scales but the gap in $\mathsf {Sp}(A)$ reveals itself at the lowest order modes. This can be seen by applying the decomposition of the state space $\mathbb {R}^d$ into the slow vector variable y, associated to $\varOmega ^s$, and the fast vector variable z, associated to $\varOmega ^f$, see the paragraph following Assumption 3.1. This grouping decomposes the first approximation of f into polynomials in y and z that equilibrate under the action of $e^{t\mathscr {L}}$ on two different time scales with gap given by the gap in $\mathsf {Sp}(A)$.