# Bias and Refinement of Multiscale Mean Field Models

Sebastian Allmeier<EMAIL_ADDRESS>0000-0003-4629-6348 and
Nicolas Gast  0000-0001-6884-8698<EMAIL_ADDRESS>InriaGrenobleFrance
Univ. Grenoble AlpesGrenobleFrance

###### Abstract.

Mean field approximation is a powerful technique which has been used in many
settings to study large-scale stochastic systems. In the case of two-timescale
systems, the approximation is obtained by a combination of scaling arguments
and the use of the averaging principle. This paper analyzes the approximation
error of this ‘average’ mean field model for a two-timescale model
$(\boldsymbol{X},\boldsymbol{Y})$, where the slow component $\boldsymbol{X}$
describes a population of interacting particles which is fully coupled with a
rapidly changing environment $\boldsymbol{Y}$. The model is parametrized by a
scaling factor $N$, e.g. the population size, which as $N$ gets large
decreases the jump size of the slow component in contrast to the unchanged
dynamics of the fast component. We show that under relatively mild conditions,
the ‘average’ mean field approximation has a bias of order $O(1/N)$ compared
to $\mathbb{E}[\boldsymbol{X}]$. This holds true under any continuous
performance metric in the transient regime, as well as for the steady-state if
the model is exponentially stable. To go one step further, we derive a bias
correction term for the steady-state, from which we define a new approximation
called the refined ‘average’ mean field approximation whose bias is of order
$O(1/N^{2})$. This refined ‘average’ mean field approximation allows computing
an accurate approximation even for small scaling factors, i.e., $N\approx
10-50$. We illustrate the developed framework and accuracy results through an
application to a random access CSMA model.

Mean Field; Refined Mean Field; Two Timescale System

††copyright: acmlicensed††journal: POMACS††journalyear: 2023††journalvolume:
7††journalnumber: 1††publicationmonth: 3††price: 15.00††doi:
10.1145/3579336††ccs: Mathematics of computing Markov processes††ccs:
Mathematics of computing Ordinary differential equations

##  1\. Introduction

The mean field approximation finds widespread application when interested in
analyzing the macroscopic behavior of large-scale stochastic systems composed
of interacting particles. Its assets lie in a reduction of the model
complexity, simplified analysis of the system due to absence of stochastic
components, and reduction of computation time compared to a stochastic
simulation. The mean field approximation can even yield closed form solutions
for the steady-state, e.g., for the well known JSQ(d) model (Mitzenmacher,
2001). The mean field approximation is generally given by a set of ordinary
differential equations which arise from the assumption that, for large system
sizes, the evolution of the particles are stochastically independent of
another. This idea works well if the number of particles is large and if the
particles can be clustered into a few groups with statistically identical
behavior (Kurtz, 1970, 1978; Benaim and Boudec, 2008; Le Boudec et al., 2007;
Gast and Gaujal, 2012). The framework established by Kurtz (Kurtz, 1970) to
derive (weak) convergence results for the stochastic system justifying the use
of the mean field approximations finds sustained attention in the literature.

More recently, the authors of (Gast, 2017; Ying, 2016) showed that for finite
system sizes the bias of the mean field approximation is of order $1/N$ when
compared to the mean behavior of the system. Here, $N$ is the scaling
parameter, which usually refers to the number of homogeneous particles in the
system. Additional works such as (Gast and Van Houdt, 2017; Gast et al., 2019)
introduced corrections, called refinement terms, which effectively increase
the rate of accuracy of the approximation and therefore the rate of
convergence. The most notable term is the first-order bias refinement, since
it offers a convincing trade-off between a significant accuracy gain and
additional computation cost.

While these classical mean field results hold for a broad class of models,
most of the results can not be transferred to systems with more intricate
dynamics such as the two-timescale case, studied for instance in (Benaim and
Boudec, 2008; Bordenave et al., 2010). A two-timescale process consists of two
coupled components, one evolving slowly compared to the other. The slowly
evolving component is often represented by a system of interacting particles,
where the state of each particle evolves as a function of the empirical
distribution of all particles but also as a function of the state of the fast
component, e.g., a fast changing environment. These types of processes and
their ‘averaged’ mean field adaptation have been of interest since the 1960s
and became increasingly relevant in the study of modern and complex systems.
We refer to (Pavliotis and Stuart, 2008) for an extensive literature
discussion. An important area of application comes from the field of computer
networks. Examples include loss networks (Hunt and Kurtz, 1994), large-scale
random access networks with interference graphs (Cecchi et al., 2021; Castiel
et al., 2021) or storage networks (Feuillet and Robert, 2014). Another recent
field of literature from which we draw inspiration are chemical reaction
networks. The works of Kang et al. (Kang and Kurtz, 2013; Kang et al., 2014)
and Ball et al. (Ball et al., 2006) establish central limit theorems for large
multiscale models motivated by biological and chemical processes. Another
application in the field of biology is given by (Robert and Vignoud, 2021) who
use the mean field idea to study neural plasticity models.

##### Contributions

The aforementioned papers underline that the mean field idea can be adapted to
two-timescale models and prove that the ‘average’ mean field approximation is
asymptotically exact as $N$ goes to infinity. They do not, however, provide
theoretical results which guarantee the accuracy or performance bounds of the
approximation for finite systems sizes. This paper aims at filling this gap.
We derive accuracy bounds for the ‘average’ mean field approximation in the
transient regime and steady-state which show that its bias is of order
$O(1/N)$, $N$ being the scaling parameter of the stochastic system. We further
derive bias correction terms for the steady-state, from which we define a new
approximation called the refined ‘average’ mean field approximation, whose
bias is of order $O(1/N^{2})$. To prove these accuracy bounds, we develop a
framework for two-timescale stochastic models whose slow component is
comparable to the concept of density dependent population processes as
introduced by Kurtz (Kurtz, 1978). Based on this representation, we utilize a
combination of generator comparison techniques, Poisson equations as well as
derivative bounds on the solution of the Poisson equation. This allows us to
bound the bias with respect to the scaling parameter $N$. To take it a step
further, we then prove the existence of correction terms which approximate the
bias of the ‘average’ mean field and allow defining the refined ‘average’ mean
field. To support the practical application of the refinements, we provide an
algorithmic way to compute the correction terms. This includes methods to
numerically solve the Poisson equation and obtain its derivatives. We
illustrate the computation of the refinement terms and confirm the accuracy of
the obtained bounds by considering a random access CSMA model. Using the
example we show that even for relatively small $N\approx 10$ the refined
‘average’ mean field approximation almost exactly indicates the steady-state
of the model.

##### Methodological Advances & Technical Challenges

To obtain our results, we build on the recent line of work on Stein’s method
(Stein, 1986). This method allows calculating the distance between two random
variables by looking at the distance between the generators of two related
systems. Recently, the method reemerged in publications within the stochastic
network community, in particular the works of Braverman et al. (Braverman et
al., 2016; Braverman and Dai, 2017; Braverman, 2022). In this paper, we use
the Poisson equation idea in two ways. First and foremost, we use a Poisson
equation, called the ‘fast’ Poisson equation, to bound the distance between a
function – in this case the drift – and its average version given by the
steady-state distribution of the fast process as introduced in Section 4.2.
This step is integral to our paper and constitutes the building block to
obtain the error bounds and closed form expressions as it allows to analyze
the distance between the coupled stochastic process and its decoupled
counterpart. The analysis however bears many technical intricacies one needs
to overcome. This includes solving the Poisson equation, stating its
derivative bounds and deducing computable expressions used to calculate the
bias term. Second, for our steady-state results, we make use of another
Poisson equation to compare the stochastic system to the equilibrium point of
the approximation. The latter is related to the methods used in (Ying, 2016;
Gast and Van Houdt, 2017) but significantly extends the ideas as the
derivation of the bias exhibits novel refinement terms which originate from
the coupling and correct the error of the ‘averaging’ method used for the mean
field. The closely linked technical challenges include obtaining numerically
feasible formulas for the ‘new’ refinement terms. Here, we utilize the derived
solution of the ‘fast’ Poisson equation to further specify the fluctuations of
the stochastic system around the equilibrium point. Carefully analyzing the
combination of the two Poisson equations enables us to obtain the new
refinement terms. The closed-form bias terms that we obtained from the steady-
state analysis are significantly more complicated than the ones of (Gast and
Van Houdt, 2017) as they further correct the error made by the averaging
method.

##### Applicability & Numerical Difficulties

In our paper, we make use of a CSMA model to demonstrate the applicability of
our results. There are several other examples captured by our framework, such
as the Michaelis-Menten enzyme model of (Kang et al., 2014) or the storage
network investigated in (Feuillet and Robert, 2014). For the latter, the
authors observe that when looking at the right timescale the loss of the
network can be characterized by a local equilibrium which is obtained using
the averaging method. The biochemical Michaelis-Menten enzyme model describes
the dynamics between three time-varying species, the enzyme, a substrate and a
product. Following the description of the model as in (Kang et al., 2014) and
by using the right scaling arguments, the model exhibits two timescales, the
fast reacting and state changing enzymes and the slower changing
concentrations of the substrate and product. Our framework can be used on both
examples to guarantee accuracy results and can be used to compute refined
approximations.

To compute the ‘average’ mean field and refinement terms, one has to overcome
numerical difficulties which arise from the averaging method. First, in order
to compute the ‘average’ drift, one needs to compute the steady-state
probabilities of the fast system, which might not be available in the closed
form as for the CSMA model. For this case we provide computational notes in
the appendix which aim to facilitate the computation. Another problem for the
refinement term is the need for the first and second derivatives of the
‘average’ drift. For the CSMA model we used symbolic representation of the
transition rates from which we define the drift and its average version. This
method allows to numerically compute the derivatives using ‘sympy’, a Python
library for symbolic computation. This method is relatively easy to implement
but has a very large computation time. This computation time could be reduced
by implementing a faster computation of the derivatives. In fact, to compute
the derivative of the matrix $K^{+}$ which is closely linked to the solution
of the ‘fast’ Poisson, we provide supplementary computational notes which
describe how to obtain an efficient implementation.

##### Roadmap

The paper is organized as follows. In Section 2, we formally introduce the
two-timescale model, its corresponding ‘average’ mean field approximation and
make regularity assumptions on the system. In Section 3 we state the main
results of this paper. Section 3.1 states the results for the transient
regime, Section 3.2 for the steady-state and Section 3.3 justifies the
existence of bias correction terms. Section 4 holds the proofs of the
aforementioned results. In Section 5 we apply our theoretical results to the
unsaturated random-access network model of (Cecchi, 2018; Cecchi et al.,
2021). Some technical lemmas and definitions are postponed to the appendix.

##### Reproducibility

The code to reproduce the paper along with all figures and the implementation
of the unsaturated random-access network model is available at
https://gitlab.inria.fr/sallmeie/bias-and-refinement-of-multiscale-mean-field-
models.

##  2\. Stochastic System and Mean Field Approximation

We consider a two-timescale, coupled, continuous time Markov chain
$(\boldsymbol{X}_{s}^{(N)},\boldsymbol{Y}_{s}^{(N)})_{s\geq 0}$ parametrized
by a scaling factor $N$, for which we study the behavior as $N$ tends to
infinity. As we will see in the examples, $N$ typically represents the number
of objects that interact together. This section introduces the precise model
and fixes notations.

###  2.1. Model

For a fixed scaling parameter $N$, the stochastic process
$(\boldsymbol{X}^{(N)}_{s},\boldsymbol{Y}^{(N)}_{s})_{s\geq 0}$ is a
continuous time Markov chain that evolves in a state-space
$\mathcal{X}^{(N)}\times\mathcal{Y}$. The set $\mathcal{Y}$ is finite and does
not depend on $N$. We further require that for all $N$, the sets
$\mathcal{X}^{(N)}$ are subsets of a convex and compact set
$\mathcal{X}\subset\mathbb{R}^{d_{x}}$. In what follows, unless it is
ambiguous in the context, we drop the dependence on $N$ to lighten the
notations.

This model has two-timescales in the sense that the size order of jumps of
$\boldsymbol{Y}_{s}$ is $N$-times larger than the ones of
$\boldsymbol{X}_{s}$. More precisely, we assume that there exists a finite
number of transitions $(\ell,\boldsymbol{y}^{\prime})\in\mathcal{T}$ with
their corresponding transition rate functions
$\alpha_{\ell,\boldsymbol{y}^{\prime}}\geq 0$, both being independent of $N$,
such that for all possible states
$(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})\in\mathcal{X}\times\mathcal{Y}$:

(1) $\displaystyle(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})\text{ jumps to
}(\boldsymbol{X}_{s}+\ell/N,\boldsymbol{y}^{\prime})$ $\displaystyle\text{ at
rate
}N\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{s},\boldsymbol{Y}_{s}).$

The above defines continuous time Markov chains with discrete state-space
whose realizations are Càdlàg, i.e., right continuous with a left limit for
every time $t$. Note that in Equation (1), we assume that the transition rates
$\alpha$ are defined for all $x\in\mathcal{X}$ (and not just for
$x\in\mathcal{X}^{(N)}$).

The notion of slow-fast system comes from the fact that the jumps of the slow
component $\mathcal{X}$ are $O(N)$ times smaller than the jumps of the fast
component $\mathcal{Y}$ while transition rates are of the same scale. As we
will see later, the different timescales imply that for large $N$, the slow
component $\boldsymbol{X}_{s}$ will ’see’ the fast component
$\boldsymbol{Y}_{s}$ as if it is stationary with distribution
$(\pi_{\boldsymbol{y}}(\boldsymbol{X}_{s}))_{\boldsymbol{y}\in\mathcal{Y}}$
which we formally define in Section 2.2.

###  2.2. Drift, Average Drift and Mean Field Approximation

The jumps of the stochastic system (1) can affect the fast and/or the slow
component. In what follows, we construct an approximation that consists (i) in
considering that the slow component is not stochastic but evolves
deterministically according to its drift (which is its average change), and
(ii) in using a time-averaging method that shows the stochastic process
$\boldsymbol{Y}$ as being in some stationary state given $\boldsymbol{x}$.
This leads us to two definitions:

  1. (i)

We call the drift of the slow system (or more concisely the drift) the average
change of $\boldsymbol{X}_{s}$. It is the sum over all possible transitions of
the rate of transition multiplied by the changes that such transitions induce
on $\boldsymbol{X}_{t}$. By the form of the transitions in (1), if the process
starts in
$(\boldsymbol{X}_{t},\boldsymbol{Y}_{t})=(\boldsymbol{x},\boldsymbol{y})$, the
drift is given by:

(2) $\displaystyle
F(\boldsymbol{x},\boldsymbol{y}):=\sum_{\ell,\boldsymbol{y}^{\prime}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})\ell\in\mathbb{R}^{d_{x}}.$

This drift function depends on the state of the fast system $\boldsymbol{y}$.

  2. (ii)

For the fast component, we define a transition kernel
$K_{\boldsymbol{y},\boldsymbol{y}^{\prime}}(\boldsymbol{x})$ that is the rate
at which the process $\boldsymbol{Y}$ jumps from $\boldsymbol{y}$ to
$\boldsymbol{y}^{\prime}\neq\boldsymbol{y}$ (divided by $N$), with the usual
convention that
$K_{\boldsymbol{y},\boldsymbol{y}}(\boldsymbol{x})=-\sum_{\boldsymbol{y}^{\prime}\neq\boldsymbol{y}}K_{\boldsymbol{y},\boldsymbol{y}^{\prime}}(\boldsymbol{x})$:

(3) $\displaystyle
K_{\boldsymbol{y},\boldsymbol{y}^{\prime}}(\boldsymbol{x})=\sum_{\ell}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y}).$

As $\mathcal{Y}$ is finite, for a fixed $\boldsymbol{x}$, $K(\boldsymbol{x})$
is a matrix that corresponds to the kernel of a continuous time Markov chain.
Our assumptions will imply that for all $\boldsymbol{x}$, the process
associated with $K(\boldsymbol{x})$ has a unique stationary distribution,
which we denote by the vector
$\pi(\boldsymbol{x})=(\pi_{\boldsymbol{y}}(\boldsymbol{x}))_{\boldsymbol{y}\in\mathcal{Y}}$.

Based on the drift (2) of the stochastic system, we define its ‘average’
version $\bar{F}$ by averaging over the stationary distribution of the fast
component. That is:

$\displaystyle\bar{F}(\boldsymbol{x}):=\sum_{y}\pi_{y}(\boldsymbol{x})F(\boldsymbol{x},\boldsymbol{y}).$

For an initial state $\boldsymbol{x}$ and $t\geq 0$, we call the mean field
approximation the solution $\phi_{t}(\boldsymbol{x})$ of the initial value
problem

(4)
$\displaystyle\frac{d}{dt}\phi_{t}(\boldsymbol{x})=\bar{F}(\phi_{t}(\boldsymbol{x})),\qquad\phi_{0}(\boldsymbol{x})=\boldsymbol{x}.$

Such an approximation is also called a fluid approximation.

###  2.3. Main assumptions

As we show later, under mild regularity conditions on the transition rate
functions $\alpha$, the mean field approximation captures the dynamics of
$\boldsymbol{X}_{t}$ well and with a decreasing bias of order $1/N$:

$\displaystyle\mathbb{E}[h(\boldsymbol{X}_{t})\mid\boldsymbol{X}_{0},\boldsymbol{Y}_{0}=\boldsymbol{x},\boldsymbol{y}]-h(\phi_{t}(\boldsymbol{x}))=\frac{C_{h}(t)}{N}+o(1/N),$

for a sufficiently regular $h$. This holds for any finite $t$ under assumption
($A_{\text{1}}$)-($A_{\text{2}}$) below. It also holds for the steady-state
regime $t=+\infty$ under the additional assumption ($A_{\text{3}}$). For the
steady-state we also show that $C_{h}$ can be computed numerically and use it
to propose a refined approximation. To obtain these results for finite time,
we will assume that:

  1. ($A_{\text{1}}$)

The set of transitions $\mathcal{T}$ is finite and for all
$\ell,\boldsymbol{y}^{\prime}\in\mathcal{T}$,
$\alpha_{\ell,\boldsymbol{y}^{\prime}}\in\mathcal{D}^{2}(\mathcal{X}\times\mathcal{Y})$,
where $\mathcal{D}^{k}(\mathcal{X}\times\mathcal{Y})$ the space of functions
from $\mathcal{X}\times\mathcal{Y}$ to $\mathbb{R}$ for which the Hölder norm
$\|h\|_{k,1}$ is finite111i.e., functions that are $k$-times differentiable
with Lipschitz continuous derivatives (see a more precise definition of this
norm in Definition 2 of Appendix 1).

  2. ($A_{\text{2}}$)

For all $\boldsymbol{x}\in\mathcal{X}$, the matrix $K(\boldsymbol{x})$ defined
in Equation (3) has a unique irreducible class.

  3. ($A_{\text{3}}$)

The ODE (4) has a unique, exponentially stable equilibrium, that we denote by
$\phi_{\text{\tiny$\infty$}}$, i.e., there exist $a,b>0$ s.t.
$\left\|\phi_{t}(\boldsymbol{x})-\phi_{\text{\tiny$\infty$}}\right\|\leq
a\exp(-bt)$ for all $\boldsymbol{x}\in\mathcal{X}$.

Assumptions ($A_{\text{1}}$) and ($A_{\text{3}}$) are classical to ensure mean
field convergence results. As stated in (Gast and Van Houdt, 2017; Ying,
2016), requiring that the transition rates are twice differentiable which is
necessary to guarantee the existence of derivatives for the drift and the
differential equation needed for the proofs of the theorems. Assumption
($A_{\text{3}}$) ensures the existence of a unique equilibrium point to which
all trajectories of the differential equation converge. This classical
assumption is essentially needed to show that the stationary distribution of
the stochastic process converges to a deterministic limit, see (Benaim and
Boudec, 2008). It guarantees the stability of both the ODE and the ’slow’
Poisson equation used in the proof.

By Assumption ($A_{\text{2}}$), we mean that for all $\boldsymbol{x}$ the
Markov chain should have a unique subset of states that is irreducible (there
can be additional states, but they should all be transient). This assumption
is equivalent to assuming the uniqueness of the stationary distribution of the
Markov chain induced by the generator matrix $K(\boldsymbol{x})$ which is
essential to define the ‘averaged’ drift and mean field approximation. For a
given $\boldsymbol{x}$, the stationary distribution $\pi(\boldsymbol{x})$ will
be non-zero for all states that are in the irreducible component
($\pi_{\boldsymbol{y}}(\boldsymbol{x})>0$ for such $\boldsymbol{y}$’s) and
will be zero for the others ($\pi_{\boldsymbol{y}}(\boldsymbol{x})=0$ for all
states that are transient for $K(\boldsymbol{x})$). This assumption is
slightly more general than assuming that $K(\boldsymbol{x})$ is irreducible
because it allows for transient states.

We will show later that the assumptions ($A_{\text{1}}$) and ($A_{\text{2}}$)
imply that $\bar{F}$ is Lipschitz continuous because they imply that
$\pi(\boldsymbol{x})$ is Lipschitz-continuous (see Lemma 1). This implies that
the mean field approximation (4) is well defined.

##  3\. Main Results

This section includes our main results which are threefold. In 3.1 we obtain
accuracy results for the mean field approximation in the transient regime. In
3.2 we obtain comparable results when the stochastic system is in its steady-
state. Lastly, in 3.3 we introduce a correction term for the steady-state,
give accuracy bounds and display closed form expressions of the corrections.

###  3.1. Transient Regime

Our model is a two timescale model, which makes it amenable to be analyzed by
time-averaging methods such as the one used in (Benaim and Boudec, 2008;
Robert and Vignoud, 2021; Cecchi et al., 2021; Castiel et al., 2021). Such
methods guarantee that the stochastic process $\boldsymbol{X}_{t}$ converges
to the solution of the ODE given by (4). Yet, most of the papers that use
averaging methods do not quantify the rate at which this convergence occurs.
Our first result, Theorem 1, states that the difference between the ‘average’
mean field approximation $\phi$ derived from the average drift $\bar{F}$
approximates the average behavior of the slow component of the two-timescale
system with a bias asymptotically equivalent to $C_{h}(t)/N$. This result is
similar to the one obtained in (Gast et al., 2019) for classical density
dependent process.

###### Theorem 1.

Consider the two-timescale stochastic system
$(\boldsymbol{X}^{(N)},\boldsymbol{Y}^{(N)})_{t\geq 0}$ as introduced in
Section 2 starting at
$\boldsymbol{X}_{0}^{(N)},\boldsymbol{Y}_{0}^{(N)}=\boldsymbol{x},\boldsymbol{y}$,
and recall that $\phi(\boldsymbol{x})$ is the solution of the ODE (4) with
initial condition $\boldsymbol{x}$. Further, assume ($A_{\text{1}}$) and
($A_{\text{2}}$). For any $h\in\mathcal{D}^{2}(\mathcal{X})$ and $t>0$, there
exists a constant $C_{h}(t)$ such that:

$\displaystyle\lim_{N\rightarrow\infty}N\left(\mathbb{E}[h(\boldsymbol{X}^{(N)}_{t})\mid\boldsymbol{X}_{0}^{(N)},\boldsymbol{Y}_{0}^{(N)}=\boldsymbol{x},\boldsymbol{y}]-h\bigl{(}\phi_{t}(\boldsymbol{x})\bigr{)}\right)=C_{h}(t).$

###### Main element of the proofs.

The full proof of this theorem is given in Section 4.3. It is decomposed into
two parts that correspond to the two approximations (scaling and averaging):

  * •

The first is to approximate the infinitesimal generator
$Lh(\boldsymbol{X}_{t},\boldsymbol{Y}_{t})$ (defined in Section 4.1) of the
stochastic system by a process whose drift is $D_{x}h(\boldsymbol{X}_{t})\cdot
F(\boldsymbol{X}_{t},\boldsymbol{Y}_{t})$. This part uses that the rate
functions $\alpha$ and $h$ are differentiable with respect to the slow
variable $\boldsymbol{x}$.

  * •

The second is that the actual drift of the mean field approximation is
$\bar{F}(\boldsymbol{x})=\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\boldsymbol{x})f(\boldsymbol{x},\boldsymbol{y})$
and not the $F(\boldsymbol{x},\boldsymbol{y})$ obtained from the stochastic
system. This leads us to bound a term of the form
$\bar{F}(\boldsymbol{x})-F(\boldsymbol{x},\boldsymbol{y})$.

The first error term is bounded by using generator techniques similar to the
ones used for classical mean field models as in (Gast, 2017). It is treated in
Section 4.3.1. To bound the second term, we use averaging ideas similar to the
ones of (Benaim and Boudec, 2008) that are related to how fast the fast
timescale converges to its stationary distribution. This is dealt with in
Section 4.3.2. ∎

###  3.2. Steady-State Results

Theorem 1 guarantees that the mean field is a good approximation for any
finite time interval. In order to obtain a similar result for the stationary
case, an almost necessary condition is that the ODE (4) has a unique fixed
point to which all trajectories converge (Benaim and Boudec, 2008). To obtain
the equivalent of Theorem 1, we assume that this unique fixed point is
exponentially stable, as is classically done to obtain steady-state guarantees
(Gast and Van Houdt, 2017; Ying, 2016). This assumption is summarized in
($A_{\text{3}}$) and leads to the following result.

###### Theorem 2.

Assume ($A_{\text{1}}$) \- ($A_{\text{3}}$), and assume that for all $N$, the
stochastic system has a stationary distribution. Denote by
$(\boldsymbol{X}_{\infty}^{(N)},\boldsymbol{Y}_{\infty}^{(N)})$ a pair of
variables having this stationary distribution. Then, for
$h\in\mathcal{D}^{2}(\mathcal{X}\times\mathcal{Y})$ there exists a constant
$C_{h}$ such that:

$\displaystyle\lim_{N\to\infty}N\left(\mathbb{E}[h(\boldsymbol{X}_{\infty}^{(N)},\boldsymbol{Y}_{\infty}^{(N)})]-\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\phi_{\text{\tiny$\infty$}})h(\phi_{\text{\tiny$\infty$}},\boldsymbol{y})\right)=C_{h}.$

###### Main elements of the proof.

A detailed proof is provided in Section 4.4. To prove the result, the first
essential step is to generate two Poisson equations:

  1. (1)

The first is used to study
$\mathbb{E}[h(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})-\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\boldsymbol{X}_{\infty})h(\boldsymbol{X}_{\infty},\boldsymbol{y})]$,
which is the difference between the true expectation and a hypothetical case
where $\boldsymbol{Y}_{\infty}$ would be independent from
$\boldsymbol{X}_{\infty}$ and distributed according to the stationary
distribution $\pi(\boldsymbol{X}_{\infty})$.

  2. (2)

The second is used to study
$\mathbb{E}[\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\boldsymbol{X}_{\infty})h(\boldsymbol{X}_{\infty},\boldsymbol{y})-\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\phi_{\text{\tiny$\infty$}})h(\phi_{\text{\tiny$\infty$}},\boldsymbol{y})]$,
which is the error of the ’slow’ process compared to the ‘average’ mean field.

The rest of the proof treats these Poisson equations to establish the constant
$C_{h}$. This is done by generator comparison techniques and application of
derivative bounds for the solutions of the Poisson equation. ∎

We note that contrary to Theorem 1, Theorem 2 allows for functions $h$ that
can depend on both the slow and the fast components: $\boldsymbol{x}$ and
$\boldsymbol{y}$. In fact, Theorem 1 would not be true if we allow functions
$h$ to depend on $\boldsymbol{Y}_{t}$ as in Theorem 2, because of the fast
transitions of $\boldsymbol{Y}_{t}$. For the steady-state, we use that
$(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})$ starts in steady-state.

In Theorem 3 below, we will show that, for the steady-state and for functions
that only depend on $\boldsymbol{x}$, we can go further and propose an almost
closed-form expression for the term $C_{h}$. This new bias term provides a
refined accuracy of order $O(1/N^{2})$. To prove this result, we will use that
Theorem 2 allows functions $h$ which depend on $\boldsymbol{X}$ and
$\boldsymbol{Y}$ to show that this new approximation is $O(1/N^{2})$-accurate.
We show below that in fact the constant $C_{h}$ can be computed by a numerical
algorithm, and can therefore be used to define a _refined approximation_ ,
similarly to what is done for classical mean field model in (Gast and Van
Houdt, 2017).

###  3.3. Steady-State Refinement

To obtain a refined approximation, we utilize ideas introduced in (Gast and
Van Houdt, 2017) and propose an almost-closed form expression for the term
$C_{h}$ of Theorem 2. As we will see later in the proofs, the bias correction
term is composed of two distinct components:

  1. (1)

The first one (terms $V$ and $W$) is the analogue of the $V$ and $W$ terms of
(Gast and Van Houdt, 2017) and corresponds to the difference between the
stochastic jumps of the slow system versus having a ODE corresponding to the
(non-averaged) drift $f(\boldsymbol{x},\boldsymbol{y})$.

  2. (2)

The second component (terms $T$, $S$ and $U$) corresponds to approximating
$\boldsymbol{Y}$ by its stationary distribution
$\pi_{\boldsymbol{y}}(\boldsymbol{x})$, and its consequence on the behavior of
the slow system $\boldsymbol{X}$.

The proofs of the expression for $V$ and $W$ are essentially those derived in
(Gast and Van Houdt, 2017). The second correction component (terms $T$, $S$
and $U$) is related to the difference between drift and its average version.
It involves studying the intricate and coupled dynamics of $\boldsymbol{X}$
and $\boldsymbol{Y}$ which, to the best of our knowledge, has not been studied
and yields novel results.

###### Theorem 3.

Assume ($A_{\text{1}}$) \- ($A_{\text{3}}$). Then, there exist vectors $V$,
$S$, $T$ and matrices $W$ and $U$ that are solutions of linear systems of
equations such that for $h\in\mathcal{D}^{2}(\mathcal{X})$, we define $C_{h}$
of Theorem 2 as:

$\displaystyle C_{h}=\sum_{i}\frac{\partial h}{\partial
x_{i}}(\phi_{\text{\tiny$\infty$}})\left(V_{i}+T_{i}+S_{i}\right)+\frac{1}{2}\sum_{i,j}\frac{\partial^{2}h}{\partial
x_{i}x_{j}}(\phi_{\text{\tiny$\infty$}})\big{(}W_{i,j}+U_{i,j}\big{)}.$

###### Proof.

This theorem is a consequence of Proposition 2 that shows the existence of the
above constants, combined with Proposition 3 which provides the linear
equations satisfied by $V$ and $W$, and Proposition 4 which provides the
equations satisfied by $T$, $U$ and $V$. ∎

We would like to emphasize that the result of Theorem 3 is only valid for
functions $h$ that do _not_ depend on $\boldsymbol{y}$. This shows that there
exists a computable constant $C_{h}$ such that, in steady-state:

$\displaystyle\mathbb{E}[h(\boldsymbol{X}_{\infty})]=\phi_{\text{\tiny$\infty$}}+\frac{C_{h}}{N}+o(1/N).$

Similar to (Gast and Van Houdt, 2017), we call
$\phi_{\text{\tiny$\infty$}}+\frac{C_{h}}{N}$ the _refined approximation_ of
$\mathbb{E}[h(\boldsymbol{X}_{\infty})]$. As we see in our numerical
experiment in Section 5, this constant is computable and can provide a more
accurate approximation than the classical mean field.

In fact, combined with Theorem 2, we can show that the $o(1/N)$ term is a
$O(1/N^{2})$ term. This shows that the refined approximation is
$O(1/N^{2})$-accurate. Note, for functions $h$ that do depend on
$\boldsymbol{y}$, the existence of such a function is guaranteed by Theorem 2
but a closed-form expression is currently out of our reach.

##  4\. Proofs

###  4.1. Stochastic Semi-Groups and Generators

Given the stochastic process of Section 2.1, we define the stochastic semi-
group operator which maps a pair of initial values
$(\boldsymbol{x},\boldsymbol{y})$ and a function $h$ to the expected value of
the system at time $t$. This semi-group $\psi_{t}$ associates to a function
$h:\mathcal{X}^{(N)}\times\mathcal{Y}\to\mathbb{R}$ the function
$\psi_{t}h:\mathcal{X}^{(N)}\times\mathcal{Y}\to\mathbb{R}$ defined as:

(5)
$\displaystyle\psi_{t}h(\boldsymbol{x},\boldsymbol{y})=\mathbb{E}[h(\boldsymbol{X}_{t},\boldsymbol{Y}_{t})\mid
X_{0},Y_{0}=\boldsymbol{x},y].$

Using the right continuity of the slow-fast system, it is easy to verify that
$\psi_{t}$ is indeed a $C_{0}$-semi-group (see a more precise definition in
Definition 1 in Appendix 1).

The infinitesimal generator of the stochastic process is the operator $L$ that
maps a function $h\in\mathcal{D}^{k}(\mathcal{X}\times\mathcal{Y})$ to
$Lh:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}$ defined by:

(6) $\displaystyle Lh(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\sum_{\ell,\boldsymbol{y}^{\prime}\in\mathcal{T}}N\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})\left(h(x+\ell/N,\boldsymbol{y}^{\prime})-h(\boldsymbol{x},\boldsymbol{y})\right).$

Note that $Lh$ is obtained by considering average infinitesimal change of the
stochastic system starting in $(\boldsymbol{x},\boldsymbol{y})$, i.e.,

$\displaystyle Lh(\boldsymbol{x},\boldsymbol{y})=\lim_{t\downarrow
0}\bigl{(}\psi_{t}h(\boldsymbol{x},\boldsymbol{y})-\psi_{0}h(\boldsymbol{x},\boldsymbol{y})\bigr{)}/t=\lim_{t\downarrow
0}\bigl{(}\mathbb{E}[h(\boldsymbol{X}_{t},\boldsymbol{Y}_{t})\mid
X_{0},Y_{0}=\boldsymbol{x},\boldsymbol{y}]-h(\boldsymbol{x},\boldsymbol{y})\bigr{)}/t.$

Similarly to the notations of the semi-group and generator of the stochastic
process, for a given function $h\in\mathcal{D}^{k}(\mathcal{X})$, we denote by
$\Phi_{t}h(\boldsymbol{x}):=h(\phi_{t}(\boldsymbol{x}))$ the $C_{0}$-semi-
group corresponding to ODE (4). For differentiable $h$, the infinitesimal
generator is given by

(7) $\displaystyle\Lambda
h(\boldsymbol{x}):=D_{x}h(\boldsymbol{x})\bar{F}(\boldsymbol{x}).$

To strengthen intuition, the time derivative of $h(\phi_{t}(\boldsymbol{x}))$
can be expressed as

(8)
$\displaystyle\frac{d}{dt}h(\phi_{t}(\boldsymbol{x}))=D_{x}h(\phi_{t}(\boldsymbol{x}))\bar{F}(\phi_{t}(\boldsymbol{x}))=\Phi_{t}\Lambda
h(\boldsymbol{x}).$

By commuting $\Lambda$ and $\Phi_{t}$ we have
$\Lambda\Phi_{t}h(\boldsymbol{x})=D_{x}(h\circ\phi_{t})(\boldsymbol{x})\bar{F}(\boldsymbol{x})$
which is equal to (8) by the results of Appendix 1. We will use this property
to prove the theorems of the following section.

##### Abuse of Notations and Dependence on the Fast Component

The semi-group and generator of the stochastic system are generally defined
for functions $h$ in $\mathcal{D}^{k}(\mathcal{X}\times\mathcal{Y})$, i.e.,
functions which depend on the slow and fast component whereas the semi-group
and generator of the ODE are defined for $h$ in $\mathcal{D}^{k}(\mathcal{X})$
that do not depend on the fast component. What we refer to as abuse of
notation is the notation we use for the mapping $\psi_{t}$ of a function
$h\in\mathcal{D}^{k}(\mathcal{X})$. $\psi_{t}h$ still depends on the fast
component even if $h$ does not (since the evolution $\boldsymbol{X}_{t}$
depends on the state of $\boldsymbol{Y}_{t}$ and the initial values
$(\boldsymbol{x},\boldsymbol{y})$):

(9)
$\displaystyle\psi_{t}h(\boldsymbol{x},\boldsymbol{y})=\mathbb{E}[h(\boldsymbol{X}_{t})\mid\boldsymbol{X}_{0},\boldsymbol{Y}_{0}=\boldsymbol{x},\boldsymbol{y}],$

and

$\displaystyle Lh(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\sum_{\ell,\boldsymbol{y}^{\prime}}N\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})\bigl{(}h(\boldsymbol{x}+\ell/N)-h(\boldsymbol{x})\bigr{)}.$

In essence, this allows to use the notion of the semi-group and generator to
functions of the slow process. For consistency, we do the same for the ODE:
for an arbitrary $\boldsymbol{y}\in\mathcal{Y}$,
$\phi_{t}(\boldsymbol{x},\boldsymbol{y}):=\phi_{t}(\boldsymbol{x})$ and
therefore
$\Phi_{t}h(\boldsymbol{x},\boldsymbol{y})=h(\phi_{t}(\boldsymbol{x},\boldsymbol{y}))=h(\phi_{t}(\boldsymbol{x}))$
which is motivated by the abuse of notation in (9). This will merely be used
in the proofs and allows focusing on integral proof ideas instead of complex
notations.

###  4.2. Generator of the Fast Process and Regularity of the Poisson Equation

The generator $L$ describes the changes induced by the transitions of the fast
and slow process. In our analysis, it will be useful to analyze the changes
due to the jumps in $Y$ only. We denote by $L_{\text{fast}}$ the generator of
the Markov chain induced by $K(\boldsymbol{x})$ that takes as input a function
$h\in\mathcal{D}^{k}(\mathcal{X}\times\mathcal{Y})$ and associates another
function $L_{\text{fast}}h$ defined by:

$\displaystyle L_{\text{fast}}h(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\sum_{\boldsymbol{y}^{\prime}}K_{\boldsymbol{y},\boldsymbol{y}^{\prime}}(\boldsymbol{x})(h(\boldsymbol{x},\boldsymbol{y}^{\prime})-h(\boldsymbol{x},\boldsymbol{y})).$

Compared to (1), $L_{\text{fast}}$ is independent of $N$ because the
transition rates are rescaled by $1/N$.

Under Assumption ($A_{\text{2}}$), for all $\boldsymbol{x}$, the matrix
$K(\boldsymbol{x})$ characterizes a Markov chain on the state space
$\mathcal{Y}$ that has a unique stationary distribution denoted by
$\pi(\boldsymbol{x})$. $\pi_{y}(\boldsymbol{x})$ is the stationary probability
of $y\in\mathcal{Y}$ of the Markov chain induced by the kernel
$K(\boldsymbol{x})$. Subsequently, we will study the distance between the
’true’ stochastic process $\boldsymbol{Y}_{t}$ and an averaged system where
the distribution of $\boldsymbol{Y}_{t}$ is replaced by the stationary
distribution $\pi(\boldsymbol{x})$. To quantify the error made when replacing
$\boldsymbol{Y}_{t}$ with the stationary distribution, we consider the
following Poisson equation:

(10) $\displaystyle
h(\boldsymbol{x},\boldsymbol{y})-\sum_{y\in\mathcal{Y}}h(\boldsymbol{x},\boldsymbol{y})\pi_{y}(\boldsymbol{x})=L_{\text{fast}}G^{\text{fast}}_{h}(\boldsymbol{x},\boldsymbol{y}).$

For a given function $h:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}^{n}$ ($n$ is
arbitrary but finite), a function
$G^{\text{fast}}_{h}:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}^{n}$ that
satisfies the above equation is called a solution to this Poisson equation. We
will have particular interest in
$h(\boldsymbol{x},\boldsymbol{y})=F(\boldsymbol{x},\boldsymbol{y})\in\mathbb{R}^{d_{x}}$
namely when the drift $F(\boldsymbol{x},\boldsymbol{y})$ is compared to its
‘average’ version
$\bar{F}(\boldsymbol{x})=\sum_{y}F(\boldsymbol{x},\boldsymbol{y})\pi(\boldsymbol{y})$.

The existence of a regular solution to this Poisson equation is guaranteed by
the following Lemma 1. Note that the solution of the above Poisson equation is
not unique: If $G^{\text{fast}}_{h}$ is a solution, then for any constant
$c\in\mathbb{R}^{n}$, a function $G^{\text{fast}}_{h}+c$ is also a solution.
Later in the proofs of the theorems, when we talk about ‘a solution of the
Poisson equation’, we refer to the solution given in Lemma 1.

###### Lemma 1.

Assume ($A_{\text{2}}$). Then for all $\boldsymbol{x}\in\mathcal{X}$:

  1. (1)

The Markov chain corresponding to $K(\boldsymbol{x})$ has a unique stationary
distribution that we denote by $\pi(\boldsymbol{x})$. We denote by
$\Pi(\boldsymbol{x})=\mathds{1}\ \pi^{T}(\boldsymbol{x})$ the matrix where
each line is equal to $\pi(\boldsymbol{x})$.

  2. (2)

The matrix $(K(\boldsymbol{x})+\Pi(\boldsymbol{x}))$ is invertible and its
inverse is a generalized inverse of $K(\boldsymbol{x})$.

  3. (3)

Define
$K^{+}(\boldsymbol{x})=(K(\boldsymbol{x})+\Pi(\boldsymbol{x}))^{-1}(I-\Pi(\boldsymbol{x}))$,
then, for all functions $h:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}^{n}$,
$G_{h}(\boldsymbol{x},\boldsymbol{y})=\sum_{\boldsymbol{y}^{\prime}}K^{+}_{\boldsymbol{y},\boldsymbol{y}^{\prime}}(\boldsymbol{x})h(\boldsymbol{x},\boldsymbol{y}^{\prime})$
is a solution of the Poisson Equation (10).

If, in addition, Assumption ($A_{\text{1}}$) holds, then
$K^{+}(\boldsymbol{x})$ is twice differentiable in $\boldsymbol{x}$.

The proof is provided in Appendix B.1. Note in particular, this result implies
that if $h$ is (twice) differentiable in $\boldsymbol{x}$ then the same holds
true for $G^{\text{fast}}_{h}$.

###  4.3. Proof of Theorem 1 \- Transient State Proof

The proof of Theorem 1 can be decomposed in two main parts. We first use a
generator transformation to show that the slow system is well approximated by
a system whose drift is $D_{x}h(\boldsymbol{X}_{t})\
F(\boldsymbol{X}_{t},\boldsymbol{Y}_{t})$ (4.3.1). This leads us to treat
terms of the form
$\bar{F}(\boldsymbol{X}_{t})-F(\boldsymbol{X}_{t},\boldsymbol{Y}_{t})$. To
study them, we use the solution of the Poisson equation for the fast system
(10). The second part is the more technical and novel. It is detailed in
4.3.2. Some technical lemmas are postponed to Appendix B.2.

####  4.3.1. Error due to replacing the stochastic jumps of
$\boldsymbol{X}_{t}$ by the drift

For $h\in\mathcal{D}^{2}(\mathcal{X})$ recall that the $C_{0}$-semi-groups of
the stochastic system and of the ODE are defined as

$\displaystyle\psi_{s}h(\boldsymbol{x}_{0},\boldsymbol{y}_{0})=\mathbb{E}[h(\boldsymbol{X}_{s})\mid\boldsymbol{X}_{0},\boldsymbol{Y}_{0}=\boldsymbol{x}_{0},\boldsymbol{y}_{0}]\quad\text{
and
}\quad\Phi_{t}h(\boldsymbol{x}_{0},\boldsymbol{y}_{0})=h(\phi_{t}(\boldsymbol{x}_{0})).$

We define
$\nu_{s}h(\boldsymbol{x},\boldsymbol{y}):=\psi_{s}\Phi_{t-s}h(\boldsymbol{x},\boldsymbol{y})$
and rewrite

(11)
$\displaystyle\mathbb{E}[h(\boldsymbol{X}_{t})-h(\phi_{t}(\boldsymbol{x}))\mid\boldsymbol{X}_{0},\boldsymbol{Y}_{0}=\boldsymbol{x}_{0},\boldsymbol{y}_{0}]=\nu_{t}h(\boldsymbol{x}_{0},\boldsymbol{y}_{0})-\nu_{0}h(\boldsymbol{x}_{0},\boldsymbol{y}_{0})=\int_{0}^{t}\frac{d}{ds}\nu_{s}h(\boldsymbol{x}_{0},\boldsymbol{y}_{0})ds.$

To show that the last equation indeed holds true, observe that

$\displaystyle\frac{d}{ds}\nu_{s}h(\boldsymbol{x}_{0},\boldsymbol{y}_{0})$
$\displaystyle=\frac{d}{ds}\psi_{s}\Phi_{t-s}h(\boldsymbol{x}_{0},\boldsymbol{y}_{0})$
(12)
$\displaystyle=L\psi_{s}\Phi_{t-s}h(\boldsymbol{x}_{0},\boldsymbol{y}_{0})-\psi_{s}\Lambda\Phi_{t-s}h(\boldsymbol{x}_{0},\boldsymbol{y}_{0}).$

By regularity assumptions on the transition rates and bounded state space the
above equation is finite for all
$(\boldsymbol{x},\boldsymbol{y})\in\mathcal{X}\times\mathcal{Y}$ and time
$s\geq 0$. The dominated convergence theorem thus justifies the interchange of
derivative and integral and validates the last equality. As pointed out in
Appendix 1, $\psi_{s}$ and $L$, the stochastic semi-group and its
infinitesimal generator, commute, i.e., $L\psi_{s}=\psi_{s}L$. Hence,

$\displaystyle\eqref{eq:generator_comparison}$
$\displaystyle=\psi_{s}(L-\Lambda)\Phi_{t-s}h(\boldsymbol{x}_{0},\boldsymbol{y}_{0})$
(13)
$\displaystyle=\mathbb{E}[(L-\Lambda)\Phi_{t-s}h(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})\mid\boldsymbol{X}_{0},\boldsymbol{Y}_{0}=\boldsymbol{x}_{0},\boldsymbol{y}_{0}],$

where the last line follows by definition of $\psi_{s}$. Let
$g(\boldsymbol{x})=h(\phi_{t-s}(\boldsymbol{x}))$. $g$ is twice differentiable
with respect to the initial condition $\boldsymbol{x}$ as it lies in
$\mathcal{D}^{2}(\mathbb{R}^{d},\mathbb{R})$ by Lemma 3. By the use of Taylor
expansion rewrite $Lg$ as

$\displaystyle Lg(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\sum_{\ell,\boldsymbol{y}^{\prime}\in\mathcal{T}}N\alpha_{\ell,y^{\prime}}(\boldsymbol{x},\boldsymbol{y})(g(\boldsymbol{x}+\frac{\ell}{N})-g(\boldsymbol{x}))$
$\displaystyle=\sum_{\ell,\boldsymbol{y}^{\prime}\in\mathcal{T}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})D_{x}g(\boldsymbol{x})\ell+\frac{1}{N}\sum_{\ell,\boldsymbol{y}^{\prime}\in\mathcal{T}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})D^{2}_{x}g(\boldsymbol{x})\cdot(\ell,\ell)+o(\frac{\left\|g\right\|_{2,1}\sup_{\ell}\left|\ell\right|^{2}}{N})$
$\displaystyle=Dg(\boldsymbol{x})F(\boldsymbol{x},\boldsymbol{y})+\frac{1}{N}\sum_{\ell,\boldsymbol{y}^{\prime}\in\mathcal{T}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})D^{2}_{x}g(\boldsymbol{x})\cdot(\ell,\ell)+o(\frac{C_{R}\left\|g\right\|_{2,1}\sup_{\ell}\left|\ell\right|^{2}}{N}),$

with
$D^{2}g(\boldsymbol{x})\cdot(\ell,\ell)=\sum_{m,n}\frac{\partial^{2}g}{\partial
x_{m}\partial x_{n}}(\boldsymbol{x})\ell_{m}\ell_{n}$. The convergence depends
further on the Hölder norm of $h(\phi(x))$, bounds on the transition rates
$C_{R}$ and jump sizes $\sup\left|\ell\right|$. The generator of the ODE
related semi-group, $\Lambda$, introduced in Equation (7), maps $g$ to
$\Lambda
g(\boldsymbol{x},\boldsymbol{y})=D_{x}g(\boldsymbol{x})\bar{F}(\boldsymbol{x})$.
Hence, we have

$\displaystyle(L-\Lambda)g(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=D_{x}g(\boldsymbol{x})(F(\boldsymbol{x},\boldsymbol{y})-\bar{F}(\boldsymbol{x}))+\frac{1}{N}\sum_{\ell,y^{\prime}\in\mathcal{T}}\alpha_{\ell,y^{\prime}}(\boldsymbol{x},\boldsymbol{y})D^{2}_{x}g(\boldsymbol{x})\cdot(\ell,\ell)+o(1/N).$

####  4.3.2. Error due to replacing the drift by the average drift

Next, we have a closer look at the first summand of the right hand side in the
above equation. Denote by $G^{\text{fast}}_{F}$ a solution of the Poisson
equation (10) where the function “$h$’ of (10) is set to the drift $F$. Since
$D_{x}g(\boldsymbol{x})$ does not depend on $\boldsymbol{y}$ we have by
definition of the Poisson equation and $L_{\text{fast}}$:

(14) $\displaystyle
D_{x}g(\boldsymbol{x})(F(\boldsymbol{x},\boldsymbol{y})-\bar{F}(\boldsymbol{x}))=D_{x}g(\boldsymbol{x})L_{\text{fast}}G^{\text{fast}}_{F}(\boldsymbol{x},\boldsymbol{y})=L_{\text{fast}}D_{x}g(\boldsymbol{x})G^{\text{fast}}_{F}(\boldsymbol{x},\boldsymbol{y}).$

Combining the above with equation (13) and plugging everything into the
integral of equation (11), we get:

$\displaystyle\mathbb{E}[h(\boldsymbol{X}_{t})-h(\phi_{t}(\boldsymbol{x}_{0}))]=$
$\displaystyle\int_{0}^{t}\mathbb{E}[L_{\text{fast}}D_{x}(h\circ\phi_{t-s})(\boldsymbol{X}_{s})\cdot
G^{\text{fast}}_{F}(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})]$ (15)
$\displaystyle+\frac{1}{N}\mathbb{E}[\sum_{\ell,\boldsymbol{y}^{\prime}\in\mathcal{T}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})D^{2}_{x}(h\circ\phi_{t-s})(\boldsymbol{X}_{s})\cdot(\ell,\ell)]ds+o(1/N),$

where we suppress the conditioning on the initial values
$\boldsymbol{X}_{0},\boldsymbol{Y}_{0}=\boldsymbol{x}_{0},\boldsymbol{y}_{0}$.
Let
$H_{s}(\boldsymbol{x},\boldsymbol{y}):=\mathbb{E}[D_{x}(h\circ\phi_{t-s})(X_{s})G^{\text{fast}}_{F}(X_{s},Y_{s})\mid\boldsymbol{X}_{0},\boldsymbol{Y}_{0}=\boldsymbol{x},\boldsymbol{y}]$.
Applying Lemma 1 with
$g_{s}(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})=D_{x}(h\circ\phi_{t-s})(X_{s})G^{\text{fast}}_{F}(X_{s},Y_{s})$
implies that

$\displaystyle 0=\frac{1}{N}H_{t}(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle-\frac{1}{N}H_{0}(\boldsymbol{x},\boldsymbol{y})-\int_{0}^{t}\mathbb{E}[\frac{1}{N}D_{x}\bigl{(}D_{x}h(\phi_{t-s}(\boldsymbol{X}_{s}))\cdot\bar{F}(\phi_{t-s}(\boldsymbol{X}_{s}))\bigr{)}\cdot
G^{\text{fast}}_{F}(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})]ds$ (16)
$\displaystyle\phantom{=\int}+\mathbb{E}[\frac{1}{N}LD_{x}(h\circ\phi_{t-s})(X_{s})G^{\text{fast}}_{F}(X_{s},Y_{s})]ds.$

As $h$ and $\phi_{t-s}$ are twice continuously differentiable, by compactness
of $\mathcal{X}$, and as $G^{\text{fast}}_{F}$ is finite,
$D_{x}\bigl{(}D_{x}h(\phi_{t-s}(\boldsymbol{X}_{s}))\bar{F}(\phi_{t-s}(\boldsymbol{X}_{s}))\bigr{)}G^{\text{fast}}_{F}(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})$
is bounded. Moreover, by Lemma 2 we have

$\displaystyle(L_{\text{fast}}-\frac{1}{N}L)D_{x}(h\circ\phi_{t-s})(X_{s})G^{\text{fast}}_{F}(X_{s},Y_{s})$
$\displaystyle=\frac{1}{N}\sum_{\ell,\boldsymbol{y}^{\prime}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})\ell\
D_{x}(D_{x}(h\circ\phi_{t-s})(\boldsymbol{X}_{s})G^{\text{fast}}_{F}(\boldsymbol{X}_{s},\boldsymbol{y}^{\prime}))+o(1/N).$

Define

$\displaystyle
C_{h}(t)=\int_{0}^{t}\mathbb{E}[\sum_{\ell,\boldsymbol{y}^{\prime}\in\mathcal{T}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})D^{2}_{x}h(\phi_{t-s}(\boldsymbol{X}_{s}))\cdot(\ell,\ell)]ds$
$\displaystyle+\int_{0}^{t}\mathbb{E}[D_{x}\bigl{(}D_{x}h(\phi_{t-s}(\boldsymbol{X}_{s}))\cdot\bar{F}(\phi_{t-s}(\boldsymbol{X}_{s}))\bigr{)}\cdot
G^{\text{fast}}_{F}(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})]ds$
$\displaystyle+\int_{0}^{t}\mathbb{E}[\sum_{\ell,\boldsymbol{y}^{\prime}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{s},\boldsymbol{Y}_{s})\ell\
D_{x}(D_{x}(h\circ\phi_{t-s})(\boldsymbol{X}_{s})G^{\text{fast}}_{F}(\boldsymbol{X}_{s},\boldsymbol{y}^{\prime}))]ds+\mathbb{E}[H_{t}(\boldsymbol{X}_{t},\boldsymbol{Y}_{t})]-H_{0}(\boldsymbol{x},\boldsymbol{y}),$

which lets us rewrite

$\displaystyle
N\left(\mathbb{E}[h(\boldsymbol{X}_{t})]-h(\phi_{t}(\boldsymbol{x}))\right)=C_{h}(t)+o(1).$

Taking the limit $N\to\infty$ concludes the proof.

###  4.4. Proof of Theorem 2 \- Steady-State Proof

Considering
$\mathbb{E}[h(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})-\sum_{\boldsymbol{y}}\pi_{y}(\phi_{\text{\tiny$\infty$}})h(\phi_{\text{\tiny$\infty$}},\boldsymbol{y})]$,
first, add and subtract the term
$\mathbb{E}[\sum_{\boldsymbol{y}}\pi_{y}(\boldsymbol{X}_{\infty})h(\boldsymbol{X}_{\infty},\boldsymbol{y})]$
into the equation, yielding

(A)
$\displaystyle\mathbb{E}[h(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})-\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\phi_{\text{\tiny$\infty$}})h(\phi_{\text{\tiny$\infty$}},\boldsymbol{y})]$
$\displaystyle=\mathbb{E}[h(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})-\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\boldsymbol{X}_{\infty})h(\boldsymbol{X}_{\infty},\boldsymbol{y})]$
(B)
$\displaystyle\phantom{=}+\mathbb{E}[\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\boldsymbol{X}_{\infty})h(\boldsymbol{X}_{\infty},\boldsymbol{y})-\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\phi_{\text{\tiny$\infty$}})h(\phi_{\text{\tiny$\infty$}},\boldsymbol{y})].$

Treating the two terms (A) and (B) separately, we define $G^{\text{fast}}_{h}$
and $G^{\text{slow}}_{h}$ as the solutions to the Poisson equations

(‘fast’ Poisson Equation) $\displaystyle
L_{\text{fast}}G^{\text{fast}}_{h}(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=h(\boldsymbol{x},\boldsymbol{y})-\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\boldsymbol{x})h(\boldsymbol{x},\boldsymbol{y}),$
(‘slow’ Poisson Equation) $\displaystyle\Lambda
G^{\text{slow}}_{h}(\boldsymbol{x})$
$\displaystyle=\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\boldsymbol{x})h(\boldsymbol{x},\boldsymbol{y})-\sum_{\boldsymbol{y}}\pi_{\boldsymbol{y}}(\phi_{\text{\tiny$\infty$}})h(\phi_{\text{\tiny$\infty$}},\boldsymbol{y}).$

Recall that the existence and regularity properties of the function
$G^{\text{fast}}_{h}$ are investigated in Lemma 1. For $G^{\text{slow}}_{h}$,
it is known (e.g., (Gast, 2017)) that the exponential stability of the unique
fixed point ($A_{\text{3}}$) along with the smoothness of $F$ guarantees that
$G^{\text{slow}}_{h}\in D^{2}(\mathcal{X})$.

Next we use that taking the expectation of the generator $L$ applied to an
arbitrary function $g$ over the stationary distribution of
$(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})$ is zero, i.e.,
$\mathbb{E}[Lg(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})]=0.$ By using
Lemma 2, it follows

$\displaystyle\eqref{eq:A}$
$\displaystyle=\mathbb{E}[L_{\text{fast}}G^{\text{fast}}_{h}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})-\frac{1}{N}LG^{\text{fast}}_{h}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})]$
$\displaystyle=\mathbb{E}[\frac{1}{N}\sum_{\ell,\boldsymbol{y}^{\prime}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})\ell\
D_{x}G^{\text{fast}}_{h}(\boldsymbol{X}_{\infty},\boldsymbol{y}^{\prime})+o(1/N)].$

For the second term, we apply the ‘slow’ Poisson Equation equation and
subtract
$\mathbb{E}[LG^{\text{slow}}_{h}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})]=0$
which yields

(17) $\displaystyle\eqref{eq:B}=\mathbb{E}[\Lambda
G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})-LG^{\text{slow}}_{h}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})].$

Note that even if $G^{\text{slow}}_{h}$ only depends on $\boldsymbol{x}$, the
function $LG^{\text{slow}}_{h}$ does depend on $\boldsymbol{x}$ and
$\boldsymbol{y}$ because the rate functions that appear in the generator $L$
do depend on $\boldsymbol{x}$ and $\boldsymbol{y}$. By using the Taylor
expansion, $LG^{\text{slow}}_{h}$ equals

$\displaystyle LG^{\text{slow}}_{h}(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\sum_{\ell,\boldsymbol{y}^{\prime}}N\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})\bigl{(}G^{\text{slow}}_{h}(x+\frac{\ell}{N})-G^{\text{slow}}_{h}(\boldsymbol{x})\bigr{)}$
$\displaystyle=\sum_{\ell,\boldsymbol{y}^{\prime}}N\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})D_{x}G^{\text{slow}}_{h}(\boldsymbol{x})\frac{\ell}{N}+\frac{1}{2}\sum_{\ell,\boldsymbol{y}^{\prime}}N\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})D_{x}^{2}G^{\text{slow}}_{h}(\boldsymbol{x})\frac{(\ell,\ell)}{N^{2}}+o(\frac{\left\|\ell\right\|^{2}}{N}).$

By application of the definition of $F(\boldsymbol{x},\boldsymbol{y})$, the
first term of the right hand side is equal to
$F(\boldsymbol{x},\boldsymbol{y})D_{x}G^{\text{slow}}_{h}(\boldsymbol{x})$.
Moreover, by definition, $\Lambda
G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})=D_{x}G^{\text{slow}}_{h}(\boldsymbol{x})\bar{F}(\boldsymbol{x})$.
This shows that

$\displaystyle\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})(\bar{F}(\boldsymbol{X}_{\infty})-F(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty}))+\sum_{\ell,\boldsymbol{y}^{\prime}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})\bigl{(}\frac{1}{2}D_{x}^{2}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})\cdot\frac{(\ell,\ell)}{N}\bigr{)}+o(\frac{\left\|\ell\right\|}{N^{2}})].$

As $G^{\text{slow}}_{h}$ has a bounded second derivative, the second term is
of order $O(1/N)$. However, to bound the first summand we need to investigate
$\bar{F}(\boldsymbol{x})-F(\boldsymbol{x},\boldsymbol{y})=\sum_{y}\pi_{\boldsymbol{y}^{\prime}}(\boldsymbol{x})F(\boldsymbol{x},\boldsymbol{y}^{\prime})-F(\boldsymbol{x},\boldsymbol{y})$.
Falling back on (‘fast’ Poisson Equation) allows two write

$\displaystyle
D_{x}G^{\text{slow}}_{h}(\boldsymbol{x})(\bar{F}(\boldsymbol{x})-F(\boldsymbol{x},\boldsymbol{y}))=-D_{x}G^{\text{slow}}_{h}(\boldsymbol{x})L_{\text{fast}}G^{\text{fast}}_{F}(\boldsymbol{x},\boldsymbol{y})=-L_{\text{fast}}D_{x}G^{\text{slow}}_{h}(\boldsymbol{x})G^{\text{fast}}_{F}(\boldsymbol{x},\boldsymbol{y}).$

The last equality holds due to the definition of $L_{\text{fast}}$, and
$D_{x}G^{\text{slow}}_{h}(\boldsymbol{x})$ which depends only on
$\boldsymbol{x}$. By adding the zero term
$\mathbb{E}[\frac{1}{N}LD_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})G^{\text{fast}}_{F}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})]$
and using Lemma 2. We see that

$\displaystyle\mathbb{E}[-L_{\text{fast}}D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})G^{\text{fast}}_{F}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})+\frac{1}{N}LD_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})G^{\text{fast}}_{F}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})]$
$\displaystyle=\mathbb{E}[\frac{1}{N}\sum_{\ell,\boldsymbol{y}^{\prime}}\alpha_{\ell,y^{\prime}}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})\ell\
D_{x}(D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})G^{\text{fast}}_{F}(\boldsymbol{X}_{\infty},\boldsymbol{y}^{\prime}))+o(1/N)].$

Put together, the error term of (B) is

(18) $\displaystyle\mathbb{E}[\frac{1}{N}\sum_{\ell,\boldsymbol{y}^{\prime}}$
$\displaystyle\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})\ell\
D_{x}(D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})G^{\text{fast}}_{F}(\boldsymbol{X}_{\infty},\boldsymbol{y}^{\prime}))+\frac{1}{2N}\sum_{\ell,\boldsymbol{y}^{\prime}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})D_{x}^{2}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})(\ell,\ell)].$

By defining $C_{h}$ as the sum of the residual term
$\mathbb{E}[\sum_{\ell,\boldsymbol{y}^{\prime}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})\ell\
D_{x}G^{\text{fast}}_{h}(\boldsymbol{X}_{\infty},\boldsymbol{y}^{\prime})]$
for (A) and equation (18) the by $N$ scaled, the proof concludes.

###  4.5. Proof of Theorem 3 (Refinement Theorem, and Closed form Expressions)

In this section we first decompose the constant $C_{h}$ in two terms in
Proposition 2. We then study these two terms in Proposition 3 and 4 where we
obtain the closed form expressions for the correction terms which allow to
numerically obtain the corrections.

We define two functions $Q(x,y)$ and $J_{h}(x,y)$ (the latter is defined for a
function $h\in\mathcal{D}^{2}(\mathcal{X})$):

$\displaystyle Q(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\sum_{\ell,\boldsymbol{y}^{\prime}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})\ell\ell^{T},$
$\displaystyle J_{h}(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=-\sum_{\ell,\boldsymbol{y}^{\prime}}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})D_{x}\bigl{(}D_{x}G^{\text{slow}}_{h}(\boldsymbol{x})G^{\text{fast}}_{F}(\boldsymbol{x},\boldsymbol{y}^{\prime})\bigr{)}\ell,$

and we denote by $\bar{J_{h}}(x)=\sum_{y}\pi_{y}(x)J(x,y)$ and
$\bar{Q}(x)=\sum_{y}\pi_{y}(x)Q(x,y)$ their “average” version. The following
proposition holds.

###### Proposition 2.

Assume ($A_{\text{1}}$)–($A_{\text{3}}$). Then, for any
$h\in\mathcal{D}^{2}(\mathcal{X})$, we have:

(19) $\displaystyle
N\bigl{(}\mathbb{E}[h(\boldsymbol{X}_{\infty})]-h(\phi_{\text{\tiny$\infty$}})\bigr{)}=\frac{1}{2}D_{x}^{2}G^{\text{slow}}_{h}(\phi_{\text{\tiny$\infty$}})\bar{Q}(\phi_{\text{\tiny$\infty$}})+\bar{J_{h}}(\phi_{\text{\tiny$\infty$}})+o(1).$

###### Proof.

Let $h\in\mathcal{D}^{2}(\mathcal{X})$ and let us denote $G^{\text{slow}}_{h}$
the solution of (‘slow’ Poisson Equation). As $h$ does not depend on $y$, this
function is such that for any $x$:
$h(x)-h(\phi_{\text{\tiny$\infty$}})=\Lambda
G^{\text{slow}}_{h}(x)=D_{x}G^{\text{slow}}_{h}(x)\bar{F}(x)$. The following
steps are similar to the ones for (B) in the proof of Theorem 2. Hence,
applying this for $\boldsymbol{X}_{\infty}$ and taking the expectation, we
get:

$\displaystyle
N\mathbb{E}[h(\boldsymbol{X}_{\infty})-h(\phi_{\text{\tiny$\infty$}})]=N\mathbb{E}[\Lambda
G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})]=N\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})\bar{F}(\boldsymbol{X}_{\infty})].$

Recall that for any bounded function
$g\in\mathcal{D}(\mathcal{X}\times\mathcal{Y})$,
$\mathbb{E}[Lg(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})]=0$. Hence,

$\displaystyle
N\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})\bar{F}(\boldsymbol{X}_{\infty})-LG^{\text{slow}}_{h}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})].$

Similarly to what we do to prove Theorem 2, we plug the definition of $L$ in
the above equation and use a Taylor expansion to show that this equals

$\displaystyle
N\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})\bar{F}(\boldsymbol{X}_{\infty})-\sum_{\ell,\boldsymbol{y}^{\prime}}N\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})\bigl{(}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty}+\frac{\ell}{N})-G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})\bigr{)}]$
$\displaystyle=N\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})\bar{F}(\boldsymbol{X}_{\infty})-\sum_{\ell,\boldsymbol{y}^{\prime}}N\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})\bigl{(}D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})\frac{\ell}{N}+\frac{1}{2}D_{x}^{2}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})\frac{(\ell,\ell)}{N^{2}}+o(\frac{\left\|\ell\right\|}{N^{2}})\bigr{)}]$
(20)
$\displaystyle=N\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})(\bar{F}(\boldsymbol{X}_{\infty})-F(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty}))-\frac{1}{2N}D_{x}^{2}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})Q(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})]+o(1).$

The rest of the proof follows using

  * •

Lemma 3 to show
$N\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})(\bar{F}(\boldsymbol{X}_{\infty})-F(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty}))]=\bar{J}(\phi_{\text{\tiny$\infty$}})+o(1)$;

  * •

Theorem 2 which states
$\mathbb{E}[\frac{1}{2}D_{x}^{2}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})Q(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})]=\frac{1}{2}D_{x}^{2}G^{\text{slow}}_{h}(\phi_{\text{\tiny$\infty$}})\bar{Q}(\phi_{\text{\tiny$\infty$}})+o(1)$.

Applying the above equations to (20) implies the statement of Proposition 2. ∎

In the rest of the section, we show that the quantity
$\frac{1}{2}D_{x}^{2}G^{\text{slow}}_{h}(\phi_{\text{\tiny$\infty$}})\bar{Q}(\phi_{\text{\tiny$\infty$}})+\bar{J}(\phi_{\text{\tiny$\infty$}})$
can be easily computed numerically by solving linear systems of equations. As
shown in Proposition 3 and 4, we obtain five different correction terms:

  * •

The first two $V\in\mathbb{R}^{d_{x}},W\in\mathbb{R}^{d_{x}\times d_{x}}$ are
closely related to the ones obtained in (Gast and Van Houdt, 2017)[Theorem
3.1] and derived from the term
$\frac{1}{2}D_{x}^{2}G^{\text{slow}}_{h}(\phi_{\text{\tiny$\infty$}})\bar{Q}(\phi_{\text{\tiny$\infty$}})$.
This term is essentially identical to the refinement term of (Gast and Van
Houdt, 2017).

  * •

In Proposition 4 we derive three other refinement terms
$S,T\in\mathbb{R}^{d_{x}}$ and $U\in\mathbb{R}^{d_{x}\times d_{x}}$ which give
closed form descriptions of $\bar{J}(\phi_{\text{\tiny$\infty$}})$. These
terms are novel and take into account the difference of the average drift
$\bar{F}$ and the actual drift $F$ obtained from the stochastic system. The
proof is based on exact expressions for the Poisson equations
$G^{\text{fast}},G^{\text{slow}}$ and relies on Lemma 1.

####  4.5.1. Correction terms $V$ and $W$

To obtain the first result, let $A$ and $B$ be the Jacobian and Hessian matrix
of the ‘average’ drift, i.e.,

(21) $\displaystyle A_{i,j}=\frac{\partial\bar{F}_{i}}{\partial
x_{j}}(\phi_{\text{\tiny$\infty$}})\quad\text{ and }\quad
B_{i,k_{1},k_{2}}=\frac{\partial^{2}\bar{F}_{i}}{\partial
x_{k_{1}}x_{k_{2}}}(\phi_{\text{\tiny$\infty$}}).$

By the exponential stability of the fixed point (Assumption ($A_{\text{3}}$)),
the matrix $A$ is invertible and the Lyapunov equation $AW+WA^{T}+\bar{Q}=0$
has a unique solution. We denote by $W$ its solution and define the vector $V$
as

$\displaystyle
V_{i}=-\frac{1}{2}\sum_{j}(A^{-1})_{i,j}\sum_{k_{1},k_{2}}B_{j,k_{1},k_{2}}W_{k_{1},k_{2}}.$

###### Proposition 3.

Assume ($A_{\text{1}}$) \- ($A_{\text{3}}$) and
$h\in\mathcal{D}^{2}(\mathcal{X})$. Then:

$\displaystyle\frac{1}{2}D_{x}^{2}G^{\text{slow}}_{h}(\phi_{\text{\tiny$\infty$}})\bar{Q}(\phi_{\text{\tiny$\infty$}})=\sum_{i}\frac{\partial
h}{\partial x_{i}}V_{i}+\frac{1}{2}\sum_{i,j}\frac{\partial^{2}h}{\partial
x_{i}x_{j}}W_{i,j}\ .$

###### Proof.

The proposition is a direct consequence of the results of (Gast and Van Houdt,
2017)[Theorem 3.1] which we apply to the function
$\bar{F}(x)=\sum_{y}\pi_{y}(x)F(x,y)$. ∎

####  4.5.2. Correction terms $T$, $S$, $U$

###### Proposition 4.

Assume ($A_{\text{1}}$) \- ($A_{\text{3}}$) and
$h\in\mathcal{D}^{2}(\mathcal{X})$. The closed form solution of $\bar{J}$ at
the equilibrium point $\phi_{\text{\tiny$\infty$}}$ is given by

$\displaystyle\bar{J}_{h}(\phi_{\text{\tiny$\infty$}})=\sum_{i}\frac{\partial
h}{\partial
x_{i}}(\phi_{\text{\tiny$\infty$}})(T_{i}+S_{i})+\sum_{i,j}\frac{\partial^{2}h}{\partial
x_{i}\partial x_{j}}(\phi_{\text{\tiny$\infty$}})U_{i,j}.$

$U$ is the unique solution to the Sylvester equation

$\displaystyle AX+XA^{T}=-O,\qquad\text{ with }$ $\displaystyle
O=\sum_{y}\pi_{y}(\phi_{\text{\tiny$\infty$}})\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\phi_{\text{\tiny$\infty$}},y)K^{+}_{y^{\prime},:}(\phi_{\text{\tiny$\infty$}})F(\phi_{\text{\tiny$\infty$}},:)\ell^{T},$

where $A,B$ are the Jacobian, Hessian of the average drift as defined in (21)
and $K^{+}$ as given in Lemma 1. $T,S$ are defined by

$\displaystyle
T_{i}:=\sum_{j}A^{-1}_{i,j}\sum_{k_{1},k_{2}}B_{j,k_{1},k_{2}}U_{k_{1},k_{2}},$
$\displaystyle
S_{i}:=\sum_{k}A^{-1}_{i,k}\sum_{y}\pi_{y}(\phi_{\text{\tiny$\infty$}})\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\phi_{\text{\tiny$\infty$}},y)\left(\sum_{y^{\prime}}F_{k}(\phi_{\text{\tiny$\infty$}},y^{\prime})\nabla^{T}_{x}K_{y^{\prime},y^{\prime}}^{+}(\phi_{\text{\tiny$\infty$}})\ell+K_{y^{\prime},y^{\prime}}^{+}(\phi_{\text{\tiny$\infty$}})\nabla^{T}_{x}F_{k}(\phi_{\text{\tiny$\infty$}},y^{\prime})\ell\right).$

The proof of this proposition is given in Appendix B.5.

##  5\. Example: CSMA Model

To illustrate our results, we consider the unsaturated CSMA random-access
networks studied in (Cecchi et al., 2021). In this paper, the authors use a
two-scale model to study the performance of a CSMA algorithm with many nodes.
The slow process corresponds to the arrival of jobs and the fast process
corresponds to the activation and deactivation of nodes. The authors of this
paper derive a mean field approximation and show that it is asymptotically
exact. With our methods, we go two steps further:

  * •

Theorem 1 and 2 show that not only the mean field approximation is
asymptotically exact but also that the error is only of order $O(1/N)$.

  * •

By using Theorem 3, we can compute a refinement term. Our numerical example
shows that, similarly to what happens for classical one-scale mean field
models (Gast and Van Houdt, 2017), this refinement term is extremely accurate.
It is much more accurate than the classical mean field approximation when the
studied system is not too large.

$1$$2$$3$ (a) The linear 3 Node Graph.

$3$$1$$5$$4$$2$ (b) The 5 Node Graph.

Figure 1. Two examples of interference graphs.

###  5.1. Model Description

We consider a model with $C$ server types, with $N$ statistically identical
servers for each class. All servers communicate through a wireless medium
using a random-access protocol and have a finite buffer size $B$. The $C$
classes form an interference graph $G=(\mathcal{C},\mathcal{E})$ with
$\mathcal{C}=\\{1,\dots,C\\}$ the classes and $\mathcal{E}$ the network
specific edges. (see for instance Figure 1 for examples of graphs with three
or five classes). This interference graph indicates that two servers cannot
transmit simultaneously if they either are of the same class or belong to an
adjacent class. For each class $c\in\mathcal{C}$, we will denote by
$y_{c}\in\\{0,1\\}$ a variable that equals $1$ if a node of class $c$ is
transmitting and $0$ otherwise. We denote by $\mathcal{Y}$ the set of possible
activation vectors $y\in\\{0,1\\}^{C}$. It is equal to the set of the
independent sets of the graph, i.e., all activity vectors for which an active
node has only inactive neighbors. For instance, for the graph with three
classes which are linearly connected as shown in Figure 1(a), the set of
feasible states is given by

$\displaystyle\mathcal{Y}=\\{(0,0,0),(1,0,0),(0,1,0),(0,0,1),(1,0,1)\\}.$

Any node can turn active if there are no neighboring active nodes. Once its
transmission is finished, it changes back to an idle state. For a given class
$c$, we define $\mathcal{Y}_{c}^{+}:=\\{y\in\mathcal{Y}:y_{c}=1\\}$ as the
subset of states for which node $c$ is active, and by
$\mathcal{Y}_{c}^{-}:=\\{y\in\mathcal{Y}:y_{c}=y_{d}=0\ \forall d\text{ s.t.
}(c,d)\in\mathcal{E}\\}$ the subset of states from which node $c$ can turn
active, i.e., all neighboring nodes are inactive. For instance, for the linear
$3$ node graph of Figure 1(a) and the class $c=1$ this yields the following
sets

$\displaystyle\mathcal{Y}_{1}^{+}=\\{(1,0,0),(1,0,1)\\},$
$\displaystyle\mathcal{Y}_{1}^{-}=\\{(0,0,0),(0,0,1)\\}.$

As in (Cecchi et al., 2021), in this non-saturated model, we consider that if
a node is in class $c$, new packets arrive to this node at rate
$\lambda_{c}>0$. If a node has a packet to transmit and no neighboring node is
transmitting, then this node becomes active at rate $\nu_{c}>0$. We assume
that a transmission from a node of class $c$ takes an exponential time of
duration $1/(N\mu_{c})>0$, after which the packet leaves the system.

###  5.2. Two-scale model representation

The model as described above fits in our two-timescale representation. To see
why, for each class $c\in\mathcal{C}$ and buffer size $b\in\\{0\dots B\\}$, we
define $X^{(N)}_{t,(c,b)}$ as the fraction of servers of class $c$ that have
at least $b$ jobs in their queue at time $t$. We denote by
$\boldsymbol{X}_{t}$ the vector of all possible $X_{t,(c,b)}$ for all
$c\in\mathcal{C}$ and $b\in\\{0\dots B\\}$. The fast component
$\boldsymbol{Y}_{t}$ is the activation at time $t$: $Y_{t,(c)}=1$ if a node of
class $c$ is transmitting at time $t$ and $0$ otherwise.

Using this representation, we characterize the possible transitions. Given a
state pair $(\boldsymbol{x},\boldsymbol{y})$, the transitions are represented
as a transition vector of the form $(\ell/N,\boldsymbol{y}^{\prime})$ and a
corresponding transition rate
$N\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})$ such
that the state $(\boldsymbol{x},\boldsymbol{y})$ jumps to
$(\boldsymbol{x}+\ell/N,\boldsymbol{y}^{\prime})$ at rate
$N\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})$. The
transitions can be distinguished into three types:

  * •

Arrival of a packet to a server of class $c\in\mathcal{C}$:

$\displaystyle\quad{\begin{array}[]{lcll}(\frac{e_{c,i}}{N},\boldsymbol{y})&\text{at
rate}&N\lambda_{c}(1-x_{c,i})&\text{for }i=1,\\\
(\frac{e_{c,i}}{N},\boldsymbol{y})&\text{at
rate}&N\lambda_{c}(x_{c,i-1}-x_{c,i})&\text{for }2\leq i\leq B.\end{array}}$

  * •

Back-off of a server of class $c\in\mathcal{C}$ with at least one packet if
the class activity vector allows the back-off, i.e.,
$\boldsymbol{y}\in\mathcal{Y}_{c}^{-}$:

$\displaystyle\quad{\begin{array}[]{lcll}(-\frac{e_{c,i}}{N},\boldsymbol{y}+e_{c})&\text{at
rate}&N\nu_{c}(x_{c,i}-x_{c,i+1})&\text{for }1\leq i<B,\\\
(-\frac{e_{c,i}}{N},\boldsymbol{y}+e_{c})&\text{at
rate}&N\nu_{c}x_{c,i}&\text{for }i=B.\\\ \end{array}}$

  * •

Transmission completion of an active node of class $c\in\mathcal{C}$, i.e.,
$\boldsymbol{y}\in\mathcal{Y}^{+}_{c}$:

$\displaystyle\quad{\begin{array}[]{lcll}(0,\boldsymbol{y}-e_{c})&\text{at
rate}&N\mu_{c}.&\\\ \end{array}}$

We denote the set of all possible transition vectors
$(\ell/N,\boldsymbol{y}^{\prime})$ from a state pair
$(\boldsymbol{x},\boldsymbol{y})$ by
$\mathcal{T}(\boldsymbol{x},\boldsymbol{y})$.

###  5.3. Steady-State distribution $\pi_{y}(x)$ and average drift

By using the above transition definitions, the matrix $K$ is given by:

  * •

If $\boldsymbol{y}\in\mathcal{Y}^{-}_{c}$, then
$K(\boldsymbol{y},\boldsymbol{y}+e_{c})=N\nu_{c}x_{c,1}$,

  * •

If $\boldsymbol{y}\in\mathcal{Y}^{+}_{c}$, then
$K(\boldsymbol{y},\boldsymbol{y}-e_{c})=N\mu_{c}$,

all other entries of the matrix being $0$.

This representation is used by the authors of (Cecchi et al., 2021) to derive
the product-form stationary distribution for a fixed server state
$\boldsymbol{x}$. This product form is closely related to the product-form
stationary distribution of saturated networks as found in (van de Ven et al.,
2010; Wang and Kar, 2005; Boorstyn et al., 1987): The quantity
$\pi_{\boldsymbol{y}}(x)$ is calculated as follows:

$\displaystyle\pi_{\boldsymbol{y}}(x)$
$\displaystyle:=\frac{Z(\boldsymbol{x},\boldsymbol{y})}{Z(\boldsymbol{x})},$
with $\displaystyle Z(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\prod_{c\in\mathcal{C}}\bigl{(}\frac{\nu_{c}}{\mu_{c}}x_{c,1}\bigr{)}^{y_{c}},$
$\displaystyle Z(\boldsymbol{x})$
$\displaystyle=\sum_{\boldsymbol{y}\in\mathcal{Y}}Z(\boldsymbol{x},\boldsymbol{y}).$

Following our definition of Section 2.2 the drift and its average version are
generically defined by:

$\displaystyle F(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\sum_{(\ell,\boldsymbol{y}^{\prime})\in\mathcal{T}(\boldsymbol{x},\boldsymbol{y})}\alpha_{\ell,\boldsymbol{y}^{\prime}}(\boldsymbol{x},\boldsymbol{y})\ell,$
and $\displaystyle\bar{F}(\boldsymbol{x})$
$\displaystyle=\sum_{w\in\mathcal{Y}}\pi_{\boldsymbol{y}}(\boldsymbol{x})F(\boldsymbol{x},\boldsymbol{y}).$

For the drift $F(\boldsymbol{x},\boldsymbol{y})$ of the random access model
this leads to the closed form expression

$\displaystyle F(\boldsymbol{x},\boldsymbol{y})$
$\displaystyle=\sum_{c\in\mathcal{C}}\bigl{(}\sum_{1<i\leq B}e_{c,i}\
\lambda_{c}(x_{c,i-1}-x_{c,i})+e_{c,1}\ \lambda_{c}(1-x_{c,1})$
$\displaystyle\qquad\quad-\mathbf{1}_{\\{y\in\mathcal{Y}_{c}^{-}\\}}(\sum_{1\leq
i<B}e_{c,i}\ \nu_{c}(x_{c,i}-x_{c,i+1})+e_{c,B}\ \nu_{c}x_{c,B})\bigr{)}.$

It should be clear that assumption ($A_{\text{1}}$) holds in our case because
the rates given in Section 5.2 are all continuous in $\boldsymbol{x}$ (in fact
they are all linear). Moreover, the model also satisfies Assumption
($A_{\text{2}}$): For a given $\boldsymbol{x}$, the set of irreducible states
for $K(\boldsymbol{x})$ contains all the feasible activation vectors $y$ such
that $y_{c}=0$ if $x_{(c,1)}=0$. The condition $x_{{c,1}}=0$ implies that the
nodes of class $c$ do not have any packets to transmit. The situation of
Assumption ($A_{\text{3}}$) is more complicated. To the best of our knowledge,
there has not been a complete stability characterization for the unsaturated
random access CSMA model. Cecchi et al. show in (Cecchi et al., 2021) that in
the case of a complete interference graph stability conditions can be derived
which assure global exponential stability. They further conjectured that
similar results hold for general interference graphs. In our analysis, we
assume that ($A_{\text{3}}$) holds.

###  5.4. Numerical results

To study the accuracy of the mean field approximation and the refined term
proposed in Theorem 2, we implemented a Python
library222https://gitlab.inria.fr/sallmeie/bias-and-refinement-of-multiscale-
mean-field-models to simulate the system and compute the mean field
approximation and the refinement term. This library is generic and can take as
an input any instance of the model which we defined above. For instance, in
the Code Cell LABEL:lst:python_code, we illustrate how to use this library to
construct a model where the interference graph is as in Figure 1(b), the rates
are $\lambda=[.5,.7,.7,.6,.4],\nu=[4,3,3,3,3],\mu=[3,3,2,4,2]$, and the buffer
size is equal to $10$. This cell shows how to initialize the 5 node model and
obtain the approximation and refinement from our implementation. We also
perform the same experiments with the linear 3 node model, for which we
provide the results in Appendix D. Note that the results are qualitatively
very similar.

⬇

# Graph structure (this is the five node example)

G = np.array([[0,1,1,1,0],

[1,0,0,0,0],

[1,0,0,0,1],

[1,0,0,0,1],

[0,0,1,1,0]])

# rates & buffer size

_lambda = np.array([.5,.7,.7,.6,.4])

nu = np.array([4,3,3,3,3])

mu = np.array([3,3,2,4,2])

buffer_size = 10

# We define the model, compute a trajectory and the refinement term.

csma = symbolic_CSMA(nu, mu, _lambda, G, buffer_size)

T, X = csma.ode(time=200) # mean field ODE

v, s, t, w, u = csma.compute_refinements(X[-1]) # steady-state refinement

In order to compute the refinement terms, the library needs to compute various
derivatives (of the drift or of the matrix $K^{+}(\boldsymbol{x})$). To
implement this, we rely on symbolic differentiation provided by the sympy
library (Meurer et al., 2017). As we see later in Table 1, the use of the
symbolic differentiation is the performance bottleneck of our implementation.
In Appendix C we furthermore show how to obtain the stationary distribution
and the derivative of $K^{+}(x)$ numerically.

##### Transient regime and illustration of Theorem 1

To illustrate the accuracy of the mean field, we use the 5 node model
described in the Code Cell LABEL:lst:python_code. We first simulate the CSMA
model and compare it with the mean field ODE. The results are reported in
Figure 2 where we plot the mean field approximation against a sample mean
$\mathbb{E}[\boldsymbol{X}_{t}]$ derived from $1000$ simulations. Initially,
all servers are idle. The plot shows the share of servers of class $c=3$ that
have at least one job, that is $\mathbb{E}[X_{t,(3,1)}]$. We compare the
results for a model with $N=10$ servers per class (top right), $N=20$ (bottom
right), or $N=50$ (left). We observe that in all cases, the evolution of the
stochastic system is very well predicted by the mean field approximation. To
quantify this more precisely, each plot contains a zoom on the trajectory
between the time $t=8$ to $t=13$. These zooms show that for $N=50$, the
quantity $\mathbb{E}[X^{(N)}_{t,(3,1)}]$ is almost indistinguishable from the
mean field approximation. For $N=10$ or $N=20$, the estimated average is
slightly above the mean field curve, but the confidence intervals remain
almost equal to the error.

Figure 2. Illustration of the transient behavior of the CSMA model. We compare
the ‘average’ mean field and stochastic simulations for three different
scaling parameters: $N=10$, $N=20$ and $N=50$.

##### Steady-state and refined accuracy

While Theorems 1 and 2 provide a guarantee on the accuracy of the mean field
approximation, Theorem 3 shows that it is possible to compute an approximation
that is more accurate than the original mean field approximation. We
illustrate this in Figure 3 where we show the steady-state average queue
lengths for the same $5$ node graph. The sample mean and confidence interval
are computed from $40$ steady-state samples which again are obtained from
independent time-averages of $7.5\times 10^{6}$ events of the Markov chain
after a warm-up of $2.5\times 10^{6}$ events. For a class $c$ and a buffer
size $b$, the quantity $\sum_{b=1}^{B}\mathbb{E}[X_{(c,b)}]$ is equal to the
steady-state average queue length of each server of class $c$. In Figure 3, we
consider different values of $N$, and calculate the average queue length by
the following three methods:

  * •

By using a stochastic simulator of the original CSMA model.

  * •

By using the fixed point of the mean field approximation:
$\sum_{b=1}^{B}(\phi_{\text{\tiny$\infty$}})_{(c,b)}$.

  * •

By computing the refinement term, $C$, of Theorem 3 and using
$\sum_{b=1}^{B}(\phi_{\text{\tiny$\infty$}})_{(c,b)}+(C)_{(c,b)}/N$.

When looking at the scale of the $y$-axis, we see that in all cases, the
accuracy of the mean field approximation is already quite good. More
importantly, we also observe that in all cases, the refined approximation
seems almost exact: For all considered cases, the refined approximation lies
within the 95 percent confidence interval of the simulations and seems to work
well even for a small number of servers, $N\approx 10,20$. This result is
similar to the one observed for one-timescale mean field models in (Gast and
Van Houdt, 2017).

Figure 3. Queue Length Distribution for the 5 Node Graph of Figure (1(b)).

##### Computation Time

While the previous figure shows that the refined approximation provides an
increase in accuracy for small values of $N$, it comes at the cost of an
increase in computation time because one needs to compute the various
derivatives of the rate functions and to solve a new linear systems of
equations. In order to quantify the additional computation time, we measure
the time taken by our implementation to compute the refinement terms which are
reported in Table 1. We compare the $5$ node model studied before and a
$3$-node model whose interference graph is as in Figure 1(a). We observe that
the time taken to compute the refinement term is significant, in particular
for the $5$-node model. Yet, when looking more carefully at what takes time,
we realize that most of the computation time is taken by the symbolic
differentiation. Indeed, to simplify our implementation, we used the automated
differentiation method of sympy. While this yields simplifications for the
implementation, we encountered that it massively slows down the refinement
computation times. Through code profiling it showed that around 95 percent of
the computing time is taken by sympy methods such as differentiation and
evaluation of symbolic expressions. For smaller interference graphs, e.g.,
linear 2 / 3 node graphs, this effect is not limiting. For larger graphs, the
differentiation turns out to be the restricting factor. In Table 1 we state
the computation times for a linear 3 node model and for the setup described
before.

We would like to emphasize that the goal of our implementation is to
illustrate the theoretical statements, and thus we did not focus on
efficiency. The table shows that if one wants to adapt our implementation to
work with larger graphs, it would be sufficient to implement a more efficient
differentiation method. For instance, this could be done by using closed form
expression of the derivatives, or by using automatic differentiation methods,
or by using finite difference methods. We believe that such methods would
probably be much faster.

| Jakobian (A) | Hessian (B) | v + w | s + t + u | Total | Sympy
---|---|---|---|---|---|---
3 Node | 6.08 (6.08) | 30.71 (30.71) | 0.19 (0) | 5.64 (3.81) | 42.63 (40.6) | 95.22%
5 Node | 97.09 (97.09) | 897.15 (897.15) | 1\. (0) | 122.36 (86.09) | 1117.6 (1080.33) | 96.66%

Table 1. Refinement computation times for Random Access CSMA model for a 3
node linear graph and the 5 node graph as in Figures 1(a), 1(b). Times are
given in seconds. In parentheses, we indicate the time taken by
differentiation and subsequent sympy methods. These numbers show that the
Sympy code takes more than $95\%$ of the computation time.

##  6\. Conclusion

In this paper we investigate the accuracy of the classical _averaging method_
that is used to study two timescale models. We study a generic two timescale
model and show that under mild regularity conditions, the bias of this
‘average’ mean field approximation is of order $O(1/N)$. This result holds for
any finite time-horizon, and extends to the steady-state regime under the
classical assumption that the system has a unique and stable fixed point. Our
results show the existence of a bias term $C_{h}(t)$ for any regular function
$h$:

$\displaystyle\mathbb{E}[h(\boldsymbol{X}_{t})]=\underbrace{\underbrace{h(\Phi_{t}(\boldsymbol{x}))}_{\text{classical
`average' mean field}}+\underbrace{\frac{1}{N}C_{h}(t)}_{\text{$O(1/N)$
expansion of the bias}}}_{\text{refined `average' mean field}}+\quad o(1/N).$

For the steady-state regime $t=+\infty$, we propose an algorithmic method to
calculate this term $C_{h}$. This correction term can be computed by solving
linear systems and is therefore easily numerically computable. We show on an
example that, similarly to what was done for classical one timescale models
(Gast and Van Houdt, 2017), the bias term leads to an approximation that is
almost exact for small values of $N$ like $N=10,20$.

An interesting open question would be to obtain a characterization of
$C_{h}(t)$ for the transient regime. Yet, it is not clear to us if those
expressions would be usable as their size grows quickly with the system size.
From an application point of view, our examples show that the new
approximation leads to very accurate estimates for CSMA models. We believe
that the same should hold for other multiscale models.

###### Acknowledgements.

This work is supported by the French National Research Agency (ANR) through
REFINO Project under Grant ANR-19-CE23-0015.

## References

  * (1)
  * Ball et al. (2006) Karen Ball, Thomas G. Kurtz, Lea Popovic, and Greg Rempala. 2006\.  Asymptotic analysis of multiscale approximations to reaction networks.  _The Annals of Applied Probability_ 16, 4 (Nov. 2006), 1925–1961.  https://doi.org/10.1214/105051606000000420
  * Benaim and Boudec (2008) M. Benaim and J. . Y. Boudec. 2008.  A class of mean field interaction models for computer and communication systems.  _Performance Evaluation_ 65 (2008).  https://doi.org/10.1016/j.peva.2008.03.005
  * Boorstyn et al. (1987) R. Boorstyn, A. Kershenbaum, B. Maglaris, and V. Sahin. 1987.  Throughput Analysis in Multihop CSMA Packet Radio Networks.  _IEEE Transactions on Communications_ 35, 3 (1987), 267–274.  https://doi.org/10.1109/TCOM.1987.1096769
  * Bordenave et al. (2010) Charles Bordenave, David R. McDonald, and Alexandre Proutière. 2010.  A particle system in interaction with a rapidly varying environment: Mean field limits and applications.  _Networks and Heterogeneous Media_ 5, 1 (2010), 31–62.
  * Braverman (2022) Anton Braverman. 2022\.  The Prelimit Generator Comparison Approach of Stein’s Method.  _Stochastic Systems_ 12, 2 (2022), 181–204.  https://doi.org/10.1287/stsy.2021.0085
  * Braverman and Dai (2017) Anton Braverman and J. G. Dai. 2017.  Stein’s method for steady-state diffusion approximations of $M/\mathit{Ph}/n+M$ systems.  _The Annals of Applied Probability_ 27, 1 (feb 2017).  https://doi.org/10.1214/16-AAP1211
  * Braverman et al. (2016) Anton Braverman, J. G. Dai, and Jiekun Feng. 2016.  Stein’s method for steady-state diffusion approximations: An introduction through the Erlang-A and Erlang-C models.  _Stochastic Systems_ 6, 2 (2016), 301 – 366.  https://doi.org/10.1214/15-SSY212
  * Castiel et al. (2021) Eyal Castiel, Sem Borst, Laurent Miclo, Florian Simatos, and Phil Whiting. 2021.  Induced idleness leads to deterministic heavy traffic limits for queue-based random-access algorithms.  _The Annals of Applied Probability_ 31, 2 (2021), 941 – 971.  https://doi.org/10.1214/20-AAP1609
  * Cecchi (2018) F. Cecchi. 2018\.  _Mean-field limits for ultra-dense random-access networks_.  Technische Universiteit Eindhoven.  Proefschrift.
  * Cecchi et al. (2021) Fabio Cecchi, Sem C. Borst, Johan S. H. van Leeuwaarden, and Philip A. Whiting. 2021.  Mean-Field Limits for Large-Scale Random-Access Networks.  _Stochastic Systems_ 11, 3 (2021), 193–217.  https://doi.org/10.1287/stsy.2021.0068 arXiv:https://doi.org/10.1287/stsy.2021.0068
  * Feuillet and Robert (2014) Mathieu Feuillet and Philippe Robert. 2014.  A Scaling Analysis of a Transient Stochastic Network.  _Advances in Applied Probability_ 46, 2 (jun 2014), 516–535.  https://doi.org/10.1239/aap/1401369705
  * Gast (2017) Nicolas Gast. 2017\.  Expected Values Estimated via Mean-Field Approximation are 1/N-Accurate.  _Proceedings of the ACM on Measurement and Analysis of Computing Systems_ 1, 1 (jun 2017), 1–26.  https://doi.org/10.1145/3084454
  * Gast et al. (2019) Nicolas Gast, Luca Bortolussi, and Mirco Tribastone. 2019\.  Size expansions of mean field approximation: Transient and steady-state analysis.  _Performance Evaluation_ 129 (Feb. 2019), 60–80.  https://doi.org/10.1016/j.peva.2018.09.005
  * Gast and Gaujal (2012) Nicolas Gast and Bruno Gaujal. 2012.  Markov chains with discontinuous drifts have differential inclusion limits.  _Performance Evaluation_ 69, 12 (2012), 623–642.  https://doi.org/10.1177/075910639203700105
  * Gast and Van Houdt (2017) Nicolas Gast and Benny Van Houdt. 2017.  A Refined Mean Field Approximation.  _Proceedings of the ACM on Measurement and Analysis of Computing Systems_ 1, 2 (dec 2017), 33:1–33:28.  https://doi.org/10.1145/3154491
  * Hunt and Kurtz (1994) P.J. Hunt and T.G. Kurtz. 1994.  Large loss networks.  _Stochastic Processes and their Applications_ 53, 2 (oct 1994), 363–378.  https://doi.org/10.1016/0304-4149(94)90071-X
  * Hunter (1982) Jeffrey J. Hunter. 1982\.  Generalized inverses and their application to applied probability problems.  _Linear Algebra Appl._ 45 (1982), 157–198.  https://doi.org/10.1016/0024-3795(82)90218-X
  * Ipsen and Meyer (1994) Ilse C. F. Ipsen and Carl D. Meyer. 1994.  Uniform Stability of Markov Chains.  15, 4 (1994), 1061–1074.  https://doi.org/10.1137/S0895479892237562
  * Kang and Kurtz (2013) Hye-Won Kang and Thomas G. Kurtz. 2013.  Separation of time-scales and model reduction for stochastic reaction networks.  _The Annals of Applied Probability_ 23, 2 (apr 2013).  https://doi.org/10.1214/12-AAP841
  * Kang et al. (2014) Hye-Won Kang, Thomas G. Kurtz, and Lea Popovic. 2014\.  Central limit theorems and diffusion approximations for multiscale Markov chain models.  _The Annals of Applied Probability_ 24, 2 (apr 2014).  https://doi.org/10.1214/13-AAP934
  * Kurtz (1970) Thomas G. Kurtz. 1970\.  Solutions of ordinary differential equations as limits of pure jump markov processes.  _Journal of Applied Probability_ 7, 1 (April 1970), 49–58.  https://doi.org/10.2307/3212147
  * Kurtz (1978) Thomas G. Kurtz. 1978\.  Strong approximation theorems for density dependent Markov chains.  _Stochastic Processes and their Applications_ 6, 3 (feb 1978), 223–240.  https://doi.org/10.1016/0304-4149(78)90020-0
  * Le Boudec et al. (2007) Jean-Yves Le Boudec, David McDonald, and Jochen Mundinger. 2007\.  A Generic Mean Field Convergence Result for Systems of Interacting Objects. In _Fourth International Conference on the Quantitative Evaluation of Systems (QEST 2007)_. IEEE, Edinburgh, Scotland, UK, 3–18.  https://doi.org/10.1109/QEST.2007.8
  * Meurer et al. (2017) Aaron Meurer, Christopher P Smith, Mateusz Paprocki, Ondřej Čertík, Sergey B Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K Moore, Sartaj Singh, et al. 2017\.  SymPy: symbolic computing in Python.  _PeerJ Computer Science_ 3 (2017), e103.
  * Mitzenmacher (2001) M. Mitzenmacher. 2001\.  The power of two choices in randomized load balancing.  _IEEE Transactions on Parallel and Distributed Systems_ 12, 10 (oct 2001), 1094–1104.  https://doi.org/10.1109/71.963420
  * Pavliotis and Stuart (2008) Grigorios A. Pavliotis and Andrew M. Stuart. 2008.  _Multiscale Methods: Averaging and Homogenization_. Texts Applied in Mathematics, Vol. 53.  Springer New York, New York, NY.
  * Pazy (1983) A. Pazy. 1983.  _Semigroups of Linear Operators and Applications to Partial Differential Equations_. Applied Mathematical Sciences, Vol. 44.  Springer New York, New York, NY.  https://doi.org/10.1007/978-1-4612-5561-1
  * Perko (2001) Lawrence Perko. 2001\.  _Differential Equations and Dynamical Systems_ (3rd ed.).  Number 7 in Texts in Applied Mathematics. Springer-Verlag, Berlin, Heidelberg.
  * Robert and Vignoud (2021) Philippe Robert and Gaëtan Vignoud. 2021.  Stochastic Models of Neural Synaptic Plasticity: A Scaling Approach.  _SIAM J. Appl. Math._ 81, 6 (2021), 2362–2386.  https://doi.org/10.1137/20M1382891
  * Stein (1986) Charles Stein. 1986\.  _Approximate Computation of Expectations_. Vol. 7.  Institute of Mathematical Statistics. i–164 pages.  http://www.jstor.org/stable/4355512
  * van de Ven et al. (2010) P. M. van de Ven, S. C. Borst, J. S. H. van Leeuwaarden, and A. Proutière. 2010.  Insensitivity and stability of random-access networks.  _Performance Evaluation_ 67, 11 (2010), 1230–1242.  https://doi.org/10.1016/j.peva.2010.08.011
  * Wang and Kar (2005) Xin Wang and K. Kar. 2005\.  Throughput modelling and fairness issues in CSMA/CA based ad-hoc networks. In _Proceedings IEEE 24th Annual Joint Conference of the IEEE Computer and Communications Societies._ (2005-03), Vol. 1. 23–34 vol. 1\.  https://doi.org/10.1109/INFCOM.2005.1497875 ISSN: 0743-166X.
  * Ying (2016) Lei Ying. 2016.  On the Approximation Error of Mean-Field Models. In _Proceedings of the 2016 ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Science - SIGMETRICS ’16_. ACM Press, Antibes Juan-les-Pins, France, 285–297.  https://doi.org/10.1145/2896377.2901463

##  Appendix A Definitions

In this section we revise some essential definitions and properties used in
the paper. As these definitions are well established, we only briefly recall
them to provide a self-contained paper.

###  A.1. $C_{0}$-Semi-Group, Hölder Norm, ODE differentiability

###### Definition 1 ($C_{0}$ Semi-Group (Pazy, 1983) Definition 2.1 ).

$T_{s}$ is called strongly continuous semi-group (or $C_{0}$-semi-group) if

(22) $\displaystyle T_{0}=Id,$ $\displaystyle T_{s+t}=T_{s}T_{t}\text{ for all
}s,t\geq 0,$ $\displaystyle\lim_{t\downarrow 0}T_{t}z=z\text{ for all }z.$

Since we will only work with $C_{0}$-semi-groups we will simply refer to them
as semi-groups.
Semi-group & Generator commutation The generator of a $C_{0}$-semi-group is
defined by

$\displaystyle Az=\lim_{t\to 0}\frac{1}{t}(T_{t}z-T_{0}z).$

A direct consequence of the definition of the generator and a standard
property is that it commutes with its defining $C_{0}$-semi-group, i.e.,

$\displaystyle AT_{s}h(z)=T_{s}Ah(z)=\frac{d}{ds}T_{s}h(z).$

Note, this therefore holds true for the semi-groups given by the stochastic
system $\Psi_{s}h(x,y)=\mathbb{E}[h(X_{s},Y_{s})\mid X_{0},Y_{0}=x,y]$ with
generator $Lh(x,y)$ as in Equation (6) and the semi-group of the ODE
$\Phi_{s}h(x)=h(\phi_{s}(x))$ with generator $\Lambda
h(x)=D_{x}h(x)\bar{F}(x)$.

###### Definition 2 (Hölder Norm and Space).

For $U\subset\mathbb{R}^{n}$ and $u\in C^{k}(U)$

$\displaystyle\left\|u\right\|_{k,\gamma}:=\sum_{k}\sup_{x\in
U}\left\|D^{k}u(x)\right\|+\sum_{k}\sup\left\\{\frac{\left\|D^{k}u(x)-D^{k}u(y)\right\|}{\left\|x-y\right\|^{\gamma}}\mid
x,y\in U,x\neq y\right\\}$

is called Hölder norm. The space of functions for which the norm is finite is
called Hölder space and denoted by $\mathcal{D}^{k}_{\gamma}(U)$. For the case
$\gamma=1$ the Hölder space encloses all functions who are $k$-times
continuously differentiable with bounded derivatives and who’s $k$-th
derivatives are Lipschitz continuous. The latter we simply denote by
$\mathcal{D}^{k}(U)$.

An important implication is that all Hölder continuous functions are uniformly
continuous.

###### Lemma 3 (Drift induced differentiability (Perko, 2001) Theorem 1 p.80).

Let $E$ be an open subset of $\mathbb{R}^{n}$ containing $x_{0}$ and assume
that $f\in C^{k}(E)$. Then there exists an $a>0$ and $\delta>0$ such that for
all $y\in
N_{\delta}(x_{0})$333$N_{\delta}(x_{0}):=\\{x\in\mathbb{R}^{n}:\|x-x_{0}\|\leq\delta\\}$
the initial value problem

$\displaystyle\dot{x}=F(x),\qquad x(0)=y,$

has a unique solution $u$ which is $k$-times continuously differentiable with
respect to the initial condition for $t\in[-a,a]$.

###### Proof.

Theorem 1 p.80-83 and Remark 1 p.83 of (Perko, 2001). ∎

##  Appendix B Technical Lemmas and Proofs

###  B.1. Proof of Lemma 1

By assumption ($A_{\text{2}}$), the transition matrix $K(x)$ has a unique
irreducible class. As pointed out in Section 2.3, the corresponding Markov
chain has a unique stationary distribution that we denote by $\pi(x)$. Let
$\Pi(x):=\mathds{1}\pi^{T}(x)$ be the matrix whose lines are all equal to
$\pi(x)$. By (Hunter, 1982)(Theorem 3.5, p.17), $\bigl{(}K(x)+\Pi(x)\bigr{)}$
is non-singular and its inverse $(K(x)+\Pi(x))^{-1}$ is a generalized inverse
to $K(x)$, which means that it satisfies $K(x)(K(x)+\Pi(x))^{-1}K(x)=K(x)$.

To obtain the solution to the Poisson equation (10) we only consider the case
where $h$ takes values in $\mathbb{R}$. The extension to a function $h$ that
takes values in $\mathbb{R}^{n}$ is straightforward as it corresponds to $n$
independent Poisson equations.

Let us suppress the dependence on $x$ for clarity. Recall that
$K^{+}=(K+\Pi)^{-1}(I-\Pi)$ and let us study the product $KK^{+}$:

$\displaystyle KK^{+}$ $\displaystyle=K(K+\Pi)^{-1}(I-\Pi)$ (by definition)
$\displaystyle=K(K+\Pi)^{-1}(I+K-K-\Pi)$
$\displaystyle=K(K+\Pi)^{-1}+K(K+\Pi)^{-1}K-K(K+\Pi)^{-1}(K+\Pi)$ (expanding
the product) $\displaystyle=K(K+\Pi)^{-1}$ (the last two terms equal $\pm K$)
$\displaystyle=(K+\Pi)(K+\Pi)^{-1}-\Pi(K+\Pi)^{-1}$ (Adding and subtracting
$\Pi(K+\Pi)^{-1}$) $\displaystyle=I-\Pi,$

where the last equality holds because $\pi^{T}\mathds{1}=1$ and therefore
$\Pi\ \Pi=\mathds{1}\pi^{T}\mathds{1}\pi^{T}=\mathds{1}\pi^{T}=\Pi$. Combined
with $\Pi\ K=0$, this shows that $\Pi=\Pi(K+\Pi)$ and therefore
$\Pi(K+\Pi)^{-1}=\Pi$.

The above computations show that if
$G^{\text{fast}}_{h}(x,y)=\sum_{y^{\prime}}K^{+}_{y,y^{\prime}}(x)h(x,y^{\prime})$,
then:

$\displaystyle K(x)G^{\text{fast}}_{h}(x,y)$
$\displaystyle=h(x,y)-\sum_{y^{\prime}}\pi_{y^{\prime}}(x)h(x,y^{\prime}),$

which shows that $G^{\text{fast}}_{h}$ is the solution of the Poisson
equation.

The differentiability of $G^{\text{fast}}_{h}$ follows from the
differentiability of $h$ and $K^{+}$: Under Assumption ($A_{\text{1}}$) $K(x)$
is continuously differentiable. By Assumption ($A_{\text{2}}$), $K(x)$ has a
unique irreducible class, this implies $\Pi(x)$ is continuously
differentiable, for which we refer to (Ipsen and Meyer, 1994), and therefore
further implies that $(I+\Pi(x))^{-1}$ and $K^{+}(x)$ are continuously
differentiable. This proves that $G^{\text{fast}}_{h}$ is continuously
differentiable in $x$.

###  B.2. Technical Lemmas used to prove Theorem 1 (Transient Regime)

###### Lemma 1.

For arbitrary but fixed $t>0$, let
$g:(s,x,y)\in[0,t]\times\mathcal{X}\times\mathcal{Y}\mapsto
g_{s}(x,y)\in\mathbb{R}$ be a continuous function that is continuously
differentiable in $s$ and let $H_{s}(x,y):=\mathbb{E}[g_{s}(X_{s},Y_{s})\mid
X_{0},Y_{0}=x,y]$. Then:

$\displaystyle
H_{t}(x,y)-H_{0}(x,y)=\int_{0}^{t}\mathbb{E}[\frac{d}{d\tau}g_{s+\tau}(X_{s},Y_{s})\mid_{\tau=0}]ds+\int_{0}^{t}\mathbb{E}[Lg_{s}(X_{s},Y_{s})]ds.$

###### Proof.

By definition of the generator $L$, the quantity
$H_{t}(x,y)=\mathbb{E}[g_{t}(X_{t},Y_{t})\mid X_{0},Y_{0}=x,y]$ is right sided
differentiable, i.e., $\frac{d^{+}}{ds}f(s)=\lim_{ds\to
0^{+}}\frac{1}{ds}\bigl{(}f(s+ds)-f(s)\bigr{)}$, with respect to time. Using
semi-group properties and bounds, its derivative is

$\displaystyle\frac{d^{+}}{ds}H_{s}(x,y)$
$\displaystyle=\mathbb{E}[\frac{d}{d\tau}g_{s+\tau}(X_{s},Y_{s})\mid_{\tau=0}]+\mathbb{E}[Lg_{s}(X_{s},Y_{s})].$

The first term corresponds to the derivation of $g_{s}$ with respect to time
and the second term to the changes of the stochastic system in
$(X_{s},Y_{s})$. The lemma therefore follows by using that
$H_{t}(x,y)-H_{0}(x,y)=\int_{0}^{t}\frac{d^{+}}{ds}H_{s}(X_{s},Y_{s})$. ∎

###### Lemma 2 (Bound for $L_{\text{fast}}g(x,y)$).

Take the two-timescale stochastic system as introduced in Section 2.1 with
generator $L$ and assume ($A_{\text{1}}$) and ($A_{\text{2}}$). Let
$L_{\text{fast}}$ be as defined in Section 4.2. Then, for
$g\in\mathcal{D}^{1}(\mathcal{X}\times\mathcal{Y})$ and
$(x,y)\in\mathcal{X}\times\mathcal{Y}$

$\displaystyle
L_{\text{fast}}g(x,y)-\frac{1}{N}Lg(x,y)=\frac{1}{N}\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(x,y)\ell\
D_{x}g(x,y^{\prime})+o(1/N).$

###### Proof.

By definition of $L_{\text{fast}}$, for a continuous function $g$ the values
of $L_{\text{fast}}g(x,y)$ and $\frac{1}{N}Lg(x,y)$ coincide in the limit. For
finite $N$, we first look at $\frac{1}{N}Lg(x,y)$ which is given by

$\displaystyle\frac{1}{N}Lg(x,y)=\frac{1}{N}\sum_{\ell,y^{\prime}}N\alpha_{\ell,y^{\prime}}(x,y)(g(x+\frac{\ell}{N},y^{\prime})-g(x,y)).$

Using the continuity of $g$ in $x$, and definition of $L_{\text{fast}}g(x,y)$
we have

$\displaystyle L_{\text{fast}}g(x,y)=\lim_{N\to\infty}\frac{1}{N}Lg(x,y)$
$\displaystyle=\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(x,y)(\lim_{N\to\infty}g(x+\frac{\ell}{N},y^{\prime})-g(x,y))$
$\displaystyle=\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(x,y)(g(x,y^{\prime})-g(x,y)).$

Using Taylor’s theorem

$\displaystyle Lg(x,y)$
$\displaystyle=N\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(x,y)(g(x+\frac{\ell}{N},y^{\prime})-g(x,y))$
$\displaystyle=N\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(x,y)(g(x,y^{\prime})-g(x,y)+\ell\
D_{x}g(x,y^{\prime})+o(\left\|\ell\right\|))$
$\displaystyle=NL_{\text{fast}}g(x,y)+\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(x,y)\ell\
D_{x}g(x,y^{\prime})+o(1).$

∎

###  B.3. Technical Lemma used to prove Theorem 3 (steady-state refinement)

The next Lemma justifies that the first term of (20) is approximated in the
limit by the ‘average’ version of $J_{h}$ as defined in Section 4.5. Second,
Lemma 4 gives stability conditions for the solution of the ‘slow’ Poisson
Equation.

###### Lemma 3.

Assume ($A_{\text{1}}$)-($A_{\text{3}}$) in particular assume
$h\in\mathcal{D}^{1}(\mathcal{X})$. Further assume that $G^{\text{fast}}_{F}$,
solution to the ‘fast’ Poisson Equation, is cont. differentiable in $x$ and
$G^{\text{slow}}_{h}$, solution to the ‘slow’ Poisson Equation, is twice cont.
differentiable in $x$. Let
$\bar{J}(\phi_{\text{\tiny$\infty$}})=\sum_{y}\pi(\phi_{\text{\tiny$\infty$}})J(\phi_{\text{\tiny$\infty$}},y)$
be the ‘average’ version of
$J_{h}(x,y)=\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(x,y)\
D_{x}\bigl{(}D_{x}G^{\text{slow}}_{h}(x)G^{\text{fast}}_{F}(x,y^{\prime})\bigr{)}\ell$
in $\phi_{\text{\tiny$\infty$}}$, then

$\displaystyle
N\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})\bigl{(}\bar{F}(\boldsymbol{X}_{\infty})-F(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})\bigr{)}]=\bar{J}(\phi_{\text{\tiny$\infty$}})+o(1).$

###### Proof.

To prove the lemma, let $G^{\text{fast}}_{F}$ be the solution to the Poisson
equation as given in (‘fast’ Poisson Equation). We use $G^{\text{fast}}_{F}$
to rewrite

$\displaystyle N$
$\displaystyle\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})(\bar{F}(\boldsymbol{X}_{\infty})-F(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty}))]=-N\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})(F(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})-\bar{F}(\boldsymbol{X}_{\infty}))]$
$\displaystyle=-N\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})(L_{\text{fast}}G^{\text{fast}}_{F}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty}))]=-N\mathbb{E}[L_{\text{fast}}D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})(G^{\text{fast}}_{F}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty}))].$

Adding
$\mathbb{E}[LD_{x}G^{\text{slow}}(\boldsymbol{X}_{\infty})(G^{\text{fast}}_{h}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty}))]=0$
and applying the steps as in the proof of Lemma 2 we see that

$\displaystyle-N\mathbb{E}[L_{\text{fast}}D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})(G^{\text{fast}}_{F}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty}))-LD_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})(G^{\text{fast}}_{F}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty}))]$
(23)
$\displaystyle=-\mathbb{E}[\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})\bigl{(}D_{x}(D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})G^{\text{fast}}_{F}(\boldsymbol{X}_{\infty},y^{\prime}))\bigr{)}\ell\
]+o(1).$

The last equality follows directly from the definition of $L_{\text{fast}}$
since $G_{h}$ only depends on $x$. Using Theorem 2 we see that
$\mathbb{E}[J_{h}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})]$ with
$J_{h}(x,y):=-\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(x,y)\ell\
D_{x}\bigl{(}D_{x}G^{\text{slow}}_{h}(x)G^{\text{fast}}(x,y^{\prime})\bigr{)}$
of equation (23) is approximated by
$\bar{J}(\phi_{\text{\tiny$\infty$}})=\sum_{y}\pi(\phi_{\text{\tiny$\infty$}})J(\phi_{\text{\tiny$\infty$}},y)$.
This concludes the proof as it implies

$N\mathbb{E}[D_{x}G^{\text{slow}}_{h}(\boldsymbol{X}_{\infty})\bigl{(}\bar{F}(\boldsymbol{X}_{\infty})-F(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})\bigr{)}]=\mathbb{E}[J_{h}(\boldsymbol{X}_{\infty},\boldsymbol{Y}_{\infty})]+o(1)=\bar{J}_{h}(\phi_{\text{\tiny$\infty$}})+o(1).$
∎

###  B.4. Stability of $G^{\text{slow}}_{h}$

###### Lemma 4 (Stability).

Assume that $\bar{F}$ and $\phi$ are k-times differentiable with uniformly
continuous derivatives and that $\phi$ has a unique exponentially stable
attractor $\phi_{\text{\tiny$\infty$}}$. Then the k-th derivative of
$G^{\text{slow}}_{h}:x\mapsto\int_{0}^{\infty}h(\phi_{s}x)-h(\phi_{\text{\tiny$\infty$}})ds$
is bounded and equal to $\int_{0}^{\infty}D^{k}_{x}(h\circ\phi_{s})(x)ds$.

###### Proof.

This is a consequence of (Gast and Van Houdt, 2017)[Lemma 3.5]. ∎

###  B.5. Proof of Proposition 4

###### Proof.

In the first part of the proof we find computable expressions of
$D_{x}G^{\text{slow}}_{h},D_{x}^{2}G^{\text{slow}}_{h}$ and
$D_{x}G^{\text{fast}}_{F}$. These expressions allow us to rewrite $J_{h}$ and
construct the closed form representation of $J_{h}$ as well as $\bar{J}_{h}$
in the following steps. By definition

$\displaystyle J_{h}(x,y)$
$\displaystyle=-\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(x,y)\
D_{x}\bigl{(}D_{x}G^{\text{slow}}_{h}(x)G^{\text{fast}}_{F}(x,y^{\prime})\bigr{)}\ell$
(24) $\displaystyle=-\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(x,y)\
\bigl{(}G^{\text{fast}}_{F}(x,y^{\prime})^{T}D^{2}_{x}G^{\text{slow}}_{h}(x)+D_{x}G^{\text{slow}}_{h}(x)D_{x}G^{\text{fast}}_{F}(x,y^{\prime})\bigr{)}\ell.$

By Lemma 1 there exists a matrix $K^{+}(x)$ such that
$G^{\text{fast}}_{F}(x,y^{\prime})$ has the form

$\displaystyle
G^{\text{fast}}_{F}(x,y^{\prime})=\sum_{y^{\prime}\in\mathcal{Y}}K^{+}_{y^{\prime},y^{\prime\prime}}(x)F(x,y^{\prime\prime}).$

Assumption ($A_{\text{1}}$) which assures differentiability of the transition
rates $\alpha$ with respect to $x$, also implies differentiability for
$G^{\text{fast}}_{F}$. Therefore,
$D_{x}G^{\text{fast}}_{F}(x,y^{\prime})=\sum_{y^{\prime\prime}\in\mathcal{Y}}F(x,y^{\prime\prime})\nabla^{T}_{x}K^{+}_{y^{\prime},y^{\prime\prime}}(x)+K^{+}_{y^{\prime},y^{\prime\prime}}(x)D_{x}F(x,y^{\prime\prime})$
with $\nabla_{x}^{T}=\left[\frac{\partial}{\partial
x_{1}},\dots,\frac{\partial}{\partial
x_{n}}\right],D_{x}f=(\frac{\partial}{\partial x_{j}}f_{i})_{i,j=1\dots n}$.
Using the results of (Gast and Van Houdt, 2017)[Lemma 3.6] with $A,B$ the
first and second derivative of $\bar{F}$ as defined in (21), it holds that for
the equilibrium point $\phi_{\text{\tiny$\infty$}}$

$\displaystyle D_{x}G^{\text{slow}}_{h}(\phi_{\text{\tiny$\infty$}})_{i}$
$\displaystyle=\sum_{j}\frac{\partial h}{\partial
x_{j}}(\phi_{\text{\tiny$\infty$}})\int_{0}^{\infty}\bigl{(}D_{x}\phi_{s}(\phi_{\text{\tiny$\infty$}})\bigr{)}_{j,i}ds=\sum_{j}\frac{\partial
h}{\partial x_{j}}(\phi_{\text{\tiny$\infty$}})(-A)^{-1}_{j,i}$

as well as

(25) $\displaystyle
D^{2}_{x}G^{\text{slow}}_{h}(\phi_{\text{\tiny$\infty$}})_{n,m}$
$\displaystyle=\sum_{i,j}\frac{\partial^{2}h}{\partial x_{i}\partial
x_{j}}(\phi_{\text{\tiny$\infty$}})\int_{0}^{\infty}\bigl{(}D_{x}\phi_{s}(\phi_{\text{\tiny$\infty$}})\bigr{)}_{j,n}\bigl{(}D_{x}\phi_{s}(\phi_{\text{\tiny$\infty$}})\bigr{)}_{i,m}ds$
$\displaystyle\qquad+\sum_{i}\frac{\partial h}{\partial
x_{i}}(\phi_{\text{\tiny$\infty$}})\int_{0}^{\infty}\bigl{(}D_{x}^{2}\phi_{s}(\phi_{\text{\tiny$\infty$}})\bigr{)}_{j,m,n}ds$
$\displaystyle=\sum_{i,j}\frac{\partial^{2}h}{\partial x_{i}\partial
x_{j}}(\phi_{\text{\tiny$\infty$}})\int_{0}^{\infty}\bigl{(}e^{As}\bigr{)}_{j,n}\bigl{(}e^{As}\bigr{)}_{i,m}ds$
$\displaystyle\qquad+\sum_{i}\frac{\partial h}{\partial
x_{i}}(\phi_{\text{\tiny$\infty$}})\bigl{(}\sum_{j}(-A)^{-1}_{i,j}\sum_{k_{1},k_{2}}B_{j,k_{1},k_{2}}\int_{0}^{\infty}(e^{As}\bigr{)}_{k_{1},n}(e^{As}\bigr{)}_{k_{2},m}ds\bigr{)}.$

Next, the above equations are used to rewrite (24). To obtain a closed form
expression for the left summand, we start by looking at the sum

(26)
$\displaystyle\sum_{y}\pi_{y}(\phi_{\text{\tiny$\infty$}})\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\phi_{\text{\tiny$\infty$}},y)\sum_{y^{\prime}}K^{+}_{y^{\prime},y^{\prime}}F_{m}(\phi_{\text{\tiny$\infty$}},y^{\prime})\ell_{n}\int_{0}^{\infty}(e^{As}\bigr{)}_{k_{1},n}(e^{As}\bigr{)}_{k_{2},m}.$

A solution to the above equation is given by the following Lyapunov equation.
To ease notations, define

$\displaystyle O_{m,n}$
$\displaystyle:=\sum_{y}\pi_{y}(\phi_{\text{\tiny$\infty$}})\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\phi_{\text{\tiny$\infty$}},y)K^{+}_{y^{\prime},:}(\phi_{\text{\tiny$\infty$}})F_{m}(\phi_{\text{\tiny$\infty$}},:)\ell_{n},$

or, equivalently in matrix notation,
$O\phantom{:}=\sum_{y}\pi_{y}(\phi_{\text{\tiny$\infty$}})\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\phi_{\text{\tiny$\infty$}},y)K^{+}_{y^{\prime},:}(\phi_{\text{\tiny$\infty$}})F(\phi_{\text{\tiny$\infty$}},:)\ell^{T}.$

If a matrix $U$ solves444As $A$ is non-singular, $A$ and $-A^{T}$ don’t share
any eigenvalues and therefore equation (27) has a unique solution. the
Sylvester equation (for $X$)

(27) $\displaystyle AX+XA^{T}=-O,$

it is equal to (26). Applying this identity and (25) to the first summand of
$\bar{J}(\phi_{\text{\tiny$\infty$}})$ which is the ‘average’ version of (24),
lets us rewrite

$\displaystyle\sum_{y}\pi_{y}(\phi_{\text{\tiny$\infty$}})\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\phi_{\text{\tiny$\infty$}},y)\
G^{\text{fast}}_{F}(\phi_{\text{\tiny$\infty$}},y^{\prime})^{T}D^{2}_{x}G^{\text{slow}}_{h}(\phi_{\text{\tiny$\infty$}})\ell$
$\displaystyle\qquad=\sum_{i,j}\frac{\partial^{2}h}{\partial x_{i}\partial
x_{j}}(\phi_{\text{\tiny$\infty$}})U_{i,j}+\sum_{i}\frac{\partial h}{\partial
x_{i}}(\phi_{\text{\tiny$\infty$}})\sum_{j}(-A)^{-1}_{i,j}\sum_{k_{1},k_{2}}B_{j,k_{1},k_{2}}U_{k_{1},k_{2}}.$

For the second major summand appearing in the definition of
$\bar{J}(\phi_{\text{\tiny$\infty$}})$, writing out the solutions to the
Poisson equations and their derivatives yields

$\displaystyle\sum_{y}\pi_{y}(\phi_{\text{\tiny$\infty$}})\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\phi_{\text{\tiny$\infty$}},y)\
D_{x}G^{\text{slow}}_{h}(\phi_{\text{\tiny$\infty$}})D_{x}G^{\text{fast}}_{F}(\phi_{\text{\tiny$\infty$}},y^{\prime})\ell$
$\displaystyle=\sum_{y}\pi_{y}(\phi_{\text{\tiny$\infty$}})\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\phi_{\text{\tiny$\infty$}},y)\times$
$\displaystyle\qquad\sum_{j,k}\left(\sum_{i}\frac{\partial h}{\partial
x_{i}}(\phi_{\text{\tiny$\infty$}})(-A)^{-1}_{i,k}\right)\left(\sum_{y^{\prime}}\frac{\partial}{\partial
x_{j}}K_{y^{\prime},y^{\prime}}^{+}(\phi_{\text{\tiny$\infty$}})F_{k}(\phi_{\text{\tiny$\infty$}},y^{\prime})+K_{y^{\prime},y^{\prime}}^{+}(\phi_{\text{\tiny$\infty$}})\frac{\partial}{\partial
x_{j}}F_{k}(\phi_{\text{\tiny$\infty$}},y^{\prime})\right)\ell_{j}.$

By using vector notation and rearranging the sums this is equal to

$\displaystyle\sum_{i}\frac{\partial h}{\partial
x_{i}}(\phi_{\text{\tiny$\infty$}})\sum_{k}(-A)^{-1}_{i,k}\sum_{y}\pi_{y}(\phi_{\text{\tiny$\infty$}})\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\phi_{\text{\tiny$\infty$}},y)\times$
(28)
$\displaystyle\qquad\left(\sum_{y^{\prime}}F_{k}(\phi_{\text{\tiny$\infty$}},y^{\prime})\nabla^{T}_{x}K_{y^{\prime},y^{\prime}}^{+}(\phi_{\text{\tiny$\infty$}})\ell+K_{y^{\prime},y^{\prime}}^{+}(\phi_{\text{\tiny$\infty$}})\nabla^{T}_{x}F_{k}(\phi_{\text{\tiny$\infty$}},y^{\prime})\ell\right).$

Lastly, define
$T_{i}:=\sum_{j}A^{-1}_{i,j}\sum_{k_{1},k_{2}}B_{j,k_{1},k_{2}}U_{k_{1},k_{2}}$
and

$\displaystyle
S_{i}:=\sum_{k}A^{-1}_{i,k}\sum_{y}\pi_{y}(\phi_{\text{\tiny$\infty$}})\sum_{\ell,y^{\prime}}\alpha_{\ell,y^{\prime}}(\phi_{\text{\tiny$\infty$}},y)\left(\sum_{y^{\prime}}F_{k}(\phi_{\text{\tiny$\infty$}},y^{\prime})\nabla^{T}_{x}K_{y^{\prime},y^{\prime}}^{+}(\phi_{\text{\tiny$\infty$}})\ell+K_{y^{\prime},y^{\prime}}^{+}(\phi_{\text{\tiny$\infty$}})\nabla^{T}_{x}F_{k}(\phi_{\text{\tiny$\infty$}},y^{\prime})\ell\right).$

This concludes the proof as by definition of $T$,$S$ and $U$,
$\bar{J}(\phi_{\text{\tiny$\infty$}})=\sum_{i}\frac{\partial h}{\partial
x_{i}}(\phi_{\text{\tiny$\infty$}})(T_{i}+S_{i})+\sum_{i,j}\frac{\partial^{2}h}{\partial
x_{i}\partial x_{j}}(\phi_{\text{\tiny$\infty$}})U_{i,j}$. ∎

##  Appendix C Computational Notes

The first note describes how to calculate the steady-state probabilities
$\pi(x)$ associated to the transition matrix $K(x)$.

###### Note 1 (Note on the computation for the stationary probabilities.).

As the finite state continuous time Markov chain with generator $K(x)$ has
unique irreducible class, there exists a non-trivial unique stationary
distribution $\pi(x)$555As the Markov chain is allowed to have transient
states the stationary distribution can take zero values.. Denote by
$\\{1,\dots,m\\}$ the states of the Markov chain. $\pi(x)$ is obtained by
solving the linear system

$\displaystyle v\ K(x)=\mathbf{0},\quad\sum_{i=1}^{m}v_{i}=1,$

with $v\geq 0$ component wise. By definition, $K(x)$ is of rank $m-1$. Using
its structure, i.e., $\sum_{y^{\prime}}K(x)_{y,y^{\prime}}=0\ \forall y$, we
can rewrite the above over-determined linear system by replacing the last
column of the generator with $\mathds{1}=[1,\dots,1]^{T}$ yielding

$\displaystyle v\
[K(x)_{:,1},\dots,K(x)_{:,m-1},\mathds{1}]=[\underbrace{0,\dots,0}_{m-1\text{
times}},1],$

for which the solution is the stationary distribution $\pi(x)$. By
$K(x)_{:,y},\ y=1,\dots,m-1$ we denote the $y$-th column of $K(x)$.

As shows in Lemma 1 the solution to the ’fast’ Poisson Equation 10 has the
form
$G^{\text{fast}}_{h}(x,y)=\sum_{y^{\prime}}K^{+}_{y,y^{\prime}}(x)h(x,y^{\prime})$.
To compute the bias correction terms, it is necessary to calculate the first
derivative of $G^{\text{fast}}_{h}$ with respect to $x$. Since the computation
of the derivative of $(K(x)+\Pi(x))^{-1}$ can be non-trivial and time
consuming, the note below elaborates how the derivative can be efficiently
obtained.

###### Note 2 (Computation of $D_{x}K^{+}(x)$).

By definition $K^{+}(x)=(K(x)+\Pi(x))^{-1}(I-\Pi(x))$. Using basic matrix
derivation rules one has

$\displaystyle\frac{\partial}{\partial x_{i}}K^{+}(x)$
$\displaystyle=\frac{\partial}{\partial
x_{i}}\bigl{(}(K(x)+\Pi(x))^{-1}(I-\Pi(x))\bigr{)}$ (29)
$\displaystyle=\frac{\partial}{\partial
x_{i}}\bigl{(}K(x)+\Pi(x)\bigr{)}^{-1}\bigl{(}I-\Pi(x)\bigr{)}-\bigl{(}K(x)+\Pi(x)\bigr{)}^{-1}\frac{\partial}{\partial
x_{i}}\Pi(x).$

Since the numerical difficulty lies only the computation of
$D_{x}\bigl{(}K(x)+\Pi(x)\bigr{)}^{-1}$ we focus solely on it. Define
$E(x)=(K(x)+\Pi(x))$ and let $I$ be the identity matrix. As pointed out in
Lemma 1 $(K(x)+\Pi(x))$ is indeed invertible and thus for the partials
derivative $\frac{\partial}{\partial x_{i}},i=0\dots d_{x}$ the following
holds:

$\displaystyle\frac{\partial}{\partial x_{i}}I$
$\displaystyle=\frac{\partial}{\partial x_{i}}(E(x)E^{-1}(x))$
$\displaystyle\Leftrightarrow$ $\displaystyle 0$
$\displaystyle=(\frac{\partial}{\partial
x_{i}}E(x))E^{-1}(x)+E(x)\frac{\partial}{\partial x_{i}}E^{-1}(x)$
$\displaystyle\Leftrightarrow$ $\displaystyle\frac{\partial}{\partial
x_{i}}E^{-1}(x)$ $\displaystyle=-E^{-1}(x)(\frac{\partial}{\partial
x_{i}}E(x))E^{-1}(x)$ $\displaystyle\Leftrightarrow$
$\displaystyle\frac{\partial}{\partial x_{i}}(K(x)+\Pi(x))^{-1}$
$\displaystyle=-(K(x)+\Pi(x))^{-1}\bigl{(}\frac{\partial}{\partial
x_{i}}(K(x)+\Pi(x))\bigr{)}(K(x)+\Pi(x))^{-1}.$

The above can now be used to compute Equation (29).

##  Appendix D Numerical Results for the 3 Node Model

For completeness, we give the numerical results of the linear 3 node model of
Graph 1(a). To obtain the results we modify the parameters of the code of Code
Cell LABEL:lst:python_code to match with the 3 node setup. The new model is
defined as in Code Cell LABEL:lst:python_code_3_node.

⬇

# Graph structure (this is the 3 node example)

G = np.array([[0,1,0],

[1,0,1],

[0,1,0]])

# rates & buffer size

_lambda = np.array([0.4,0.2,0.5])

nu = np.array([1.2,2.,1.5])

mu = np.array([1.4,1.3,1.7])

buffer_size = 10

As seen in Figure 4, we obtain similar results as for the 5 node model when
considering steady-state average queue length values. The sample mean and
confidence interval are computed from $40$ steady-state samples in the same
manner as for the 5 node example. The refined approximation almost exactly
indicates the stationary value of the stochastic process even for small $N$
while the mean field approximation gets more accurate as $N$ grows.

Figure 4. Stationary Queue Length Distribution for the 3 Node Graph of Figure
(1(a)).