Convergence and stability of a micro-macro acceleration method:linear slow-fast stochastic differential equations with additive noise

TL;DR

This paper analyzes the convergence and stability of a micro-macro acceleration method for simulating stiff linear stochastic differential equations with time-scale separation, focusing on additive noise and extrapolation of the slow mean.

Contribution

It provides a rigorous convergence proof for Gaussian initial distributions and a stability analysis for non-Gaussian laws in the context of the micro-macro acceleration algorithm.

Findings

Proves convergence to microscopic dynamics for Gaussian initial distributions.

Establishes stability conditions for non-Gaussian initial laws.

Analyzes the method's behavior on linear SDEs with additive noise.

Abstract

We analyse the convergence and stability of a micro-macro acceleration algorithm for Monte Carlo simulations of stiff stochastic differential equations with a time-scale separation between the fast evolution of the individual stochastic realizations and some slow macroscopic state variables of the process. The micro-macro acceleration method performs a short simulation of a large ensemble of individual fast paths, before extrapolating the macroscopic state variables of interest over a larger time step. After extrapolation, the method constructs a new probability distribution that is consistent with the extrapolated macroscopic state variables, while minimizing Kullback-Leibler divergence with respect to the distribution available at the end of the Monte Carlo simulation. In the current work, we study the convergence and stability of this method on linear stochastic differential equations with additive noise, when only extrapolating the mean of the slow component. For this case, we prove convergence to the microscopic dynamics when the initial distribution is Gaussian and present a stability result for non-Gaussian initial laws.

Equations

$$d{X}_t=A{X}_{t}dt+Bd{W}_t$$

$${\mathbb{R}}^d={\mathbb{R}}^{d_s}\oplus {\mathbb{R}}^{d-d_s},A={\Pi }^s{A}^s\oplus {\Pi }^f{A}^f.$$

$$\mathcal{D}(Q||P)={\mathbb{E}}_Q\left[\ln \left(\frac{dQ}{dP}\right)\right],$$

$$\mathcal{M}({\bar{\mu }}^s,P)=\underset{Q\in \mathcal{P}}{\text{arg min}}\mathcal{D}(Q||P),\text{s.t.}{\mathbb{E}}_Q[{\Pi }^s]={\bar{\mu }}^s.$$

$$\mathchoice{\frac{\operatorname{d}\!\bar{Q}}{\operatorname{d}\!P}}{\operatorname{d}\!\bar{Q}/\!\operatorname{d}\!P}{\operatorname{d}\!\bar{Q}/\!\operatorname{d}\!P}{\operatorname{d}\!\bar{Q}/\!\operatorname{d}\!P}(x)=\exp\!\big{(}\bar{\lambda}^{s}\cdot x-A(\bar{\lambda}^{s},P)\big{)},$$

$${\nabla }_{{\lambda }_s}A({\bar{\lambda }}^s,P)={\bar{\mu }}_s.$$

$${\pi }_Q(x)=\exp \left({\bar{\lambda }}^s\cdot x-A\left({\bar{\lambda }}^s,P\right)\right){\pi }_P(x).$$

$$\frac{d\bar{Q}}{dP}(y,z)=\frac{d{\bar{Q}}^s}{d{P}^s}(y),y\in {\mathbb{R}}^{d_s},z\in {\mathbb{R}}^{d_f},$$

$${\mathbb{E}}_{\bar{Q}}[f(Z)]={\int }_{{\mathbb{R}}^{d_s}}{\mathbb{E}}_P[f(Z)|Y=y]{\bar{Q}}^s(dy).$$

$$\displaystyle\begin{aligned} A(\lambda,P)&=\ln\int_{\mathbb{R}^{d}}\exp\!\big{(}\lambda\mathbin{\mathbf{\cdot}}\Pi^{s}x\big{)}\,P(\operatorname{d}\!{x})\\[1.99997pt] &=\ln\!\Big{\{}\int_{\mathbb{R}^{d_{s}}}\exp\!\big{(}\lambda\mathbin{\mathbf{\cdot}}y\big{)}\int_{\mathbb{R}^{d_{f}}}P^{f|s}(\operatorname{d}\!{z}|y)\,P^{s}(\operatorname{d}\!{y})\Big{\}}\\[1.99997pt] &=\ln\mathbb{E}_{P^{s}}\big{[}\exp(\lambda\mathbin{\mathbf{\cdot}}\Pi^{s})]=A(\lambda,P^{s}),\end{aligned}$$

$${\mathbb{E}}_{\bar{Q}}[f(Z)]={\int }_{{\mathbb{R}}^d}f(z)\bar{Q}(dy,dz)={\int }_{{\mathbb{R}}^d}f(z)\frac{d\bar{Q}}{dP}(y,z)P(dy,dz).$$

$$\displaystyle\begin{aligned} \mathbb{E}_{\bar{Q}}[f(Z)]&=\int_{\mathbb{R}^{d_{s}}}\int_{\mathbb{R}^{d_{f}}}f(z)\mathchoice{\frac{\operatorname{d}\!\bar{Q}^{s}}{\operatorname{d}\!P^{s}}}{\operatorname{d}\!\bar{Q}^{s}/\!\operatorname{d}\!P^{s}}{\operatorname{d}\!\bar{Q}^{s}/\!\operatorname{d}\!P^{s}}{\operatorname{d}\!\bar{Q}^{s}/\!\operatorname{d}\!P^{s}}(y)\,P^{f|s}(\operatorname{d}\!{z}|y)\,P^{s}(\operatorname{d}\!{y})\\ &=\int_{\mathbb{R}^{d_{s}}}\Big{[}\int_{\mathbb{R}^{d_{f}}}f(z)\,P^{f|s}(\operatorname{d}\!{z}|y)\Big{]}\mathchoice{\frac{\operatorname{d}\!\bar{Q}^{s}}{\operatorname{d}\!P^{s}}}{\operatorname{d}\!\bar{Q}^{s}/\!\operatorname{d}\!P^{s}}{\operatorname{d}\!\bar{Q}^{s}/\!\operatorname{d}\!P^{s}}{\operatorname{d}\!\bar{Q}^{s}/\!\operatorname{d}\!P^{s}}(y)P^{s}(\operatorname{d}\!{y}),\end{aligned}$$

$$X_{n,k+1}=\left(I+A\delta t\right)X_{n,k}+\sqrt{\delta t}B\delta W_{n,k},$$

$${\mu }_{n,k}^s=\mathbb{E}[{\Pi }^s{P}_{n,k}],k=0,\dots ,K.$$

$${\mu }_{n+1}^s={\mu }_n^s+\frac{\Delta t}{K\delta t}\left({\mu }_{n,K}^s-{\mu }_n^s\right).$$

$$P_{n+1}=\mathcal{M}\left({\mu }_{n+1}^s,P_{n,K}\right),$$

$$\mu =\left[\begin{array}{c}{\mu }^s\\ {\mu }^f\end{array}\right],\Sigma =\left[\begin{array}{cc}{\Sigma }^s& C\\ C^T& {\Sigma }^f\end{array}\right].$$

$$\bar{\mu}^{f}=\int_{\mathbb{R}^{d_{s}}}\big{\{}\mu^{f}+C^{T}(\Sigma^{s})^{-1}(y-\mu^{s})\big{\}}\,\bar{Q}^{s}(dy)=\mu^{f}+C^{T}(\Sigma^{s})^{-1}(\bar{\mu}^{s}-\mu^{s}).$$

$${\bar{\Sigma }}^f={\int }_{{\mathbb{R}}^{d_s}}{\mathbb{E}}_P[(Z-{\bar{\mu }}^f)(Z-{\bar{\mu }}^f{)}^T|Y=y]{\bar{Q}}^s(dy).$$

$$\text{centering}{\bar{\Sigma }}^f\text{@add@centering}$$

$$\text{centering}\overline{\Sigma }{}^f\text{@add@centering}$$

$$={\Sigma }^f-C^T{\left({\Sigma }^s\right)}^{-1}C+C^T{\left({\Sigma }^s\right)}^{-1}{\Sigma }^s{\left({\Sigma }^s\right)}^{-1}C={\Sigma }^f,$$

$$A=\left[\begin{array}{cc}{A}^s& V\\ W& A^f\end{array}\right].$$

$${\mu }_n=\left[\begin{array}{c}{\mu }_n^s\\ {\mu }_n^f\end{array}\right].$$

$$X_{n,1}=(I+A\delta t)X_n+\sqrt{\delta t}B{\xi }_n,{\xi }_n\sim \mathcal{N}(0,1),$$

$${\mu }_{n,1}=(I+A\delta t){\mu }_n,{\Sigma }_{n,1}=(I+\delta tA){\Sigma }_n(I+\delta tA{)}^T+\delta tB{B}^T,$$

$$\begin{array}{rl}{\mu }_{n+1}^s& ={\mu }_n^s+\frac{\Delta t}{\delta t}\left({\mu }_{n,1}^s-{\mu }_n^s\right)\\ & ={\mu }_n^s+\frac{\Delta t}{\delta t}\left({\mu }_n^s+\delta t{A}^s{\mu }_n^s+\delta tV{\mu }_n^f-{\mu }_n^s\right)\\ & =(I^s+\Delta t{A}^s){\mu }_n^s+\Delta tV{\mu }_n^f.\end{array}$$

$$\begin{array}{rl}{\mu }_{n+1}^f& ={\mu }_{n,1}^f+C_{n,1}^T{\left({\Sigma }_{n,1}^s\right)}^{-1}\left({\mu }_{n+1}^s-{\mu }_{n,1}^s\right)\\ & =\delta tW{\mu }_n^s+(I^f+\delta t{A}^f){\mu }_n^f+C_{n,1}^T{\left({\Sigma }_{n,1}^s\right)}^{-1}\left((\Delta t-\delta t)A^s{\mu }_n^s+(\Delta t-\delta t)V{\mu }_n^f\right)\\ & =\left[\delta tW+(\Delta t-\delta t)C_{n,1}^T{\left({\Sigma }_{n,1}^s\right)}^{-1}{A}^s\right]{\mu }_n^s+\left[I^f+\delta t{A}^f+(\Delta t-\delta t)C_{n,1}^T{\left({\Sigma }_{n,1}^s\right)}^{-1}V\right]{\mu }_n^f.\end{array}$$

$$\left[\begin{array}{c}{\mu }_{n+1}^s\\ {\mu }_{n+1}^f\end{array}\right]=\left[\begin{array}{cc}{I}^s+\Delta t{A}^s& \Delta tV\\ \delta tW+(\Delta t-\delta t)C_{n,1}^T{\left({\Sigma }_{n,1}^s\right)}^{-1}{A}^s& I^f+\delta t{A}^f+(\Delta t-\delta t)C_{n,1}^T{\left({\Sigma }_{n,1}^s\right)}^{-1}V\end{array}\right]\left[\begin{array}{c}{\mu }_n^s\\ {\mu }_n^f\end{array}\right].$$

$$\underset{\delta t\to 0}{\lim }\underset{\Delta t\to \delta t}{\lim }\mathcal{D}(P_{n_{\Delta t}(T)}||P_T)=0.$$

$$\mathcal{D}(P_{n_{\Delta t}(T)}||P_T)$$

Full text

*Keywords and phrases: * micro-macro acceleration methods, stiff stochastic differential equations, entropy minimization, Kullback-Leibler divergence, convergence & stability

1 Introduction

Applications with multiple time scales arise in many domains, such as nanoscience [1], fluid dynamics [2], material science [3] and life sciences [4]. Still, the design and analysis of efficient numerical schemes for stiff stochastic differential equations (SDEs) remains challenging. On the one hand, due to a small stability domain, explicit schemes require too many time steps to reach the end of the simulation. On the other hand, while implicit methods are successful for ordinary differential equations, they fail to compute the correct invariant distribution for SDEs [5]. Therefore, stiff SDEs require new dedicated numerical multiscale methods. A lot of work has already been done over the years, and we refer to the heterogeneous multiscale method [6, 7], equation-free techniques [8, 9], and S-ROCK [10, 11] in particular, as starting points in the literature.

Recently, a micro-macro acceleration algorithm was introduced to accelerate the Monte Carlo simulation of SDEs with a time-scale separation between the fast individual stochastic paths and some slow macroscopic state variables of interest, which we define as expectations of some quantities of interest over the microscopic distributions [12]. Micro-macro acceleration connects the two levels of description of the process: individual paths of the SDE model on the microscopic level and the macroscopic state variables at the macroscopic level. The method alleviates the computational cost of a direct Monte Carlo simulation by interleaving a short bursts of Monte Carlo simulation with extrapolation of the macroscopic state variables over a larger time step. After extrapolation, the method constructs a new probability distribution by matching the last available distribution after the microscopic simulation with the extrapolated macroscopic state variables. Matching minimally perturbs the last available distribution after the microscopic simulation to make it consistent with the extrapolated state variables. Thus, in this approach, matching is inherently an optimization problem. There are many ways of measuring the difference between probability distributions. Following [13], we choose relative entropy or Kullback-Leibler divergence in the present paper, based on information theoretic considerations; for some other strategies see [12].

A few fundamental properties of micro-macro acceleration have already been investigated. First, in [13], convergence was studied for general SDEs and for any time-scale separation. The analysis requires fixing an infinite hierarchy of macroscopic state variables that are used for extrapolation. It is then shown that convergence not only depends on taking the micro and extrapolation time steps to zero, but also on extrapolating an increasing number of macroscopic state variables as the extrapolation time step decreases. Second, in [14], the asymptotic numerical stability of the micro-macro acceleration method was studied on a linear system of SDEs with Gaussian initial conditions. The stability criterion that was employed checks if the distributions obtained from the micro-macro acceleration method reach the equilibrium Gaussian distribution of the underlying microscopic integrator as the number of micro-macro time-steps tends to infinity. The analysis reveals that, when the slow and fast components of the system are decoupled, the maximal extrapolation time step is independent of the time-scale separation.

Convergence and stability are concerned with two limiting situations: convergence studies the method’s behaviour at a fixed moment in time as the extrapolation time step tends to zero, whereas stability studies the method’s behaviour for large time steps and long time horizons. Once both convergence and stability have been established, one can look at the appropriate selection of the extrapolation time step and the number of macroscopic state variables for accuracy and efficiency. In [15], we recently investigated the accuracy and efficiency of the micro-macro acceleration on slow-fast systems, showing numerically that micro-macro acceleration can simultaneously take larger time steps than the microscopic time integrator, while obtaining a smaller error than approximate macroscopic models for the slow component of the system.

In this work, we expand the convergence and stability study of the micro-macro acceleration method on linear slow-fast SDEs with additive noise. For these equations, we prove convergence to the microscopic dynamics for Gaussian initial conditions when only extrapolating the mean of the slow component. To this end, we look at the propagation of the mean and variance of the full microscopic state throughout the micro-macro acceleration scheme, only extrapolating the mean of the slow component of the SDE. We prove that, when the extrapolation time step goes to zero, these two first moments converge to the corresponding ones produced by the underlying Euler-Maruyama scheme. Dragging the microscopic time step to zero, we further obtain convergence to the exact dynamics of the slow-fast SDE. In contrast to the convergence result from [13], this analysis does not require any hierarchy of macroscopic state variables; both the mean of the fast component and the variance of the SDE are never extrapolated.

We also present a stability result for non-Gaussian initial laws that complements the previous findings in the Gaussian framework [14]. For this part, we consider a class of (non-Gaussian) initial conditions having Gaussian tails. Inspired by the case of Gaussian initial conditions analyzed in [14], where the explicit formulas are available, we prove that the stability of the micro-macro acceleration method hinges on the stability of the mean obtained from the matching procedure. In this case, however, due to the non-Gaussianity of distributions, we have to look at the propagation of all higher moments throughout the method.

These findings illustrate a certain robustness of the matching procedure: solely using the first moment, matching reconstructs the remaining higher moments, so that we converge to the microscopic scheme as extrapolation vanishes, and recover its invariant distribution as the number of fixed extrapolations grows to infinity.

Slow-fast linear SDEs with additive noise

The objects of interest in this manuscript are linear SDEs with additive noise, or Ohrnstein-Uhlenbeck processes, of the form

[TABLE]

with a square drift matrix $A\in\mathbb{R}^{d\times d}$, a rectangular diffusion matrix $B\in\mathbb{R}^{d\times m}$ and the Wiener process $W_{t}\in\mathbb{R}^{m}$. There are a few reasons why linear SDEs with additive noise are useful to study. First, the dynamics of linear systems is well understood, and we can derive stronger convergence results of micro-macro acceleration on such systems. Second, Ohrnstein-Uhlenbeck processes have an invariant distribution. We can then investigate whether micro-macro acceleration converges in distribution to the correct equilibrium distribution, as was done in [14] as a function of the time-scale separation and the extrapolation step size. Third, linear systems are popular in the context of ordinary differential equations to determine the stability of deterministic methods by an eigenvalue analysis of the drift matrix $A$. Although the concept of linearization is ambiguously defined in the stochastic case, linear SDEs with additive noise are useful to study in their own right.

In this manuscript, we are concerned with linear SDEs with a time-scale separation between some slow and some fast variables. A spectral gap in the drift matrix $A$ is a good indication of a time-scale separation present in the linear system. The larger the gap, the more stiff the linear system (1) becomes. In the context of micro-macro acceleration, we are mainly interested in the evolution of some moments of the slow components of (1). Using the spectral decomposition theorem, we introduce the orthogonal projections $\Pi^{s}\colon\mathbb{R}^{d}\to\mathbb{R}^{d_{s}}$ and $\Pi^{f}\colon\mathbb{R}^{d}\to\mathbb{R}^{d-d_{s}}$ that map the full state space onto the ‘slow’ $\mathbb{R}^{d_{s}}$ and ‘fast’ $\mathbb{R}^{d-d_{s}}$ state spaces, respectively, which correspond to the gap in the spectrum of $A$. Such a procedure is also called ‘coarse-graining’ [14]. The decomposition is such that we can express the full state space and the drift matrix as

[TABLE]

We aim at approximating the moments of the projected process $\Pi^{s}X_{t}$ as well as possible to compute the exact evolution of these moments with a reasonable accuracy. For the remainder of the manuscript, a superscript ‘s’ denote the slow components and a superscript ‘f’ the fast.

Outline of the paper

The paper is organized as follows: in Section 2, we introduce the micro-macro acceleration algorithm, specifically in our context of linear SDEs with additive noise. Section 3 contains the proof of convergence of micro-macro acceleration to the complete microscopic dynamics with only slow mean extrapolation. Section 4 investigates the stability of the method when the initial condition has Gaussian tails, when only extrapolating slow mean. In Section 5, we illustrate the theoretical results on convergence and stability with numerical examples. We consider a system of linear SDEs with additive noise, with an extra periodic force on the slow component. Section 6 contains a concluding discussion.

2 The micro-macro acceleration algorithm

One cycle of the micro-macro acceleration consists of four parts: (i) a microscopic simulation of the (stiff) stochastic differential equation over a short time interval, discretized with small time steps $\delta t$; (ii) restriction or computing an estimate of the macroscopic state variables based on the microscopic ensembles from (i); (iii) extrapolation of the restricted macroscopic state variables over a larger time step $\Delta t\gg\delta t$; (iv) matching the extrapolated state variables onto a probability distribution that perturbs the final distribution from (i) minimally.

Matching is the hardest step of the micro-macro acceleration algorithm. During matching, we build a new probability distribution that is consistent with the extrapolated state variables. This problem is often ill-posed, since there can be many probability distributions consistent with a given set of macroscopic state variables. Therefore, it was proposed in [12] to find the distribution that minimizes the divergence with respect to a prior distribution $P$. Such a prior is naturally available as the final distribution from the simulation step of the micro-macro acceleration method. In this work, we use matching introduced in [12, 13] and based on minimizing the Kullback-Leibler divergence (also called relative entropy)

[TABLE]

over all distributions $Q$ that are consistent with the extrapolated states.

In this section, we first derive some explicit formulas for this matching procedure for linear slow fast SDEs (Section 2.1), after which we present the complete micro-macro acceleration method in full mathematical detail (Section 2.2).

2.1 Matching with the slow mean

In this paper, we are particularly interested in the matching procedure that reconstructs a full microscopic distribution based only on the slow mean. For more general matching operators, see [12, 13]. In this case, denoting by $\bar{\mu}^{s}$ the mean of the slow component, the matching reads

[TABLE]

The distribution $\bar{Q}$ solving (2) is always absolutely continuous with respect to the prior $P$ and its density has exponential shape given by

[TABLE]

where the normalization constant (log-partition function) is $A(\lambda^{s},P)=\ln\mathbb{E}_{P}\left[\exp\left(\lambda^{s}\cdot\Pi^{s}\right)\right]$. The optimal Lagrange multipliers $\bar{\lambda}^{s}\in\mathbb{R}^{d_{s}}$ are unique and fulfill [13]

[TABLE]

In particular, when the prior $P$ has density $\pi_{P}$ with respect to the Lebesgue measure on $\mathbb{R}^{d}$, so does $\bar{Q}$ and its density reads

[TABLE]

Moreover, when the prior $P$ is Gaussian, there is a closed expression for $\bar{Q}$, as will become clear in Lemma 1 (Section 3.1). When $P$ is not Gaussian, we need to resort to numerical methods to solve (3) for the Lagrange multipliers, see, e.g., [12].

The following result connects matching of the full microscopic distribution $P$ with a given slow mean $\bar{\mu}^{s}$ to the corresponding matching procedure using the slow marginal $P^{s}$ of the microscopic distribution as prior.

Proposition 1**.**

Let $\bar{Q}=\mathcal{M}(\bar{\mu}^{s},P)$ be the solution to (2) and $\bar{Q}^{s}=\mathcal{M}(\overline{\mu}^{s},P^{s})$ be the solution to the matching of slow prior marginal $P^{s}$. Then, the matching densities satisfy

[TABLE]

and, in particular, all slow observables of $\bar{Q}$ equal the corresponding observables of $\bar{Q}^{s}$. Moreover, for any function $f$ on $\mathbb{R}^{d_{f}}$, the fast observable of $\bar{Q}$ generated by $f$ is given as

[TABLE]

Proof.

To see the first identity, consider the log-partition function $A(\lambda,P)$. Employing the decomposition of $P$ into its marginal and conditional [16, Thm. 10.2.1], we compute

[TABLE]

where $A(\lambda,P^{s})$ is the log-partition function of the marginal prior $P^{s}$. Since both log-partition functions agree, the vector of Lagrange multipliers $\bar{\lambda}\in\mathbb{R}^{d_{s}}$ of $\bar{Q}$ corresponds exactly to the one of $\bar{Q}^{s}$. Therefore, we can write \mathchoice{\frac{\operatorname{d}\!\bar{Q}}{\operatorname{d}\!P}}{\operatorname{d}\!\bar{Q}/\!\operatorname{d}\!P}{\operatorname{d}\!\bar{Q}/\!\operatorname{d}\!P}{\operatorname{d}\!\bar{Q}/\!\operatorname{d}\!P}(y,z)=\exp\big{(}\overline{\lambda}\mathbin{\mathbf{\cdot}}y-A(\bar{\lambda},P^{s})\big{)}, which proves the first identity.

To unwrap the fast mean of $\bar{Q}$, let us first write

[TABLE]

Using the identity between the densities and [16, Thm. 10.2.1] once more, we arrive at

[TABLE]

which concludes the proof of (4). ∎

2.2 The complete micro-macro acceleration method

In this section, we describe the four steps of the micro-macro acceleration algorithm in detail. We first present the time discretization of the linear SDE (1). We further introduce the restriction operator together with linear extrapolation of the macroscopic state variables. Finally, we use the matching operators as discussed in the previous section.

Let $P_{n}$ be the probability distribution at time $t_{n}=n\Delta t$. The micro-macro acceleration algorithm advances the distribution $P_{n}$ to a distribution $P_{n+1}$ at time $t_{n+1}=t_{n}+\Delta t$ in four stages:

Step 1: Microscopic time integration

In the first step, we perform a simulation of (1) over a time window of size $\delta\tau$. The computational cost of time propagation is usually high and we choose $\delta\tau$ to be of the order of the stiffest part of (1). In practice, we usually take $K$ time steps of size $\delta t$, such that $\delta\tau=K\delta t$. In this text, we use the Euler-Maruyama scheme to discretize (1), reading

[TABLE]

for $k=1,\dots,K$, and where $\delta W_{n,k-1}$ are Brownian increments. The random variables $X_{n,k},\ k=0,\dots,K$ have probability distributions $P_{n,k}$ and we denote the initial distribution as $P_{n,0}=P_{n}$.

Step 2: Restriction

Second, to transition from the full microscopic description to the macroscopic state variables, we compute the mean of the slow component of the process (1), reading $\mathbb{E}_{P}\left[\Pi^{s}\right]$. We restrict the slow mean at every microscopic time step, generating a sequence of values of the slow mean:

[TABLE]

Step 3: Extrapolation

In the third step, we perform time integration on the macroscopic level over a time interval of size $\Delta t$. Given the slow means $\mu^{s}_{n,k},\ k=0,\dots,K$ at times $t_{n}+k\delta t$ from the previous step, we compute the slow mean $\mu^{s}_{n+1}$ at time $t_{n+1}=t_{n}+\Delta t$ by linear extrapolation

[TABLE]

Note that we only use the state variables at time $t_{n}$ and $t_{n}+K\delta t$.

Step 4: Matching

Finally, we construct a new probability distribution that is consistent with $\mu^{s}_{n+1}$. To this end, we employ the matching operator (2) and define

[TABLE]

to obtain a new microscopic distribution $P_{n+1}$ at time $t_{n+1}=t_{n}+\Delta t$. The prior distribution $P_{n,K}$ is the final distribution computed during Step 1.

3 Convergence of the micro-macro acceleration method with slow mean extrapolation

In this section, we prove that the micro-macro acceleration method of Section 2 converges to the exact dynamics of the linear SDE (1), when only extrapolating the mean of the slow process and when the initial condition is Gaussian (Theorem 1 in Section 3.2). Before proceeding to the proof, we need an intermediate result, that explicitly describes the evolution of the mean and variance of the full microscopic system under the micro-macro acceleration method. This intermediate result is the subject of Section 3.1. Theorem 1 differs from the main convergence result in [13], as the latter requires a hierarchy of macroscopic state variables to form a complete description of the density it represents. The slow mean by itself never forms such a complete description of the underlying density. An extension of Theorem 1 to non-linear SDEs or non-Gaussian initial conditions is highly non-trivial.

3.1 An iterative formula for slow mean-only extrapolation

The proof of the convergence result in Theorem 1 relies on an iterative formula that describes how the complete mean and variance propagate through one step of the micro-macro acceleration scheme. As the micro-macro method preserves the Gaussianity of the initial condition, which we show below, the knowledge of the mean and variance suffices to control the distribution throughout the whole simulation. The derivation here assumes only one Euler-Maruyama inner step of size $\delta t$ for simplicity, but can easily be extended to $K$ inner steps.

We start with a lemma, proven also in [14] but by different means, that gives the matched distribution when the prior is Gaussian and we only match with the slow mean.

Lemma 1**.**

Suppose $P$ is the prior Gaussian distribution with mean $\mu$ and covariance matrix $\Sigma$,

[TABLE]

The distribution $\bar{Q}=\mathcal{M}(\bar{\mu}^{s},P)$, which solves (2), is also Gaussian with the same variance and mean $\bar{\mu}=[\bar{\mu}^{s},\bar{\mu}^{f}]^{T}$ where $\bar{\mu}^{f}=\mu^{f}+C^{T}\left(\Sigma^{s}\right)^{-1}(\bar{\mu}^{s}-\mu^{s})$.

Proof.

Let $\bar{Q}^{s}=\mathcal{M}(\bar{\mu}^{s},P^{s})$, where $P^{s}$ is the slow marginal of $P$. Since $P^{s}=\mathcal{N}(\mu^{s},\Sigma^{s})$, by the Gaussianity of $P$, a standard result for Kullback-Leibler minimization states that $\bar{Q}^{s}$ is also Gaussian with the mean $\bar{\mu}^{s}$ and the same variance $\Sigma^{s}$. That $\bar{Q}$ is Gaussian and its slow variance equals $\Sigma^{s}$ follows directly from the expression connecting the Radon-Nikodym derivatives in Proposition 1.

To compute the fast mean $\bar{\mu}^{f}$ and variance $\bar{\Sigma}^{f}$ of $\bar{Q}$, we use the second part of Proposition 1. Focusing on $\bar{\mu}^{f}$ first, employing formula (4) with $f(z)=z$ and a well-known expression for the conditional mean for Gaussian distributions, we can write the matched fast mean as

[TABLE]

Similarly, we can express the fast matched variance by choosing $f(z)=(z-\bar{\mu}^{f})(z-\bar{\mu}^{f})^{T}$ in (4)

[TABLE]

By adding and subtracting the fast conditional mean $\mu^{f|s}(y)=\mu^{f}+C^{T}\left(\Sigma^{s}\right)^{-1}(y-\mu^{s})$ to each $Z-\bar{\mu}^{f}$, we get

[TABLE]

The first summand under the integral represents the fast conditional variance of $P$ and, since $P$ is Gaussian, it is independent of $y$ and equals $\Sigma^{f}-C^{T}\left(\Sigma^{s}\right)^{-1}C$. The second summand can be expanded using the expressions for $\mu^{f|s}(y)$ and $\bar{\mu}^{f}$. Writing it out, we obtain

[TABLE]

where we used the fact that the slow variance of $\bar{Q}^{s}$ equals $\Sigma^{s}$. ∎

With Lemma 1, we are armed to obtain a closed expression for the time-discrete evolution of the mean of the slow-fast SDE (1), as generated by the micro-macro acceleration method of Section 2. Suppose that the drift matrix $A$ is given in block form

[TABLE]

and that at time $t_{n}=n\Delta t$ the distribution is Gaussian with mean $\mu_{n}$, and covariance $\Sigma_{n}$ with

[TABLE]

We can then write each of the four steps of the algorithm in explicit form. First, we consider the microscopic simulation step. During one Euler-Maruyama step, a Monte Carlo particle $X_{n}$ is propagated as

[TABLE]

and the distribution of $X_{n,1}$ is also Gaussian. Next, we perform the restriction step. By taking expectations, the mean $\mu_{n,1}$ and the covariance $\Sigma_{n,1}$ after the Euler-Maruyama step read

[TABLE]

after which we extrapolate the slow mean as

[TABLE]

Now that we have the extrapolated slow mean, we can use Lemma 1 to explicitly obtain the result of matching. According to Lemma 1, the matched distribution is also Gaussian when only extrapolating the slow mean. Furthermore, the covariance matrix is not affected by matching, i.e., $\Sigma_{n+1}=\Sigma_{n,1}$, and the fast mean is given by

[TABLE]

Bundling the propagation of the slow (8) and fast matched mean (9) in one vector $\mu_{n+1}$ gives

[TABLE]

To conclude, the time-discrete evolution of the mean of the slow-fast SDE (1), as generated by the micro-macro acceleration method is given by the time-dependent linear system (10), with initial condition equal to the mean of the initial distribution of (1).

3.2 Convergence theorem

All elements are now in place to prove convergence of the micro-macro acceleration method of Section 2 that only extrapolates the slow mean of the process. The proof makes use of the iterative formula above, and holds for general linear SDEs with additive noise.

Theorem 1**.**

Given a linear SDE with a Gaussian initial distribution, consider the micro-macro acceleration algorithm with relative-entropy matching and slow-mean extrapolation. Also, fix an end time $T>0$. Denote by $P_{T}$ the exact distribution of the linear SDE at time $T$, and by $P_{n_{\Delta t}(T)}$ the distribution obtained using $n_{\Delta t}(T)$ steps of the micro-macro acceleration scheme of Section 2 with $K=1$, where $n_{\Delta t}(T)=\lfloor T/\Delta t\rfloor$. Then,

[TABLE]

As a consequence, the distributions $P_{n_{\Delta t}(T)}$ obtained by micro-macro acceleration at time $T$ converge in total variation to the exact microscopic distribution $P_{T}$ in the same limits.

Proof.

Since the initial condition is Gaussian, all intermediate distributions of the exact solution, the Euler-Maruyama method, and micro-macro acceleration are Gaussian too. By a standard expression for the Kullback-Leibler divergence between two Gaussian distributions [17], the divergence in (11) becomes

[TABLE]

where $\mu_{T}$ and $\Sigma_{T}$ are the mean and variance of $X_{T}$, and $\mu_{n_{\Delta t}(T)}$ and $\Sigma_{n_{\Delta t}(T)}$ are the mean and variance of $X_{n_{\Delta t}(T)}$.

First, we fix $\delta t\leq\Delta t$ and let $\Delta t$ decrease to $\delta t$. If we perform back-substitution in equation (10) to write the mean at time $T$ as a function of the initial mean vector, we obtain a product of $n_{\Delta t}(T)$ different matrices. The number of matrices increases to $n_{\delta t}(T)$ as $\Delta t$ decreases to $\delta t$, but there always remain a finite number of matrices because $\delta t>0$. The contribution of the largest off-diagonal term in (10) also reduces to zero and as a result, the mean vector $\mu_{n_{\Delta t}(T)}$ approaches the respective mean $\mu_{n_{\delta t}(T)}$ of the Euler-Maruyama scheme. We obtain

[TABLE]

where $\mu_{0}^{s}$ and $\mu_{0}^{f}$ are the slow and fast mean respectively of the initial condition. Similarly, since the variance stays constant during matching, $\Sigma_{n_{\Delta t}(T)}$ converges to $\Sigma_{n_{\delta t}(T)}$, the variance of the Euler-Maruyama scheme. Hence, the limit (11) reduces to

[TABLE]

which is the relative entropy between the Euler-Maruyama scheme and the exact solution at time $T$. This expression converges to zero because, as $\delta t$ decreases to zero, the mean and variance of the Euler-Maruyama method converge to their respective values of the exact solution at time $T$. Finally, by Pinsker’s inequality, $P_{n_{\Delta t}(T)}$ converges to $P_{T}$ in total variation. ∎

Theorem 1 might be surprising, since it does not require the number of macroscopic state variables to increase to infinity as the extrapolation time step $\Delta t$ decreases; using only the slow mean as a macroscopic state variable is sufficient. However, the result only holds for linear SDEs with Gaussian initial conditions, as the proof relies heavily on iteration (10). At the moment, no proof exists on convergence for non-linear SDE with any (fixed) finite number of macroscopic state variables.

4 Stability of micro-macro acceleration with initial condition with Gaussian tails

In this section, we study the stability of micro-macro acceleration when applied to (1), i.e., the convergence of the laws it generates, in the limit as the number of extrapolation with fixed step size $\Delta t$ goes to infinity, to the invariant distribution of the underlying Euler-Maruyama scheme. When $\delta t$ denotes the microscopic step, the invariant distribution of the Euler-Maruyama scheme is the zero-mean normal distribution $\mathcal{N}_{0,V_{\infty}^{\delta t}}$, with variance

[TABLE]

as can be seen by repeatedly applying the recursion for the mean and variance in (7).

The stability question was analyzed in [14] in the Gaussian setting, where the distributions that are generated by the micro-macro acceleration can be computed explicitly. The main result there focuses on a simple diagonal case and reads:

Theorem 2**.**

When applying the micro-macro acceleration method to the linear SDE (1) with block-diagonal drift matrix

[TABLE]

the mean $\mu_{n}$ and the covariance matrix $\Sigma_{n}$ of the resulting Gaussian law at the $n$th step satisfy

[TABLE]

whenever

[TABLE]

Here, $\rho(\cdot)$ denotes the spectral radius of a matrix. Condition (13) is necessary to stabilize the extrapolation of the slow mean $\mu_{n}^{s}$, by bounding the values of $\Delta t$, and to stabilize the Euler-Maryuama stage, by bounding the values of $\delta t$. For the proof of Theorem 2, it suffices to establish the proper asymptotic behaviour of the means $\mu_{n}$ and variances $\Sigma_{n}$ of Gaussian distributions generated by the micro-macro acceleration scheme.

In this section, we go beyond the Gaussian case and work within the larger class of probability measures that have Gaussian tails. As a consequence, we do not have explicit formulas for the distributions generated by the scheme and we have to control all moments to show stability. Rather than trying to obtain stability bounds on $\Delta t$ directly, we concentrate on showing that convergence of $\mu^{n}$ to [math], i.e., asymptotic stability of the mean, already yields convergence of the distributions to $\mathcal{N}_{0,V_{\infty}^{\delta t}}$. The relation between these two notions was illuminated in the Gaussian setting in [14].

The main stability result of this paper, Theorem 3 in Section 4.2, gives weak convergence of distributions produced by the micro-macro acceleration method to $\mathcal{N}_{0,V_{\infty}^{\delta t}}$, as the number of extrapolation steps goes to infinity. To prove it, we explore the properties of cumulant generating functions (CGFs) to produce a recursion formula for the laws generated by the method (Section 4.1).

4.1 Micro-macro step in terms of the cumulant generating function

Let us first define the cumulant generating function as follows.

Definition 1**.**

For any probability distribution $P\in\mathscr{P}^{d}$, we define the cumulant generating function of $P$

[TABLE]

where $\Pi$ is the identity on $\mathbb{R}^{d}$, and the effective domain reads $\Theta_{P}=\{\theta\in\mathbb{R}^{d}:\ \mathcal{K}_{P}(\theta)<+\infty\}$.

When $X\sim P$ we also write $\mathcal{K}_{X}$ instead of $\mathcal{K}_{P}$, and $\Theta_{X}$ instead of $\Theta_{P}$. We summarize the basic properties of CGFs in A.

Example 1**.**

If $X\sim\mathcal{N}_{\mu,\Sigma}$, then

[TABLE]

To effectively use the CGFs to describe the micro-macro acceleration procedure, we assume that for the initial random variable $X_{0}$ it holds $\Theta_{0}\doteq\Theta_{X_{0}}=\mathbb{R}^{d}$. In view of Proposition 5, the CGF of the matched distribution results from shifting and translating the CGF of the prior based on the current value of the Lagrange multipliers. Having priors with full effective domain avoids the issue of falling outside the effective domain while shifting the CGF – a clear sign that the matching is impossible.

Remark 1** (What does $\Theta_{0}=\mathbb{R}^{d}$ mean?).**

Let us fix $\theta\in\mathbb{R}^{d}$ and $r>0$. According to the Chernoff’s bound [18, p. 392], we have

[TABLE]

for all $s\geq 0$. From Proposition 4(iii) in the Appendix, applied with $l=1$ and $M=\theta^{\!\mathsf{T}}$, we have $\mathcal{K}_{\theta\mathbin{\mathbf{\cdot}}X_{0}}(s)=\mathcal{K}_{X_{0}}(s\theta)$. Thus, taking logarithms on both sides of (14), we can equivalently write

[TABLE]

where we denote $\mathcal{K}_{0}=\mathcal{K}_{X_{0}}$. Since $\mathcal{K}_{0}(s\theta)$ is finite for all $s\geq 0$, dividing by $r$ and taking the limit gives

[TABLE]

Therefore, because $s$ can be arbitrarily large, the log-tail function of $\theta\mathbin{\mathbf{\cdot}}X_{0}$ is superlinear at $+\infty$. The same holds at $-\infty$ by repeating the argument for $\mathbb{P}(\theta\mathbin{\mathbf{\cdot}}X_{0}\leq-r)$.

To simplify the notation, we again use only one micro step for each extrapolation in the micro-macro acceleration procedure. In the Proposition below, we consider the micro-macro acceleration method as applied to the linear slow-fast SDE (1).

Proposition 2**.**

Let $n\geq 1$ and assume that for a random variable $X_{n-1}$ with cumulant generating function $\mathcal{K}_{X_{n-1}}=\mathcal{K}_{n-1}$ we have $\Theta_{n-1}=\mathbb{R}^{d}$. Then, if $X_{n}$ is obtained from the micro-macro procedure with extrapolation of the (slow) $s$-marginal mean (as described in Section 2.2, with $K=1$), its CGF $\mathcal{K}_{n}=\mathcal{K}_{X_{n}}$ has effective domain $\Theta_{n}=\mathbb{R}^{d}$ and satisfies

[TABLE]

where $\lambda^{s}_{n}$ is a vector of Lagrange multipliers.

Proof.

According to Section 2, the law of $X_{n}$ is given by $\mathcal{M}(\mu^{s}_{n},\mathrm{Law}(X^{\delta t}_{n-1,1}))$, where $\mu^{s}_{n}$ is extrapolated as in (6) and $X^{\delta t}_{n-1,1}$ is obtained from $X_{n-1}$ by one Euler-Maruyama step over $\delta t$. Therefore, for the cumulant, Proposition 5 yields

[TABLE]

where $\lambda^{s}_{n}$ is the Lagrange multiplier associated to the extrapolated marginal mean $\mu^{s}_{n}$.

In the current notation, the recursive formula (5) reads

[TABLE]

Because $X_{n-1}$ and $\delta W_{n,1}$ are independent, we obtain from Proposition 4(v) that

[TABLE]

The law of $B\delta W_{n,1}$ is $\mathcal{N}_{0,\delta tB\!B^{\!\mathsf{T}}}$ thus, according to Example 1, $\mathcal{K}_{B\delta W_{n,1}}(\theta)=\delta t/2\,\theta^{\!\mathsf{T}}\!BB^{\!\mathsf{T}}\theta$ and, using Proposition 4(iii), we obtain

[TABLE]

Combining (16) with (17) results in (15). Since the effective domain of $X_{n-1}$ was equal to $\mathbb{R}^{d}$, the right-hand side of (15) is finite for all $\theta$. This implies $\Theta_{n}=\mathbb{R}^{d}$. ∎

4.2 Convergence to the equilibrium with stable mean extrapolation

The main result of this section depends on the following assumption on the initial random variable for the micro-macro acceleration method. Recall from Proposition 4(iv) that the cumulant generating function is always analytic on the interior of its effective domain. Here and in what follows, for any constant $C\geq 0$ and functions $f,g\colon\mathbb{R}_{+}\to\mathbb{R}$ with $g$ positive in a neighbourhood of infinity, $f(r)\sim Cg(r)$ means $\lim_{r\to+\infty}f(r)/g(r)=C$.

Assumption 1**.**

The CGF $\mathcal{K}_{0}$, of an initial random variable $X_{0}$, has full effective domain (i.e., $\Theta_{0}=\mathbb{R}^{d}$) and for every $\theta\in\mathbb{R}^{d}$ it satisfies

[TABLE]

where $w_{0}\colon\mathbb{R}^{d}\to(0,+\infty)$ is continuous and homogeneous of order $2$.

Assumption 1 derives from the theory of regular variation. Let us first discuss its connection with the tails of random variables.

Remark 2** (On log-quadratic tails).**

The asymptotic relation $\mathcal{K}_{0}(r\theta)\sim w_{0}(\theta)\,r^{2}$ and Proposition 4(iii) leads to $\mathcal{K}_{\theta\mathbin{\mathbf{\cdot}}X_{0}}(r)\sim w_{0}(\theta)\,r^{2}$. By the Kasahara-Tauberian Theorem [19, Thm. 4.12.7], the last relation is equivalent to the property

[TABLE]

That is, for every $\theta\in\mathbb{R}^{d}$, the random variable $\theta\mathbin{\mathbf{\cdot}}X_{0}$ has regularly varying log-quadratic tail decay, a feature shared by all Gaussian laws. Therefore, Assumption 1 sharpens the superlinear behaviour of tail functions that resulted from assuming $\Theta_{0}=\mathbb{R}^{d}$, see Remark 1. The inclusion of derivatives in (18) is related to the notion of higher-order regular variation(compare [20] and [19, p. 44]).

In the proof of Theorem 3 below, we employ the following technical result.

Lemma 2**.**

Let $\mathcal{K}$ be continuously differentiable function on $\mathbb{R}^{d}$ that for every $\theta\in\mathbb{R}^{d}$ satisfies $\mathcal{K}(r\theta)\sim w(\theta)\,r^{2}$ and

[TABLE]

Then, for every $\theta,\theta_{0}\in\mathbb{R}^{d}$, it holds $\mathcal{K}(r\theta+\theta_{0})\sim 2\,w(\theta)r^{2}$.

Proof.

Express the ratio $\mathcal{K}(r\theta+\theta_{0})/r^{2}$ as the following sum

[TABLE]

The second summand converges to $2w(\theta)$ as $r\to+\infty$ from the asymptotic property of $\mathcal{K}$. That the first fraction disappears can be seen by applying the mean value inequality

[TABLE]

where $r^{\prime}\in[0,1]$. Since

[TABLE]

stays bounded as $r\to+\infty$, due to (19), the right-hand side of the previous inequality converges to zero in this limit. ∎

Theorem 3**.**

Suppose that the CGF $\mathcal{K}_{0}$ of an initial random variable $X_{0}$ satisfies Assumption 1, both $\mathcal{K}_{0}(\theta)$, and $\nabla_{\!}\mathcal{K}_{0}(\theta)\cdot\theta$ fulfill (19), and $\mathcal{K}_{0}$ has all derivatives of order $3$ or higher bounded. If $\rho({I}+\delta tA)<1$, and the mean $\mu_{n}$ of random variables $X_{n}$ obtained from the micro-macro acceleration method is stable, that is $\lim_{n\to+\infty}\mu_{n}=0$, the laws of $X_{n}$ converge weakly to $\mathcal{N}_{0,V_{\infty}^{\delta t}}$.

The stability of the mean $\mu_{n}$ is closely related to the stability bounds on $\Delta t$ that guarantee that the slow marginal mean $\mu_{n}^{s}$ satisfies $\lim_{n\to\infty}\mu_{n}^{s}=0$. The behaviour of fast marginal $\mu_{n}^{f}$ is influenced at each step by both extrapolation over $\Delta t$ and the matching. In the Gaussian setting, Lemma 1 provides the exact formula for $\mu_{n}^{f}$ in terms of $\mu_{n}^{s}$ and the first two moments of distributions produced by the Euler-Maryuama step. Using this formula, we can obtain explicit stability bounds on $\Delta t$ and $\delta t$, exemplified in (13). In the non-Gaussian setting, we can only say that the influence of matching on $\mu_{n}^{f}$ is encoded in the nonlinear procedure to obtain Lagrange multipliers $\lambda_{n}$. Having no explicit formulas for the multipliers we work under the assumption that the full mean is stable when using the micro-macro acceleration method.

Proof.

To establish stability, we employ the recursive relation (15). However, due to the presence of the Lagrange multipliers $\lambda^{s}_{n}$ in the argument of $\mathcal{K}_{n}$ and in the additional last term of (15), we cannot immediately pass to the limit as $n$ goes to $+\infty$. These multipliers exist for all $n$, a consequence of $\Theta_{n}=\mathbb{R}^{d}$, but we do not have any a priori estimates that would allow to control $\lambda^{s}_{n}$ as $n$ increases. Therefore, our strategy in proving the convergence is to look at the recurrences for the tails of $\mathcal{K}_{n}$, which do not contain $\lambda^{s}_{n}$ any more. Before we look at $\mathcal{K}_{n}$ itself, let us use the boundedness of its higher order derivatives to show the convergence of higher cumulants.

First note that the last term in formula (15) is of second order in $\theta$, so it disappears after differentiating this identity three times. More precisely, we have, for any $j\geq 3$,

[TABLE]

where we treat the $j$th order derivative $\operatorname{D}\!^{j}\mathcal{K}_{n}(\theta)$ at $\theta\in\mathbb{R}^{d}$ as a symmetric multilinear mapping on $(\mathbb{R}^{d})^{j}$ and denote $\operatorname{D}\!^{j}\mathcal{K}_{n}(\theta)[\theta^{\prime}]\doteq\operatorname{D}\!^{j}\mathcal{K}_{n}(\theta)[\theta^{\prime},\dotsc,\theta^{\prime}]$. By back-substituting we obtain

[TABLE]

Since $\operatorname{D}\!^{j}\mathcal{K}_{0}$ is bounded on $\mathbb{R}^{d}$ for all $j\geq 3$, we can estimate from the above relation that

[TABLE]

where $\|\operatorname{D}\!^{j}\mathcal{K}_{n}(\theta)\|_{\mathrm{mult}}\doteq\sup|\operatorname{D}\!^{j}\mathcal{K}_{n}(\theta)(\theta^{1},\dotsc,\theta^{j})|/\|\theta^{1}\|\dotsm\|\theta^{j}\|$ and the supremum is taken over all $\theta^{i}\neq 0$. Since $\rho({I}+\delta tA)<1$, this demonstrates that $\operatorname{D}\!^{j}\mathcal{K}_{n}(\theta)$ converges to zero as $n\to+\infty$, uniformly in $\theta$. In particular, $\lim_{n\to+\infty}\operatorname{D}\!^{j}\mathcal{K}_{n}(0)=0$ for all $j\geq 3$.

Let us now return to (15). We show by induction that for all $\theta\in\mathbb{R}^{d}$ the function $r\mapsto\mathcal{K}_{n}(r\theta)$ is asymptotically quadratic in $r$ and find the recursion for the corresponding constants $w_{n}(\theta)$. Suppose that $\mathcal{K}_{n-1}(r\theta)\sim w_{n-1}(\theta)r^{2}$ with some function $w_{n-1}$. Note that Assumption 1 guarantees that this holds for $\mathcal{K}_{0}$. The recursive relation (15) gives

[TABLE]

Note that simultaneously, by differentiating (15), the gradient of $\mathcal{K}_{n}$ satisfies

[TABLE]

thus property (19) propagates from $\mathcal{K}_{0}$ throughout all $\mathcal{K}_{n}$. Therefore, we can use Lemma 2, together with the inductive assumption, to conclude that the limit as $r$ goes to $+\infty$ on the right hand side exists and is equal to 2w_{n-1}\big{(}(I+\delta tA^{\!\mathsf{T}})\theta\big{)}+\delta t\theta^{\!\mathsf{T}}BB^{\!\mathsf{T}}\theta. Denoting by $2w_{n}(\theta)$ the limit of the left hand side, we have the recursion w_{n}(\theta)=w_{n-1}\big{(}(I+\delta tA^{\!\mathsf{T}})\theta\big{)}+(\delta t/2)\theta^{\!\mathsf{T}}\!BB^{\!\mathsf{T}}\!\theta. By back-substituting, we obtain

[TABLE]

Note also that, by a similar reasoning as for $\mathcal{K}_{n}$, but using (21) this time, we can inductively establish the relation $\operatorname{d}\!/\!\operatorname{d}\!{r}\mathcal{K}_{n}(r\theta)\sim 2w_{n}(\theta)r$. In particular, after differentiating, this relation yields

[TABLE]

with $w_{n}$ satisfying (22).

Since the micro time step $\delta t$ is stable, we can now take the limit as $n\to\infty$ on the right-hand side of (22). In consequence, we know that the point-wise limit $\omega_{\infty}\doteq\lim_{n\to\infty}\omega_{n}$ exists and, by (12), results in the following expression for the limiting function

[TABLE]

Having established the limiting behaviour of higher cumulants $\operatorname{D}\!^{j}\mathcal{K}_{n}(0)$, for $j\geq 3$, we will now use identity (24) to demonstrate that the second cumulants $\nabla^{2}_{\!}\mathcal{K}_{n}(0)$ converge to $V_{\infty}^{\delta t}$. Via the fundamental theorem of calculus we have

[TABLE]

Dividing by $r$ and changing variable of integration $t\mapsto t/r$ gives

[TABLE]

by the asymptotic property (23) of $r\mapsto\nabla_{\!}\mathcal{K}_{n}(r\theta)\cdot\theta$. Denoting \mathop{}\mkern-3.0mu\mathchoice{\hbox to0.0pt{\displaystyle\vbox{\hbox to7.00002pt{\hss$\textstyle{-}$\hss}}\hss}}{\hbox to0.0pt{\textstyle\vbox{\hbox to7.00002pt{\hss$\scriptstyle{-}$\hss}}\hss}}{\hbox to0.0pt{\scriptstyle\vbox{\hbox to5.0pt{\hss$\scriptscriptstyle{-}$\hss}}\hss}}{\hbox to0.0pt{\scriptscriptstyle\vbox{\hbox to5.0pt{\hss$\scriptscriptstyle{-}$\hss}}\hss}}\mkern-3.0mu\int_{0}^{\infty}\doteq\lim_{r\to+\infty}\frac{1}{r}\int_{0}^{r}, we have from (24)

[TABLE]

Because the Taylor expansion of $t\mapsto\mathcal{K}_{n}(t\theta)$ around $t=0$ gives $\nabla^{2}_{\!}\mathcal{K}_{n}(t\theta)=\nabla^{2}_{\!}\mathcal{K}_{n}(0)+\operatorname{D}\!^{3}\mathcal{K}_{n}(\theta^{\prime}_{t})$, with some $\theta^{\prime}_{t}$, and $\operatorname{D}\!^{3}\mathcal{K}_{n}(\theta^{\prime}_{t})$ converges to zero as $n$ goes to $+\infty$, uniformly in $\theta^{\prime}_{t}$, we get

[TABLE]

for any $\theta\in\mathbb{R}^{d}$. This proves the convergence of the second cumulants.

In conclusion, since $\nabla_{\!}\mathcal{K}_{n}(0)=\mu_{n}$ converges to [math], by assumption, $\nabla^{2}_{\!}\mathcal{K}_{n}(0)$ converges to $V_{\infty}^{\delta t}$, and all higher order derivatives to [math], we obtain the following limiting sequence of cumulants as $n$ goes to $+\infty$:

[TABLE]

This sequence uniquely determines the normal distribution with mean zero and covariance matrix $V_{\infty}^{\delta t}$. Thus, from the Frechét-Shohat theorem [21, p. 307], the laws of $X_{n}$ converge weakly to $\mathcal{N}_{0,V_{\infty}^{\delta t}}$. ∎

5 Numerical illustration: a periodically driven linear slow-fast system

In this section, we numerically illustrate the convergence and stability results of the previous sections on an academic example. As in [5, 15], we define a linear slow-fast SDE with additive noise, in which we add a periodic forcing to the slow component

[TABLE]

The parameter $\varepsilon$ is the time-scale separation between $X$ and $Y$ and controls the stiffness in the system. Introducing a periodic forcing allows to easily measure errors between the exact and numerical solution by computing the $L^{2}$ difference between both curves over one period. At the same time, the driving force does not alter stability and convergence properties of micro-macro acceleration.

We first look at some numerical convergence tests (Section 5.1), before moving on to computational experiments on stability (Section 5.2).

5.1 Convergence to the microscopic time integrator

By taking expectations of (26), one can show that the means $\mu_{X}$ (resp. $\mu_{Y}$) of the slow (resp. fast) component of the exact solution are given by

[TABLE]

where $M=\begin{pmatrix}-2&-2\\ \ 10&-10\end{pmatrix}$ and $\begin{pmatrix}\mu_{X_{0}}&\mu_{Y_{0}}\end{pmatrix}^{T}$ is the mean of the initial condition of (26). The constants $A,B,C$ and $D$ are the solution of the linear system

[TABLE]

To illustrate convergence (Theorem 1), we compute the error in the slow means obtained micro-macro acceleration against both the exact solution (27) and the numerical result obtained by the Euler-Maruyama integrator. We are mostly interested in the error between micro-macro acceleration and the Euler-Maruyama method because we want to understand the effect of extrapolation. As parameters we take $\varepsilon=0.5$ and $\varepsilon=0.05$ for the time-scale separation and we perform computations to an end time $T=6$. We choose a small time step $\delta t=\varepsilon/20$ for the microscopic time integrator and many values for the extrapolation time step $\Delta t$. The error is computed as the $L_{2}$ norm of the difference between two curves and averaged over 10 independent runs. The convergence results are depicted on Figure 1.

Figure 1 shows that the micro-macro acceleration error lowers as $\Delta t$ decreases to $\delta t$, as proven in Theorem 1. For a large $\varepsilon=0.5$ the error decreases linearly, while for a small $\varepsilon=0.05$ the decrease is quadratically. The order of convergence of micro-macro acceleration is not well-understood yet. In Figure 1 we also see that the error computed against the analytic solution and the Euler-Maruyama method is almost the same for small $\varepsilon=0.05$ since the microscopic integrator is very accurate. For larger $\varepsilon$, the error between micro-macro acceleration and the Euler-Maruyama method is smaller than with the exact solution because the Euler-Maruyama method makes a non-negligible error.

5.2 Stability with initial conditions having Gaussian tails

In the next experiment, we also consider the periodically driven linear system, but look at large extrapolation steps to study stability properties. The derivation in Section 4 does not give a threshold on the extrapolation step above which micro-macro acceleration becomes unstable and below which the algorithm is stable. To determine that a simulation was unstable, we will therefore rely on an alternative strategy, also proposed in [14] for Gaussian initial conditions. In [14], it was shown that instability of the micro-macro acceleration technique unavoidably leads to so-called matching failures, even before the solution blows up to infinity. A matching failure occurs when there exists no probability distribution that is consistent with the given macroscopic state variables $m$. In other words, the pair $(m,\mu)$ does not lie in the domain of the matching operator $\mathcal{M}$ for any prior distribution $\mu$, unless $\mu$ is consistent with $m$.

In practice, we employ a Newton-Raphson method to compute the Lagrange multipliers in (3) and detect a matching failure when the iterative method does not converge. Specifically, we must solve

[TABLE]

with a Newton-Raphson procedure, where we compute the integral using a Monte-Carlo representation of the prior distribution $P$. When the Newton-Raphson solver fails to reach the extrapolated slow mean $\bar{\mu}^{s}$ within an absolute tolerance of $10^{-11}$ in 50 steps, we mark a matching failure.

For the numerical experiment, we simulate the periodically driven linear system (26) with $\epsilon=1$ for different pairs of step sizes $(\delta t,\Delta t)$ up to the end time $T=100$. When a matching failure occurs, we mark the pair of parameters as unstable, and stable otherwise. The initial condition is Gaussian with mean zero and unit variance, which fits in the framework of Section 4. The number of Monte Carlo replicas is $N=10^{5}$ and we use $K=1$ step in the microscopic time integrator. The numerical results are summarized in Figure 2.

Since $M$ is neither diagonal nor lower-triangular, there exist no stability bounds on the extrapolation step $\Delta t$ yet [14]. As a proxy, we first compute the deterministic stability bound, as if there were no Brownian motion in (26). The eigenvalues of $M$ are $-6\pm 2i$, implying the maximal deterministic time step is $\Delta t_{\text{max}}=0.3$. The stability domain of micro-macro acceleration is V-shaped and for every $\delta t$, the maximal extrapolation step before instability is always greater than $0.3$. Micro-macro acceleration thus has good stability properties.

6 Conclusion

We presented a micro-macro acceleration scheme, based on a combination of microscopic simulation and extrapolation of some macroscopic quantities of interest. We demonstrated that using only the slow mean during extrapolation results in a convergent and stable algorithm. The proofs hold for linear stochastic differential equations with additive noise. We complemented the analysis with numerical results that indicate that the error of micro-macro acceleration decreases to zero when both the extrapolation and microscopic step size decrease to zero. For stability, we investigated for which pairs microscopic and extrapolation time steps micro-macro acceleration is stable, and compared the numerical results to the deterministic stability bounds in the case without Brownian motion. Empirically, the stability domain is V-shaped, and for every value of the microscopic time step, the maximal extrapolation step is above the deterministic stability bound.

Acknowledgements

The authors thank Kristian Debrabant for proofreading the material in Section 4. The authors acknowledge the support of the Research Council of the University of Leuven through grant ‘PDEOPT’ and of the Research Foundation – Flanders (FWO – Vlaanderen) under grant G.A003.13.

Appendix A Properties of cumulant generating functions

First, we mention the convexity of $\mathcal{K}_{P}$.

Proposition 3** ([22, Thm. A.4]).**

Let $P\in\mathscr{P}^{d}$. Then

1. (i)

The set $\Theta_{P}$ is convex. 2. (ii)

$\mathcal{K}_{P}$* is a convex function on $\Theta_{P}$, and strictly convex if and only if $P$ is not concentrated in a single point.*

When $X\sim P$ we also write $\mathcal{K}_{X}$ instead of $\mathcal{K}_{P}$, and $\Theta_{X}$ instead of $\Theta_{P}$.

Proposition 4** ([22, Thm. A.1 & A.7], [23, 31]).**

Assume that $X\sim P$, where $P\in\mathscr{P}^{d}$. Then

1. (i)

$-\infty<\mathcal{K}_{P}(\theta)\leq+\infty$* for $\theta\in\mathbb{R}^{d}$.* 2. (ii)

$\mathcal{K}_{P}(0)=0$. 3. (iii)

If $M\in\mathbb{R}^{l\times d}$ and $c\in\mathbb{R}^{l}$

[TABLE] 4. (iv)

If $0\in\operatorname{int}\Theta_{P}$, $\mathcal{K}_{P}$ is analytic on $\operatorname{int}\Theta_{P}$ with Taylor expansion around [math]

[TABLE]

*where $\operatorname{D}\!^{j}\mathcal{K}_{P}(0)$ is the *$j$th cumulant of $P$, considered as the $j$-linear mapping on $\mathbb{R}^{d}$, and $\operatorname{D}\!^{j}\mathcal{K}_{P}(0)[\theta]\doteq\operatorname{D}\!^{j}\mathcal{K}_{P}(0)[\theta,\dotsc,\theta]$. In particular, $X$ has vector mean and covariance matrix

[TABLE] 5. (v)

If $X,Y$ are independent, $\mathcal{K}_{X+Y}=\mathcal{K}_{X}+\mathcal{K}_{Y}$ and $\Theta_{X+Y}=\Theta_{X}\cap\Theta_{Y}$.

Proposition 5**.**

If $Q=\mathcal{M}(\overline{\mu}^{s},P)$ is the solution to (2) with a (slow) marginal mean $\overline{\mu}^{s}\in\mathbb{R}^{d_{s}}$ and a prior $P\in\mathscr{P}^{d}$, then

[TABLE]

where $\overline{\lambda}^{s}\in\mathbb{R}^{d_{s}}$ is a vector of Lagrange multipliers corresponding to $\overline{\mu}^{s}$ and $0^{f}=(0,\dotsc,0)\in\mathbb{R}^{d_{f}}$.

Proof.

Using formula (2) we compute

[TABLE]

It remains to note that by the definitions of log-partition function and the marginal distribution

[TABLE]

Corollary 1**.**

For any vectors $\overline{\mu},\mu\in\mathbb{R}^{d}$ and a symmetric, non-negative definite matrix $\Sigma\in\mathbb{R}^{d\times d}$, $\mathcal{M}(\overline{\mu},\mathcal{N}_{\mu,\Sigma})=\mathcal{N}_{\overline{\mu},\Sigma}$.

Proof.

According to (3), the vector of Lagrange multipliers $\overline{\lambda}$ corresponding to matching with $\overline{\mu}$ satisfies $\nabla_{\!\lambda}A(\overline{\lambda},\mathcal{N}_{\mu,\Sigma})=\overline{\mu}$. When matching with mean only, the log-partition function coincides with the cumulant generating function, that is $A(\lambda,\mathcal{N}_{\mu,\Sigma})=\mathcal{K}_{\mathcal{N}_{\mu,\Sigma}}(\lambda)$, and using formula from Example 1 we obtain

[TABLE]

Employing Proposition 5 with $Q=\mathcal{M}(\overline{\mu},\mathcal{N}_{\mu,\Sigma})$ we compute

[TABLE]