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

Przemysław Zieliński  KU Leuven, Department of Computer Science, NUMA Section,
Celestijnenlaan 200A box 2402, 3001 Leuven, Belgium  Hannes Vandecasteele  KU
Leuven, Department of Computer Science, NUMA Section, Celestijnenlaan 200A box
2402, 3001 Leuven, Belgium  Giovanni Samaey  KU Leuven, Department of Computer
Science, NUMA Section, Celestijnenlaan 200A box 2402, 3001 Leuven, Belgium

(August 27, 2024)

###### 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.

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

$dX_{t}=AX_{t}dt+BdW_{t}$ (1)

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

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

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)

$\mathcal{D}(Q||P)=\mathbb{E}_{Q}\\!\left[\ln\left(\mathchoice{\frac{\operatorname{d}\\!Q}{\operatorname{d}\\!P}}{\operatorname{d}\\!Q/\\!\operatorname{d}\\!P}{\operatorname{d}\\!Q/\\!\operatorname{d}\\!P}{\operatorname{d}\\!Q/\\!\operatorname{d}\\!P}\right)\right],$

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

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

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

$\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{)},$

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]

$\nabla_{\lambda_{s}}A(\bar{\lambda}^{s},P)=\bar{\mu}_{s}.$ (3)

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

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

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

$\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)=\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),\quad
y\in\mathbb{R}^{d_{s}},z\in\mathbb{R}^{d_{f}},$

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

$\mathbb{E}_{\bar{Q}}[f(Z)]=\int_{\mathbb{R}^{d_{s}}}\mathbb{E}_{P}[f(Z)|Y=y]\,\bar{Q}^{s}(\operatorname{d}\\!{y}).$
(4)

###### 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

$\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}$

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

$\mathbb{E}_{\bar{Q}}[f(Z)]=\int_{\mathbb{R}^{d}}f(z)\,\bar{Q}(\operatorname{d}\\!{y},\operatorname{d}\\!{z})=\int_{\mathbb{R}^{d}}f(z)\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)\,P(\operatorname{d}\\!{y},\operatorname{d}\\!{z}).$

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

$\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}$

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

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

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:

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

#### 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

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

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

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

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$,

$\mu=\begin{bmatrix}\mu^{s}\\\ \mu^{f}\end{bmatrix},\ \
\Sigma=\begin{bmatrix}\Sigma^{s}&C\\\ C^{T}&\Sigma^{f}\end{bmatrix}.$

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

$\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}).$

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

$\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).$

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

$\displaystyle\centering\bar{\Sigma}^{f}\@add@centering$
$\displaystyle=\int_{\mathbb{R}^{d_{s}}}\big{\\{}\mathbb{E}_{P}[(Z-\mu^{f|s}(y))(Z-\mu^{f|s}(y))^{T}|Y=y]+(\mu^{f|s}(y)-\bar{\mu}^{f})(\mu^{f|s}(y)-\bar{\mu}^{f})^{T}\big{\\}}\,\bar{Q}^{s}(dy).$

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

$\displaystyle\centering\mkern
1.5mu\overline{\mkern-1.5mu\Sigma\mkern-1.5mu}\mkern 1.5mu^{f}\@add@centering$
$\displaystyle=\Sigma^{f}-C^{T}\left(\Sigma^{s}\right)^{-1}C+\int_{\mathbb{R}^{d_{s}}}C^{T}\left(\Sigma^{s}\right)^{-1}(y-\bar{\mu}^{s})(y-\bar{\mu}^{s})^{T}\left(\Sigma^{s}\right)^{-1}C\,\bar{Q}^{s}(dy)$
$\displaystyle=\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},$

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

$A=\begin{bmatrix}A^{s}&V\\\ W&A^{f}\end{bmatrix}.$

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

$\mu_{n}=\begin{bmatrix}\mu^{s}_{n}\\\ \mu^{f}_{n}\end{bmatrix}.$

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

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

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

$\mu_{n,1}=(I+A\delta t)\mu_{n},\ \ \Sigma_{n,1}=(I+\delta
tA)\Sigma_{n}(I+\delta tA)^{T}+\delta tBB^{T},$ (7)

after which we extrapolate the slow mean as

$\begin{split}\mu^{s}_{n+1}&=\mu^{s}_{n}+\frac{\Delta t}{\delta
t}\left(\mu_{n,1}^{s}-\mu_{n}^{s}\right)\\\ &=\mu^{s}_{n}+\frac{\Delta
t}{\delta t}\left(\mu^{s}_{n}+\delta tA^{s}\mu^{s}_{n}+\delta
tV\mu^{f}_{n}-\mu^{s}_{n}\right)\\\ &=(I^{s}+\Delta tA^{s})\mu_{n}^{s}+\Delta
tV\mu^{f}_{n}.\end{split}$ (8)

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

$\begin{split}\mu^{f}_{n+1}&=\mu^{f}_{n,1}+C^{T}_{n,1}\left(\Sigma_{n,1}^{s}\right)^{-1}\left(\mu^{s}_{n+1}-\mu^{s}_{n,1}\right)\\\
&=\delta tW\mu^{s}_{n}+(I^{f}+\delta
tA^{f})\mu^{f}_{n}+C^{T}_{n,1}\left(\Sigma_{n,1}^{s}\right)^{-1}\left((\Delta
t-\delta t)A^{s}\mu^{s}_{n}+(\Delta t-\delta t)V\mu^{f}_{n}\right)\\\
&=\left[\delta tW+(\Delta t-\delta
t)C^{T}_{n,1}\left(\Sigma_{n,1}^{s}\right)^{-1}A^{s}\right]\mu^{s}_{n}+\left[I^{f}+\delta
tA^{f}+(\Delta t-\delta
t)C^{T}_{n,1}\left(\Sigma_{n,1}^{s}\right)^{-1}V\right]\mu^{f}_{n}.\end{split}$
(9)

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

$\begin{bmatrix}\mu^{s}_{n+1}\\\
\mu^{f}_{n+1}\end{bmatrix}=\begin{bmatrix}I^{s}+\Delta tA^{s}&\Delta tV\\\
\delta tW+(\Delta t-\delta
t)C^{T}_{n,1}\left(\Sigma_{n,1}^{s}\right)^{-1}A^{s}&I^{f}+\delta
tA^{f}+(\Delta t-\delta
t)C^{T}_{n,1}\left(\Sigma_{n,1}^{s}\right)^{-1}V\end{bmatrix}\begin{bmatrix}\mu^{s}_{n}\\\
\mu^{f}_{n}\end{bmatrix}.$ (10)

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,

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

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

$\displaystyle\mathcal{D}(P_{n_{\Delta t}(T)}\,||\,P_{T})$
$\displaystyle=\frac{1}{2}\Big{(}\ln\frac{\left|\Sigma_{T}\right|}{\left|\Sigma_{n_{\Delta
t}(T)}\right|}-d+\text{Tr}\big{(}\Sigma_{T}^{-1}\Sigma_{n_{\Delta
t}(T)}\big{)}+(\mu_{T}-\mu_{n_{\Delta
t}(T)})^{T}\Sigma_{T}^{-1}(\mu_{T}-\mu_{n_{\Delta t}(T)})\Big{)},$

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

$\displaystyle\lim\limits_{\Delta t\to\delta
t}\begin{bmatrix}\mu^{s}_{n(\Delta t)}\\\ \mu^{f}_{n(\Delta
t)}\end{bmatrix}=\begin{bmatrix}I^{s}+\delta tA^{s}&\delta tV\\\ \delta
tW&I^{f}+\delta tA^{f}\end{bmatrix}^{n(\delta t)}\begin{bmatrix}\mu^{s}_{0}\\\
\mu^{f}_{0}\end{bmatrix},$

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

$\displaystyle\underset{\delta t\to 0}{\lim}\
\frac{1}{2}\Big{(}\ln\frac{\left|\Sigma_{T}\right|}{\left|\Sigma_{n_{\delta
t}(T)}\right|}-d+\text{Tr}\big{(}\Sigma_{T}^{-1}\Sigma_{n_{\delta
t}(T)}\big{)}+(\mu_{T}-\mu_{n_{\delta
t}(T)})^{T}\Sigma_{T}^{-1}(\mu_{T}-\mu_{n_{\delta t}(T)})\Big{)},$

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

$V_{\infty}^{\delta t}=\delta t\sum_{j=0}^{\infty}(I+\delta
tA)^{j}BB^{\\!\mathsf{T}}(I+\delta tA^{\\!\mathsf{T}})^{j},$ (12)

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

$A=\begin{bmatrix}A^{s}&0\\\ 0&A^{f}\end{bmatrix}$

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

$\lim_{n\to\infty}\mu_{n}=0,\quad\lim_{n\to\infty}\Sigma_{n}=V_{\infty}^{\delta
t},$

whenever

$\rho(I^{s}+\Delta tA^{s})<1\quad\text{and}\quad\rho(I^{s}+\delta tA^{f})<1.$
(13)

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 $0$, 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$

$\mathcal{K}_{P}(\theta)=\ln\mathbb{E}_{P}\big{[}e^{\theta\mathbin{\mathbf{\cdot}}\Pi}\big{]},\quad\theta\in\Theta_{P},$

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

$\mathcal{K}_{X}(\theta)=\mu\mathbin{\mathbf{\cdot}}\theta+\frac{1}{2}\theta^{\\!\mathsf{T}}\Sigma\theta.$

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

$\mathbb{P}(\theta\mathbin{\mathbf{\cdot}}X_{0}\geq r)\leq
e^{-sr+\mathcal{K}_{\theta\mathbin{\mathbf{\cdot}}X_{0}}(s)},$ (14)

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

$-\ln\mathbb{P}(\theta\mathbin{\mathbf{\cdot}}X_{0}\geq r)\geq
sr-\mathcal{K}_{0}(s\theta),$

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

$\lim_{r\to+\infty}\frac{-\ln\mathbb{P}(\theta\mathbin{\mathbf{\cdot}}X_{0}\geq
r)}{r}\geq s.$

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

$\displaystyle\begin{aligned}
\mathcal{K}_{n}(\theta)&=\mathcal{K}_{n-1}\big{(}(I+\delta
tA^{\\!\mathsf{T}})(\theta+\lambda^{s}_{n}\oplus
0^{f})\big{)}-\mathcal{K}_{n-1}\big{(}({I}+\delta
tA^{\\!\mathsf{T}})(\lambda^{s}_{n}\oplus 0^{f})\big{)}\\\\[1.99997pt]
&\hphantom{=\ }+\frac{\delta
t}{2}\big{[}\theta^{\\!\mathsf{T}}\\!BB^{\\!\mathsf{T}}\theta+(\lambda^{s}_{n}\oplus
0^{f})^{\\!\mathsf{T}}\\!BB^{\\!\mathsf{T}}\theta+\theta^{\\!\mathsf{T}}\\!BB^{\\!\mathsf{T}}(\lambda^{s}_{n}\oplus
0^{f})\big{]},\end{aligned}$ (15)

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

$\mathcal{K}_{n}(\theta)=\mathcal{K}_{n-1,1}(\theta+\lambda^{s}_{n}\oplus
0^{f})-\mathcal{K}_{n-1,1}(\lambda^{s}_{n}\oplus 0^{f}),$ (16)

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

$X^{\delta t}_{n-1,1}=({I}+\delta tA)X_{n-1}+B\delta W_{n,1}.$

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

$\mathcal{K}_{n-1,1}=\mathcal{K}_{(I+\delta tA)X_{n-1}}+\mathcal{K}_{B\delta
W_{n,1}}.$

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

$\mathcal{K}_{n-1,1}(\theta)=\mathcal{K}_{n-1}\big{(}({I}+\delta
tA^{\\!\mathsf{T}})\theta\big{)}+\frac{\delta
t}{2}\theta^{\\!\mathsf{T}}BB^{\\!\mathsf{T}}\theta.$ (17)

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

$\mathcal{K}_{0}(r\theta)\sim
w_{0}(\theta)\,r^{2},\quad\frac{\operatorname{d}\\!}{\operatorname{d}\\!{r}}\mathcal{K}_{0}(r\theta)\sim
2\,w_{0}(\theta)\,r,$ (18)

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

$-\ln\mathbb{P}(\theta\mathbin{\mathbf{\cdot}}X_{0}\geq
r)\sim-\ln\mathbb{P}(\theta\mathbin{\mathbf{\cdot}}X_{0}\leq-r)\sim
w_{0}(\theta)^{-1}r^{2},\qquad\text{as}\ r\to+\infty.$

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

$\sup_{\theta\neq
0}\frac{\|\nabla_{\\!}\mathcal{K}(\theta)\|}{\|\theta\|}<+\infty.$ (19)

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

$\frac{\mathcal{K}(r\theta+\theta_{0})-\mathcal{K}(r\theta)}{r^{2}}+\frac{\mathcal{K}(r\theta)}{r^{2}}.$

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

$\frac{|\mathcal{K}(r\theta+\theta_{0})-\mathcal{K}(r\theta)|}{r^{2}}\leq\frac{\|\nabla_{\\!}\mathcal{K}(r\theta+r^{\prime}\theta_{0})\|\|\theta_{0}\|}{r^{2}},$

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

$\frac{\|\nabla_{\\!}\mathcal{K}(r\theta+r^{\prime}\theta_{0})\|}{r}=\frac{\|\nabla_{\\!}\mathcal{K}(r\theta+r^{\prime}\theta_{0})\|}{\|r\theta+r^{\prime}\theta_{0}\|}\cdot\big{\|}\theta+\frac{r^{\prime}}{r}\theta_{0}\big{\|}$

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$,

$\operatorname{D}\\!^{j}\mathcal{K}_{n}(\theta)[\theta^{\prime}]=\operatorname{D}\\!^{j}\mathcal{K}_{n-1}\big{(}(I+\delta
tA^{\\!\mathsf{T}})(\theta+\lambda^{s}_{n}\oplus 0^{f})\big{)}\big{[}(I+\delta
tA^{\\!\mathsf{T}})\theta^{\prime}\big{]},$

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

$\operatorname{D}\\!^{j}\mathcal{K}_{n}(\theta)[\theta^{\prime}]=\operatorname{D}\\!^{j}\mathcal{K}_{0}\Big{(}(I+\delta
tA^{\\!\mathsf{T}})^{n}\theta+\sum_{k=1}^{n}(I+\delta
tA^{\\!\mathsf{T}})^{k}(\lambda^{s}_{k}\oplus 0^{f})\Big{)}\big{[}(I+\delta
tA^{\\!\mathsf{T}})^{n}\theta^{\prime}\big{]}.$

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

$\|\operatorname{D}\\!^{j}\mathcal{K}_{n}(\theta)\|_{\mathrm{mult}}\leq\max_{\theta^{\prime}\in\mathbb{R}^{d}}\|\operatorname{D}\\!^{j}\mathcal{K}_{0}(\theta^{\prime})\|_{\mathrm{mult}}\|I+\delta
tA^{\\!\mathsf{T}}\|^{jn},$

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

$\displaystyle\begin{aligned}
\frac{\mathcal{K}_{n}(r\theta)}{r^{2}}&=\frac{\mathcal{K}_{n-1}\big{(}r(I+\delta
tA^{\\!\mathsf{T}})\theta+(I+\delta tA^{\\!\mathsf{T}})(\lambda^{s}_{n}\oplus
0^{f})\big{)}}{r^{2}}-\frac{\mathcal{K}_{n-1}\big{(}(I+\delta
tA^{\\!\mathsf{T}})(\lambda^{n}_{n}\oplus 0^{f})}{r^{2}}\\\\[1.99997pt]
&\hphantom{=\ }+\frac{\delta
t}{2}\frac{r^{2}\theta^{\\!\mathsf{T}}BB^{\\!\mathsf{T}}\theta+r(\lambda^{s}_{n}\oplus
0^{f})^{\\!\mathsf{T}}BB^{\\!\mathsf{T}}\theta+r\theta^{\\!\mathsf{T}}BB^{\\!\mathsf{T}}(\lambda^{s}_{n}\oplus
0)}{r^{2}}.\end{aligned}$ (20)

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

$\nabla_{\\!}\mathcal{K}_{n}(\theta)=(I+\delta
tA)\nabla_{\\!}{\mathcal{K}_{n-1}}\big{(}(I+\delta
tA^{\\!\mathsf{T}})(\theta+\lambda^{s}_{n}\oplus 0^{f})\big{)}+\delta
t(\theta+\lambda^{s}_{n}\oplus 0^{f})^{\\!\mathsf{T}}\\!BB^{\\!\mathsf{T}},$
(21)

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

$w_{n}(\theta)=w\big{(}(I+\delta
tA^{\\!\mathsf{T}})^{n}\theta\big{)}+\frac{\delta
t}{2}\theta^{\\!\mathsf{T}}\Big{[}\sum_{j=0}^{n-1}(I+\delta
tA)^{j}BB^{\\!\mathsf{T}}(I+\delta tA^{\\!\mathsf{T}})^{j}\Big{]}\theta.$ (22)

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

$\lim_{r\to+\infty}\frac{\nabla_{\\!}\mathcal{K}_{n}(r\theta)\cdot\theta}{r}=2w_{n}(\theta),$
(23)

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

$w_{\infty}(\theta)=\frac{1}{2}\theta^{\\!\mathsf{T}}\Big{[}\delta
t\sum_{j=0}^{\infty}(I+\delta tA)^{j}BB^{\\!\mathsf{T}}(I+\delta
tA^{\\!\mathsf{T}})^{j}\Big{]}\theta=\frac{1}{2}\theta^{\\!\mathsf{T}}V_{\infty}^{\delta
t}\theta.$ (24)

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

$\int_{0}^{1}r\theta^{\\!\mathsf{T}}\nabla^{2}_{\\!}\mathcal{K}_{n}(tr\theta)\theta\,\operatorname{d}\\!{t}=\nabla_{\\!}\mathcal{K}_{n}(r\theta)\cdot\theta-\nabla_{\\!}\mathcal{K}_{n}(0)\cdot\theta.$

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

$\lim_{r\to+\infty}\frac{1}{r}\int_{0}^{r}\theta^{\\!\mathsf{T}}\nabla^{2}_{\\!}\mathcal{K}_{n}(t\theta)\theta\,\operatorname{d}\\!{t}=\lim_{r\to+\infty}\Big{(}\frac{\nabla_{\\!}\mathcal{K}_{n}(r\theta)\cdot\theta}{r}-\frac{\nabla_{\\!}\mathcal{K}_{n}(0)\cdot\theta}{r}\Big{)}=2w_{n}(\theta),$

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)

$\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}\theta^{\\!\mathsf{T}}\nabla^{2}_{\\!}\mathcal{K}_{n}(t\theta)\theta\,\operatorname{d}\\!{t}=2w_{\infty}(\theta)=\theta^{\\!\mathsf{T}}V^{\delta
t}_{\infty}\theta.$

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

$\lim_{n\to+\infty}\theta^{\\!\mathsf{T}}\nabla^{2}_{\\!}\mathcal{K}_{n}(0)\theta=\lim_{n\to+\infty}\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}\theta^{\\!\mathsf{T}}\nabla^{2}_{\\!}\mathcal{K}_{n}(t\theta)\theta\,\operatorname{d}\\!{t}=\theta^{\\!\mathsf{T}}V^{\delta
t}_{\infty}\theta,$

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
$0$, by assumption, $\nabla^{2}_{\\!}\mathcal{K}_{n}(0)$ converges to
$V_{\infty}^{\delta t}$, and all higher order derivatives to $0$, we obtain
the following limiting sequence of cumulants as $n$ goes to $+\infty$:

$(0,V_{\infty}^{\delta t},0,\dotsc).$ (25)

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

$\begin{cases}dX=-2(X+Y)dt+\sin(2\pi t)dt+dW_{x}\\\
dY=\frac{1}{\varepsilon}(X-Y)dt+\frac{1}{\sqrt{\varepsilon}}dW_{y}.\end{cases}$
(26)

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

$\begin{pmatrix}\mu_{X}(t)\\\
\mu_{Y}(t)\end{pmatrix}=e^{tM}\begin{pmatrix}\mu_{X_{0}}-A\\\
\mu_{Y_{0}}-C\end{pmatrix}+\begin{pmatrix}A\\\ C\end{pmatrix}\cos(2\pi
t)+\begin{pmatrix}B\\\ D\end{pmatrix}\sin(2\pi t),$ (27)

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

$\begin{pmatrix}-2\pi&2&0&2\\\ 2&2\pi&2&0\\\ 0&-10&-2\pi&10\\\
-10&0&10&2\pi\end{pmatrix}\begin{pmatrix}A\\\ B\\\ C\\\
D\end{pmatrix}=\begin{pmatrix}1\\\ 0\\\ 0\\\ 0\end{pmatrix}.$

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: The error in the slow mean of micro-macro acceleration as a function
of the extrapolation step size $\Delta t$, computed against the analytic
solution (27) (blue) and against the numerical result obtained by the Euler-
Maruyama method with time step $\delta t$ (orange) for $\varepsilon=0.5$
(left) and $\varepsilon=0.05$ (right). We clearly see that the micro-macro
acceleration error decreases when $\Delta t$ decreases, as given by Theorem 1.
Moreover, for small $\varepsilon$ there is almost no difference between the
error computed against the analytic solution and the microscopic time
integrator as the latter is very accurate. This difference is higher for
larger $\varepsilon$.

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

$\int_{G}x\exp\left(\bar{\lambda}^{s}\cdot
x-A(\bar{\lambda}^{s},\mu)\right)dP^{s}(x)=\bar{\mu}^{s},$

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.

Figure 2: The $(\delta t,\Delta t)$ stability plane of micro-macro
acceleration on the linear driven system (26). A blue dot indicates stability
and a red dot instability. When at least one matching failure occurs, we mark
an instability at the corresponding time step sizes. The stability domain has
a V-shaped domain, and for every microscopic time step, the maximal
extrapolation time step is larger than the deterministic stability bound of
0.3.

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.

## References

  * [1] J. Fish, Bridging the scales in nano engineering and science, Journal of Nanoparticle Research 8 (5) (2006) 577–594.
  * [2] M. O. Steinhauser, Computational Multiscale Modeling of Fluids and Solids, Springer, 2017.
  * [3] Y. W. Kwon, D. H. Allen, R. Talreja, Multiscale modeling and simulation of composite materials and structures, Vol. 47, Springer, 2008.
  * [4] T. S. Deisboeck, G. S. Stamatakos, Multiscale cancer modeling, CRC Press, 2010.
  * [5] T. Li, A. Abdulle, W. E, Effectiveness of implicit methods for stiff stochastic differential equations, Communications in Computational Physics 3 (2) (2008) 295–307.
  * [6] W. E, B. Engquist, Multiscale modeling and computation, Notices of the AMS 50 (9) (2003) 1062–1070.
  * [7] A. Abdulle, W. E, B. Engquist, E. Vanden-Eijnden, The heterogeneous multiscale method, Acta Numerica 21 (2012) 1–87.
  * [8] I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidis, O. Runborg, C. Theodoropoulos, Equation-free, coarse-grained multiscale computation: Enabling mocroscopic simulators to perform system-level analysis, Communications in Mathematical Sciences 1 (4) (2003) 715–762.
  * [9] I. G. Kevrekidis, G. Samaey, Equation-free multiscale computation: Algorithms and applications, Annual review of physical chemistry 60 (2009) 321–344.
  * [10] A. Abdulle, S. Cirilli, Stabilized methods for stiff stochastic systems, Comptes Rendus Mathematique 345 (10) (2007) 593–598.
  * [11] A. Abdulle, S. Cirilli, S-ROCK: Chebyshev methods for stiff stochastic differential equations, SIAM Journal on Scientific Computing 30 (2) (2008) 997–1014.
  * [12] K. Debrabant, G. Samaey, P. Zieliński, A micro-macro acceleration method for the Monte Carlo simulation of stochastic differential equations, SIAM Journal on Numerical Analysis 55 (6) (2017) 2745–2786.
  * [13] T. Lelièvre, G. Samaey, P. Zieliński, Analysis of a micro-macro acceleration method with minimum relative entropy moment matching, arXiv preprint arXiv:1801.01740.
  * [14] K. Debrabant, G. Samaey, P. Zieliński, Study of micro-macro acceleration schemes for linear slow-fast stochastic differential equations with additive noise, arXiv:1805.10219.
  * [15] H. Vandecasteele, P. Zieliński, G. Samaey, Efficiency of a micro-macro acceleration method for scale-separated stochastic differential equations, In preparation.
  * [16] R. M. Dudley, Real Analysis and Probability, Vol. 47, Cambridge University Press, 2002.
  * [17] J. R. Hershey, P. A. Olsen, Approximating the kullback leibler divergence between gaussian mixture models, in: 2007 IEEE International Conference on Acoustics, Speech and Signal Processing - ICASSP ’07, Vol. 4, 2007, pp. IV–317–IV–320.  doi:10.1109/ICASSP.2007.366913.
  * [18] T. M. Cover, J. A. Thomas, Elements of information theory, John Wiley & Sons, 2012\.
  * [19] N. H. Bingham, C. M. Goldie, J. L. Teugels, Regular Variation, Encyclopedia of Mathematics and its Applications, Cambridge University Press, 1987.  doi:10.1017/CBO9780511721434.
  * [20] A. Granata, The theory of higher-order types of asymptotic variation for differentiable functions. part i: Higher-order regular, smooth and rapid variation, Advances in Pure Mathematics 6 (12) (2016) 776.
  * [21] K. B. Athreya, S. N. Lahiri, Measure theory and probability theory, Springer Science & Business Media, 2006.
  * [22] B. Jørgensen, R. Labouriau, Exponential families and theoretical inference, Ph.D. thesis, University of Aarhus (1995).
  * [23] R. W. Keener, Theoretical statistics: Topics for a core course, Springer, 2011.

##  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}$

$\displaystyle\begin{aligned}
\Theta_{MX+c}&=\\{s\in\mathbb{R}^{l}:M^{\\!\mathsf{T}}s\in\Theta_{X}\\}\\\\[1.99997pt]
\mathcal{K}_{MX+c}(s)&=\mathcal{K}_{X}(M^{\\!\mathsf{T}}s)+s\mathbin{\mathbf{\cdot}}c,\quad
s\in\mathbb{R}^{l}.\end{aligned}$

  4. (iv)

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

$\mathcal{K}_{P}(\theta)=\sum_{j=1}^{+\infty}\frac{1}{j!}\operatorname{D}\\!^{j}\mathcal{K}_{P}(0)[\theta],$

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

$\mathbb{E}[X]=\nabla_{\\!}\mathcal{K}_{P}(0),\quad\mathbb{V}(X)=\nabla^{2}_{\\!}\mathcal{K}_{P}(0).$

  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

$\mathcal{K}_{Q}(\theta)=\mathcal{K}_{P}(\theta+\overline{\lambda}^{s}\\!\oplus
0^{f})-\mathcal{K}_{P}(\overline{\lambda}^{s}\\!\oplus 0^{f}),$

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

$\displaystyle\mathcal{K}_{Q}(\theta)$
$\displaystyle=\ln\mathbb{E}_{P}\Big{[}e^{\theta\mathbin{\mathbf{\cdot}}\Pi}\,\mathchoice{\frac{\operatorname{d}\\!Q}{\operatorname{d}\\!P}}{\operatorname{d}\\!Q/\\!\operatorname{d}\\!P}{\operatorname{d}\\!Q/\\!\operatorname{d}\\!P}{\operatorname{d}\\!Q/\\!\operatorname{d}\\!P}\Big{]}$
$\displaystyle=\ln\mathbb{E}_{P}\Big{[}e^{\theta\mathbin{\mathbf{\cdot}}\Pi+\overline{\lambda}^{s}\\!\mathbin{\mathbf{\cdot}}\Pi^{s}}\,e^{-A(\overline{\lambda}^{s},P^{s})}\Big{]}$
$\displaystyle=\ln\mathbb{E}_{P}\Big{[}e^{(\theta+\overline{\lambda}^{s}\\!\oplus
0^{f})\mathbin{\mathbf{\cdot}}\Pi}\Big{]}-A(\overline{\lambda}^{s},P^{s}).$

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

$A(\overline{\lambda}^{s},P^{s})=\ln\mathbb{E}_{P^{s}}\Big{[}e^{\overline{\lambda}^{s}\mathbin{\mathbf{\cdot}}\Pi^{s}}\Big{]}=\ln\mathbb{E}_{P}\Big{[}e^{(\overline{\lambda}^{s}\\!\oplus
0^{f})\mathbin{\mathbf{\cdot}}\Pi}\Big{]}.\qed$

######  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

$\Sigma\overline{\lambda}=\overline{\mu}-\mu.$

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

$\displaystyle\mathcal{K}_{Q}(\theta)$
$\displaystyle=\mu\mathbin{\mathbf{\cdot}}\theta+\frac{1}{2}\big{(}\theta^{\\!\mathsf{T}}\Sigma\overline{\lambda}+\overline{\lambda}^{\\!\mathsf{T}}\Sigma\theta\big{)}+\frac{1}{2}\theta^{\\!\mathsf{T}}\Sigma\theta$
$\displaystyle=\mu\mathbin{\mathbf{\cdot}}\theta+(\overline{\mu}-\mu)\mathbin{\mathbf{\cdot}}\theta+\frac{1}{2}\theta^{\\!\mathsf{T}}\Sigma\theta$
$\displaystyle=\overline{\mu}\mathbin{\mathbf{\cdot}}\theta+\frac{1}{2}\theta^{\\!\mathsf{T}}\Sigma\theta.\qed$