S. Allmeier and N. Gast Stochastic Systems

Accuracy of the Graphon Mean Field Approximation

Accuracy of the Graphon Mean Field Approximation for Interacting Particle
Systems

Sebastian Allmeier

Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, 38000 Grenoble, France
<EMAIL_ADDRESS>
Nicolas Gast

Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, 38000 Grenoble, France
<EMAIL_ADDRESS>
We consider a system of $N$ particles whose interactions are characterized by
a (weighted) graph $G^{N}$. Each particle is a node of the graph with an
internal state. The state changes according to Markovian dynamics that depend
on the states and connection to other particles. We study the limiting
properties, focusing on the dense graph regime, where the number of neighbors
of a given node grows with $N$. We show that when $G^{N}$ converges to a
graphon $G$, the behavior of the system converges to a deterministic limit,
the graphon mean field approximation. We obtain convergence rates depending on
the system size $N$ and cut-norm distance between $G^{N}$ and $G$. We apply
the results for two subcases: When $G^{N}$ is a discretization of the graph
$G$ with individually weighted edges; when $G^{N}$ is a random graph obtained
through edge sampling from the graphon $G$. In the case of weighted
interactions, we obtain a bound of order $\mathcal{O}(1/N)$. In the random
graph case, the error is of order $\mathcal{O}(\sqrt{\log(N)/N})$ with high
probability. We illustrate the applicability of our results and the numerical
efficiency of the approximation through two examples: a graph-based load-
balancing model and a heterogeneous bike-sharing system.

graphon, mean field approximation, interacting particle systems, random graph,
bias bound

##  1 Introduction

Mean field approximations are ubiquitously used in the study of large scale
stochastic systems. Its mathematical foundations date back to the
${}^{\prime}60$s and ${}^{\prime}70$s with the seminal papers of Kurtz (1970,
1971), McKean (1967) and Norman (1972). While originating from statistical
physics, the mean field methodology has found application in many areas, such
as communication networks Le Boudec et al. (2007), load balancing Mitzenmacher
(2001), and the study of epidemics Roy et al. (2023). The fundamental idea of
the mean field method is to represent the particle process by a Markovian
state descriptor based on averaged quantities of the system. A quantity
commonly used in this context is state occupancy processes, e.g., the averaged
load of stations for a bike sharing system Fricker and Gast (2016) or the
fraction of servers having at least a certain queue length for load balancing
systems Mitzenmacher (2001). To apply the ‘classical’ mean field method, it is
crucial that the particles of the system are homogeneous and therefore
interchangeable, ensuring the Markov property for the mentioned aggregate
quantities. Its applicability to real world scenarios is often hindered by the
intrinsically heterogeneous behavior of realistic models as heterogeneity
breaks the quintessential exchangeability assumption. Therefore, it is
necessary to keep track of the evolution of the whole set of particles in the
system which makes the application of ‘classical’ mean field methods
complicated, if not infeasible.

To overcome these difficulties, mean field models which use graphons as a
description of diversity of individuals and the variation of their interaction
have been proposed Budhiraja et al. (2019), Bet et al. (2024), Vizuete et al.
(2020), Keliger (2023), Bayraktar et al. (2023). One example for where this
approach finds widespread application are epidemic models for which graphons
characterize the epidemic propagation Delmas et al. (2023), Decreusefond et
al. (2012), Aurell et al. (2022), Roy et al. (2023). For many settings,
results show that deterministic graphon-based models are limiting processes
for sequences of graph-based interacting and densely connected particle
systems. While the use of the graphon mean field model is well justified, the
bias and convergence speed in its application can vary immensely depending on
the convergence properties of the considered graph sequence. Surprisingly, the
literature concerning the issue of convergence rates and accuracy is
relatively sparse. This paper contributes to filling this gap as we derive
generic bounds on the bias for the graphon approximation for finite-sized and
densely connected interacting particle systems.

###  1.1 Contributions

We provide bias bounds for the graphon mean field approximation for finite-
sized systems consisting of $N\in\mathbb{N}$ interacting particles for which
the graph $G^{N}$ models the connection of the population. Our results show
that it is possible to derive bias bounds which largely depend on the
convergence properties of the graph sequence $G^{N}$ and of its limiting
graphon $G$. To be more precise, we start from a stochastic interacting
particle model of finite size $N$, where each individual $k$ is characterized
by a time-varying state $S_{k}(t)$. The connection of the particle to the
population is given by the edges $(k,l)$ of a graph $G^{N}$. Based on this
description, we construct a (deterministic) integro-differential equation
based on the graphon $G$, and show that it has a unique solution $\BFx^{G}(t)$
that we call the graphon mean field approximation. This differential equation
is constructed such that for fixed $N$, $x^{G}_{k/N,s}(t)$ approximates the
probability of particle $k\in[N]$ to be in state $s$ at time $t$. Denoting
$\mathbb{P}(S^{G^{N}}_{k}(t)=s)$ this probability, our main result shows that

$\displaystyle\mathbb{P}(S^{G^{N}}_{k}(t)=s)=x^{G}_{k/N,s}(t)+\mathcal{O}(\frac{1}{N}+\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}),$

where
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}$
is the $L_{2}$ distance of the step graphon representation of $G^{N}$ and the
graphon $G$, which is equivalent to the distance implied by graphon typical
cut norm, see Theorem 4.2 for a precise description.

To show the extent of the result, we consider two specific graph sampling
strategies for $G^{N}$. In the first case (which we denote _deterministic
sampling_), $G^{N}$ is the discretization of the limiting graphon $G$, with
$G^{N}_{kl}:=G_{k/N,l/N}\in(0,1]$ being the strength of interaction between
two particles. We consider this case as it illustrates how our result can be
used to model In the second case (that we call _stochastic sampling_), $G^{N}$
is a random graph generated from $G$, where an edge is present between two
nodes $k$ and $\ell$ with probability $G(k/N,l/N)$. The second case
corresponds to a generalization of Erdős–Rényi graph and stochastic block
models to possibly non-uniform probabilities.

Imposing some mild assumptions such as a piecewise Lipschitz condition on the
graphon $G$ or symmetry of $G$, in the case of stochastic sampling, we can
bound the distance between the graph $G^{N}$ and the graphon $G$ by:

$\displaystyle\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}=\begin{cases}\mathcal{O}(N^{-1}).&\text{
(Deterministic Sampling)}\\\
\mathcal{O}(\sqrt{\frac{\log(N)}{N}})~{}\text{w.h.p.,}&\text{ (Stochastic
Sampling)}\end{cases}$

Here, with high probability (w.h.p) means that the right-hand side in case of
the stochastic sampling holds with probability at least $1-2/N$. The precise
assumptions and results are given in the Corollaries 4.4 and 4.3. We point out
that the main result applies to the considered cases but is not limited to
those. It is possible to consider other graph sampling strategies that satisfy
our assumptions and for which the $L_{2}$ distance between the step graphon
representation of $G^{N}$ and the graphon $G$ can be bounded. To illustrate
our results and emphasize their applicability, we provide two examples, one
each for the stochastic and deterministic sampling method. Our first example
considers a load-balancing system with stochastically drawn connections
between servers. Jobs in the load balancing system arrive at a server site,
where each server similarly acts as a dispatcher and keeps or forwards the job
according to the JSQ($2$) policy based on its connected neighbors. The second
example illustrates the application of deterministic sampling for a bike
sharing system, with particles being stations and the graphon determining the
popularity of the stations. In both cases, a simple discretized version of the
integro-differential equation already yields precise estimations of the system
dynamics while the numerical complexity of the approximation only slightly
increases compared to the homogeneous case.

###  1.2 Organization

The paper is organized as follows. In Section 2 we introduce the heterogeneous
particle model. Section 3 defines graphon, related sampling methods, and
related preliminary results. Our definition of the approximation, the main
results, and the proofs are displayed in Section 4. Last, in Section 5 we
present the two numerical examples.

###  1.3 Related Work

##### Dynamical Systems on Random Graphs

In recent years there has been growing interest in the behavior of interacting
particles that are interconnected by an underlying graph topology, for
example, Bayraktar and Wu (2021), Bet et al. (2024), Bhamidi et al. (2019),
Abbe (2018) and previously mentioned references. For a general introduction to
the topic of random graph networks and limiting graphon functions for dense
graphs, we refer to the works Van Der Hofstad (2017), Lovász (2012). The
majority of papers focus on dense graphs, having edges of order $N^{2}$ such
as Erdős–Rényi type graphs or graphs generated by stochastic block models.
Additional related work can be found in the game theoretic setting for graphon
mean field games, see for example Caines and Huang (2021), Aurell et al.
(2022). In the ‘not-so-dense’ setting available results are more limited, with
some newer references being Bayraktar et al. (2023), Delmas et al. (2023). In
the case of sparse graphs, such as $d$-regular graphs or random geometric
graphs, the typical mean field methods break down as they fail to capture the
importance of the spatial graph structure and its implications on the local
dynamics of particles. Some recent works in this setting include Ganguly and
Ramanan (2022), Ramanan (2022), Ganguly (2022).

##### Load Balancing on Graphs

Our load balancing example is inspired by the recent works of Rutten and
Mukherjee (2023), Zhao and Mukherjee (2024), Zhao et al. (2024), Budhiraja et
al. (2019). Here, the authors study a variety of load balancing systems with
dynamics based on compatibility or locality constraints, which give rise to
intricate connectivity between dispatchers and servers. While not directly
transferable into the framework of this paper, the authors similarly deal with
graph-based systems and limit approximations that are strongly related to the
ones our framework suggests. Note that the techniques developed in these
papers are very model-specific and allow for transitions to depend on the
states of multiple queues, whereas our approach aims at deriving results for a
more general framework for transition rates with the restriction to the case
of pair interactions.

##### Generator and Stein’s Method

For our proofs, we adapt techniques used in Gast and Van Houdt (2017),
Allmeier and Gast (2022), Gast (2017), which in turn rely on the use of
Stein’s method Stein (1986). The method is used to estimate and bound the
distance between two random variables through their respective generators.
Since the works Braverman et al. (2017), Braverman and Dai (2017) Stein’s
method has seen an increase in the stochastic network community and is an
actively evolving area.

##  2 Heterogeneous Network Particle System

###  2.1 The Interaction Model

We consider particle systems with $N\in\mathbb{N}$ interacting particles.
Connections between particles are characterized by a (possibly weighted)
adjacency matrix $G^{N}\in[0,1]^{N\times N}$, with $G^{N}_{kl}$ indicating the
connection strength between particle $k$ and $l$. Each of the particles has a
finite state space $\mathcal{S}$ where the state of the $k$-th object at time
$t\geq 0$ is denoted by $S^{G^{N}}_{k}(t)$. As we see later, $G^{N}$ can
correspond to a random graph for which $G^{N}_{kl}\in\\{0,1\\}$ indicates the
presence or absence of an undirected edge between the particles $k$ and $l$ or
can be an arbitrary weight matrix, see Section 3. The state of the whole
system is denoted by
$\BFS^{G^{N}}(t):=\bigl{(}S^{G^{N}}_{1},\dots,S^{G^{N}}_{N}\bigr{)}(t)\in\mathcal{S}^{N}$.
We assume that the process $\BFS^{G^{N}}:=(\BFS^{G^{N}}(t))_{t\geq 0}$ is a
continuous time Markov chain (CTMC) with the dynamics of the system described
as in the following. Each particle $k\in[N]$ changes its state from $s_{k}$ to
$s_{k}^{\prime}\in\mathcal{S}$ in one of the two ways:

  1. (Unilateral)

The change of state occurs at rate $r_{k,s_{k}\to
s_{k}^{\prime}}^{N,\text{uni}}$ independently of the other particles.

  2. (Pairwise)

The change of state is triggered by another particle $l\in[N]$ that is in
state $s_{l}\in\mathcal{S}$. This occurs at rate $r_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}^{N,\text{pair}}G^{N}_{kl}/N$.

Note that the rate functions are assumed to be heterogeneous, i.e., they can
depend on the items $k$ and $l$. The rates have a $1/N$ factor, as each
particle can potentially interact with all $N-1$ remaining particles. Hence,
our condition implies that the total rate of transitions of the particle
system is of order $O(N)$ and that the transition rates for all particles are
of the same order. As we will further discuss in Example 5.1, our results can
also be used if the scaling factor depends not on the system size $N$ but on
the node degrees. We further want to point out that we restrict our framework
to the interaction of two particles, which can be utilized to model many
relevant interacting particle systems on graphs such as _e.g._ , epidemic
spreading, power-of-two-choices load-balancing, or bike sharing systems. It is
nonetheless possible to use the same underlying approach to extend the
framework and results to interactions of higher order. This, however, comes at
the cost of increasingly cumbersome notations and limited added theoretical
insight.

###  2.2 The Binary State Representation

In order to ease computations and definitions, we will use a binary
representation of the state based on indicator functions. We denote the new
representation by
$\BFX^{G^{N}}=(X^{G^{N}}_{k,s}(t))_{\begin{subarray}{c}k\in[N]\
s\in\mathcal{S}\\\ t\geq 0\end{subarray}}$, where

$\displaystyle
X^{G^{N}}_{k,s}(t):=\mathds{1}_{\\{S^{G^{N}}_{k}(t)=s\\}}:=\begin{cases}1&\text{if
object $k$ is in state $s$ at time $t$},\\\ 0&\text{otherwise.}\end{cases}$

The space of attainable states is denoted by
$\mathcal{X}^{N}\subset\\{0,1\\}^{N\times\mathcal{S}}$.

While this representation is less compact than the original, it allows for an
easier definition of transition rates as well as the definition of the mean
field approximation. Denote by $\BFe^{N}_{k,s}$ a matrix of size
$N\times\left|\mathcal{S}\right|$ whose $(k,s)$ component is equal to 1, all
other entries being zero. For each $k\in[N]$, and
$s_{k},s_{k}^{\prime}\in\mathcal{S}$, $X^{G^{N}}$ jumps to
$X^{G^{N}}+\BFe_{k,s_{k}}^{N}-\BFe_{k,s_{k}^{\prime}}^{N}$ at rate

$\displaystyle X_{k,s_{k}}r_{k,s_{k}\to
s_{k}^{\prime}}^{N,\text{uni}}+X_{k,s_{k}}\sum_{l\in[N]}\sum_{s_{l}\in\mathcal{S}}r_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}^{N,\text{pair}}\frac{G^{N}_{kl}}{N}X_{l,s_{l}}.$ (1)

In the above equation, the first term of the rate corresponds to the
unilateral transition of the particle $k\in[N]$ changing its state from
$s_{k}$ to $s_{k}^{\prime}$, as this transition occurs at intensity
$r_{k,s_{k}\to s_{k}^{\prime}}$ if particle $k$ is in state $s_{k}$ (_i.e._ ,
$X_{k,s_{k}}=1$). The second term describes the pairwise transitions leading
to the state change of particle $k$ from $s_{k}$ to $s_{k}^{\prime}$. Similar
to the unilateral one, the transition can only happen if particle $k$ is in
state $s_{k}$, represented in the rate by the prefactor $X_{k,s_{k}}$. The
remaining factor corresponds to the interaction with other particles,
expressed by the weighted sum over all other particles and their states. The
intensity of the transition is scaled by $r_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}^{N,\text{pair}}$ and the connectivity of the particles
is given by the connectivity matrix $G^{N}\in[0,1]^{N\times N}$.

###  2.3 Drift of the System of Size $N$

By using the state representation $\BFX^{G^{N}}$, we define what we call the
_drift_ of the system of $N$ particles as the expected change for
$\BFX^{G^{N}}$ in state $\BFX\in\mathcal{X}^{(N)}$. It is equal to the sum of
all possible transitions of the changes induced by this transition times the
rate at which the transition occurs. We denote this quantity as
$\BFF^{G^{N}}(\BFX)$. By using Equation (1), it is equal to:

$\displaystyle\BFF^{G^{N}}(\BFX)$
$\displaystyle=\sum_{k\in[N],s_{k},s_{k}^{\prime}\in\mathcal{S}}(\BFe^{N}_{k,s_{k}^{\prime}}-\BFe^{N}_{k,s_{k}})\left[X_{k,s_{k}}r_{k,s_{k}\to
s_{k}^{\prime}}^{N,\text{uni}}+X_{k,s_{k}}\sum_{l\in[N]}\sum_{s_{l}\in\mathcal{S}}r_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}^{N,\text{pair}}\frac{G^{N}_{kl}}{N}X_{l,s_{l}}\right]$

The quantity $\BFF^{G^{N}}(\BFX)$ is a vector-valued function of $\BFX$. By
reorganizing the above sum, $F^{G^{N}}_{k,s}(\BFX)$ –its $(k,s)$ component– is
equal to

$\displaystyle
F^{G^{N}}_{k,s_{k}}(\BFX)=\sum_{s^{\prime}\in\mathcal{S}}\left[X_{k,s^{\prime}}r_{k,s_{k}^{\prime}\to
s_{k}}^{N,\text{uni}}-X_{k,s}r_{k,s_{k}\to
s_{k}^{\prime}}^{N,\text{uni}}+\sum_{l\in[N]}\sum_{s_{l}\in\mathcal{S}}(X_{k,s_{k}^{\prime}}r_{k,l,s_{k}^{\prime}\to
s_{k},s_{l}}^{N,\text{pair}}-X_{k,s}r_{k,l,s_{\to}s_{k}^{\prime},s_{l}}^{N,\text{pair}})\frac{G^{N}_{kl}}{N}X_{l,s_{l}}\right].$

In the following, it will be convenient to replace the above sum by matrix
multiplications. To do so, let us denote by $X_{k}$ the vector
$(X_{k,s})_{s\in\mathcal{S}}$. The above equation can be written as:

$\displaystyle
F^{G^{N}}_{k,s}(\BFX)=\BFR^{N,\text{uni}}_{k,s}\BFX_{k}+\BFX^{T}_{k}\sum_{l\in[N]}\BFR^{N,\text{pair}}_{k,l,s}\BFX_{l}\frac{G^{N}_{kl}}{N},$
(2)

where $\BFR^{N,\text{uni}}_{k,s}$ is a row vector whose $s^{\prime}$ component
is $r_{k,s^{\prime}\to s}^{N,\text{uni}}$ if $s^{\prime}\neq s$ and
$-\sum_{\tilde{s}}r_{k,s\to\tilde{s}}^{N,\text{uni}}$ for $s=s^{\prime}$; and
$\BFR^{N,\text{pair}}_{k,l,s}$ is a matrix whose $(s^{\prime},s_{l})$
component is $r_{k,l,s^{\prime}\to s,s_{l}}^{N,\text{pair}}$ if
$s^{\prime}\neq s$ and
$-\sum_{\tilde{s}}r_{k,l,s\to\tilde{s},s_{l}}^{N,\text{pair}}$ if
$s=s^{\prime}$.

###  2.4 Representation of $\BFX^{G^{N}}$, $G^{N}$ and $\BFF^{G^{N}}$ as
Functions from $(0,1]$

To study the limit as $N$ goes to infinity, it will be convenient to view the
functions $\BFX^{G^{N}}$ not as a vector with $N$ components
$(X^{G^{N}}_{k})_{k\in[N]}$ but as a step function $(X^{N}_{u})_{u\in(0,1]}$,
where for any state $s\in\mathcal{S}$, we set

$\displaystyle X^{G^{N}}_{u,s}:=X^{G^{N}}_{k,s}\in\\{0,1\\}\text{ for
}u\in((k-1)/N,k/N].$

By abuse of notation, we do not introduce separate notations for the discrete
and continuous variables but make the distinction by reserving the subscript
letters $k,l\in[N]$ for the discrete case and $u,v\in(0,1]$ for the continuous
variable.

Similarly, we also write $G^{N}_{uv}=G^{N}_{kl}$ for $u\in((k-1)/N,k/N]$ and
$v\in((l-1)/N,l/N]$, $F^{G^{N}}_{u,s}=F^{G^{N}}_{k,s}$, and
$\BFR^{N,\mathrm{uni}}_{k,s}=\BFR^{N,\mathrm{uni}}_{u,s}$. By using this
notation, the sum of $k\in[N]$ –for instance in (2)– can be replaced by an
integral, i.e.,
$F^{G^{N}}_{k,s}(\BFX)=F^{G^{N}}_{u,s}(\BFX)=\BFR^{N,\text{uni}}_{u,s}\BFX_{u}+\BFX^{T}_{u}\int_{0}^{1}\BFR^{N,\text{pair}}_{u,v,s}\BFX_{v}G^{N}_{uv}dv$
for $u\in((k-1)/N,k/N]$.

###  2.5 Notations

Throughout the paper, matrices, and vectors are written in bold letters, i.e.
$\BFX,\BFx,\dots$, and regular letters are used to denote scalars like
$X_{u,s},r_{u,s\to s^{\prime}}^{\text{uni}}$. The indices
$s,s_{k},s_{k}^{\prime},s_{l},...$ are reserved for the states, $k,l,...$ are
reserved for particles, and $u,v$ are reserved for values in the unit
interval. When we write that a quantity $h$ is of order $\mathcal{O}(1/N)$ or
$h=\mathcal{O}(1/N)$ equivalently, this means that there exists a constant $C$
such that $h\leq\frac{C}{N}$. Most of our results will be expressed in terms
of $L_{2}$ norm (see the definition in Section 3.3). The space $L_{2}(0,1]$
refers to the quotient space of the square Lebesgue integrable functions.
Throughout we will deal with vectors of $L_{2}(0,1]$ functions, i.e., $f,g\in
L_{2}(0,1]^{\mathcal{S}}$. We for the vectors $g,f$ we denote the scalar
product by $\langle
f,g\rangle_{L_{2}(0,1]}=\sum_{s\in\mathcal{S}}\int_{0}^{1}f_{u,s}g_{u,s}\ du$
for and induced norm
$\left\|f\right\|_{L_{2}}=\sqrt{\sum_{s\in\mathcal{S}}\int_{0}^{1}f_{u}^{2}du}$.
For a function $G:(0,1]^{2}\to\mathbb{R}$ denote by
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}=\sup_{\\{f\in
L_{2}(0,1],\left\|f\right\|_{L_{2}(0,1]}\leq
1\\}}\sqrt{\int_{0}^{1}(\int_{0}^{1}G_{uv}f_{v}\ dv)^{2}du}$ the $L_{2}$
operator norm.

##  3 Limiting Graph and Graphon

In this section, we specify the properties that the interaction graph $G^{N}$
needs to satisfy as $N$ goes to infinity. To do so, we introduce the notion of
graphon and define the associated cut-norm. We also introduce two sampling
methods that can be used to generate a graph with $N$ nodes from a graphon.
This section only reviews the material that is necessary for our results, and
we refer to the famous book Lovász (2012) for a detailed introduction to
graphons.

###  3.1 Graphons

In this paper, what we call a graphon is a measurable function111In the
literature, graphons are often restricted to symmetric functions, which
correspond to undirected graphs. This restriction is not needed for our case.
$G:(0,1]^{2}\to(0,1]$. The notion of graphon can be viewed as a generalization
of the notion of graph. Indeed, for any $N$, a weighted graph $G^{N}$ can be
viewed as a piecewise constant function defined on $(0,1]^{2}$, where the
value of this function at a point $(u,v)\in(0,1]^{2}$ is equal to $G^{N}_{kl}$
whenever $u\in((k-1)/N,k/N]$ and $v\in((l-1)/N,l/N]$. Hence, a finite graph is
a graphon that has a special structure. The notion of graphon generalizes the
notion of finite graph by allowing $G$ to be any measurable function. We
provide an illustration of a graphon and of a finite graph viewed as a graphon
on Figure 1.

Throughout the paper, we consider piecewise Lipschitz continuous graphons,
which are defined as follows.

######  Definition 3.1 (Piecewise Lipschitz Graphon)

A graphon $G$ is called piecewise Lipschitz if there exists a constant $L_{G}$
and a finite partition of $(0,1]$ of non-overlapping intervals
$\mathcal{A}_{k}=(a_{k-1},a_{k}]$ with $0=a_{0}<a_{1}<...<a_{K_{G}}=0$ for a
finite $K_{G}\in\mathbb{N}$, such that for any $k_{1},k_{2}\in[K+1]$ and pairs
$(u,v),(u^{\prime},v^{\prime})\in\mathcal{A}_{k_{1}}\times\mathcal{A}_{k_{2}}$

$\displaystyle\left|G_{uv}-G_{u^{\prime}v^{\prime}}\right|\leq
L_{G}(\left|u-u^{\prime}\right|+\left|v-v^{\prime}\right|).$

A particular case of a piecewise Lipschitz continuous graphon is the case of a
step graphon, that is such that the Lipschitz constant $L_{P}=0$. This implies
that the function $G$ is constant on all intervals
$\mathcal{A}_{i}\times\mathcal{A}_{j}$. For instance, all finite graphs can be
seen as step graphons by using the partition such that
$\mathcal{A}_{k}=(k-1/N,k/N]$ with $N$ being the number of nodes in the graph.

Stochastic Graph Sampling, Graphon discretization, Graphon Illustrations

(a) Graphon corresponding to $G_{uv}=\frac{1}{2}(1-\left|u-v\right|)$

(b) Step graphon associated to a sampled graph with $N=50$ nodes.

(c) Step graphon obtained from graphon discretization for $N=16$.

Figure 1: Exemplary connectivity functions sampled after the methods of
Section 3.2 and their corresponding graphon.

###  3.2 Generation of a Finite Graph $G^{N}$ from a Graphon $G$

Based on a given graphon $G$, we consider two distinct sampling methods to
generate a finite graph $G^{N}$ from $G$:

  * •

In the first case, termed _deterministic sampling_ , $G^{N}$ is a
discretization of $G$ on $N^{2}$ uniformly sampled points, i.e.,
$G^{N}_{kl}=G(k/N,l/N)$ for $k,l\in[N]$.

  * •

For the second case, termed _stochastic sampling_ , and under the preliminary
assumption that the graphon is symmetric, the values $G^{N}$ are drawn
according to independent Bernoulli random variables, i.e.,
$G^{N}_{kl}=\mathrm{Bernoulli}\bigl{(}G(k/N,l/N)\bigr{)}$, and
$G^{N}_{lk}=G^{N}_{kl}$ for all $k<l$.

In the case that the graphon is constant for all $u,v\in(0,1]$, the stochastic
sampling method is equivalent to sampling an Erdős–Rényi graph. If the graphon
is block-wise constant, the sampling is similar to the stochastic block model;
see Van Der Hofstad (2017).

###  3.3 Graphon Distances and Convergence

To measure the distance between two graphs (and in particular to quantify how
fast does a finite graph $G^{N}$ converge to $G$), we will use the $L_{2}$
operator norm. More precisely, for a measurable function
$f:(0,1]\to\mathbb{R}$, we denote by
$\left\|f\right\|_{L_{2}(0,1]}=\sqrt{\int_{0}^{1}f_{u}^{2}du}$ the $L_{2}$
norm of $f$. For a graphon $G$, we denote by
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|$
the operator norm of $G$ in $L_{2}$, that is:

$\displaystyle\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}=\sup_{\\{f\in
L_{2}(0,1]\text{ such that }\left\|f\right\|_{L_{2}(0,1]}\leq
1\\}}\left\|Gf\right\|_{L_{2}}.$

The distance between two graphons (or finite graphs) $G$ and $G^{\prime}$ is
measured as
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G-G^{\prime}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|$.

To bound the $L_{2}$ distance between a piecewise Lipschitz graphon and the
function associated with a randomly sampled graph, we utilize the results
displayed in Avella-Medina et al. (2018). The authors show that for a large
enough (see Definition 3.3) number $N$ of graph nodes, in our case particles,
the distance between the graphon and the graph function can be upper bounded
as follows:

######  Lemma 3.2 (Theorem 1 from Avella-Medina et al. (2018))

Let $G$ be a symmetric piecewise Lipschitz graphon, and let $G^{N}$ be a
stochastic sampling of it as defined in 3.2. Then, for large enough $N$ (as in
Definition 3.3) with probability as least $1-\delta$ the distance in the
$L_{2}$ operator norm between the graphon $G$ and the sampled graph $G^{N}$ is
bounded by

$\displaystyle\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G-G^{N}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}\leq\sqrt{\frac{4\log{(2N/\delta)}}{N}}+2\sqrt{\frac{(L^{2}_{G}-K^{2}_{G})}{N^{2}}+\frac{K_{G}}{N}}=:\psi_{\delta,G}(N).$
(3)

######  Definition 3.3 (Large Enough $N$ Avella-Medina et al. (2018))

Given a piecewise Lipschitz graphon $G$ with partition of $(0,1]$ of non-
overlapping intervals $\mathcal{A}_{k}=(a_{k-1},a_{k}]$ as in Definition 3.1
and $\delta<e^{-1}$. The quantity $N$ is called ‘large enough’, if

$\displaystyle\frac{2}{N}<\min_{k\in\\{0..K_{G}\\}}(a_{k}-a_{k-1}),$ (4a)
$\displaystyle\frac{1}{N}\log\left(\frac{2N}{\delta}\right)+\frac{2K_{G}+3L_{G}}{N}<\sup_{u\in(0,1]}\int_{0}^{1}G_{u,v}\
dv\quad\text{and}$ (4b) $\displaystyle Ne^{-N/5}<\delta.$ (4c)

Implied by Theorem 3.2, by setting $\delta_{N}=\frac{2}{N}$, we see that the
distance between a piecewise Lipschitz graphon and a sampled graph is with
high probability of order
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G-G^{N}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}(0,1]}=\mathcal{O}\left(\sqrt{\frac{\log(N)}{N}}\right)$.
The meaning of ‘with high probability’ in this context is that the right-hand
side holds with probability at least $1-2/N$.

##  4 Main Results

###  4.1 Assumptions

In Section 2, we constructed a model that depends on a scaling parameter $N$.
We state the necessary assumptions that will be used to state the accuracy of
the graphon mean field approximation.

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

The state space $\mathcal{S}$ is finite.

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

The graphon $G$ is piecewise Lipschitz continuous.

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

There exists bounded and piecewise Lipschitz-continuous rate functions for
$r_{u,s_{u}\to s_{u}^{\prime}}^{\text{uni}}$ and $r_{u,v,s_{u}\to
s_{u}^{\prime},s_{v}}^{\text{pair}}$ for $u,v\in(0,1]$ and
$s_{u},s_{u}^{\prime},s_{v}\in\mathcal{S}$ such that the rates function of the
original $N$ particle systems have the relation:

$\displaystyle r_{k/N,s_{k}\to s_{k}^{\prime}}^{\text{uni}}=r_{k,s_{k}\to
s_{k}^{\prime}}^{N,\text{uni}}\quad\text{ and }\quad r_{k/N,l/N,s_{k}\to
s_{k}^{\prime},s_{l}}^{\text{pair}}=r_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}^{N,\text{pair}}\quad\text{ for }k,l\in[N].$ (5)

The first assumption is technical and simplifies the definition of the $L_{2}$
space. The second assumption is very classical regularity assumption when
studying sequences of graphs that converge to graphons Lovász (2012). In
practice, our bounds will depend on the distance between the original graph
$G^{N}$ and the graphon $G$. The last assumption ensures that the particle
transition rates for the finite and the graphon system are equal. We point out
that assumption ($A_{\text{3}}$) can be generalized by replacing the equality
in (5) with bounds for the distance of $r^{N,\cdot}$ and $r^{\cdot}$ therefore
modifying the bounds of the theorem to include terms of the form
$\left\|r^{N,\cdot}-r^{\cdot}\right\|$.

###  4.2 The Graphon Mean Field Approximation

The graphon mean field approximation aims at approximating the dynamics of the
original system $\BFX$. We define the graphon related drift similarly to the
drift of the particle system in Equation (2). For $u,v\in(0,1]$ the graphon
based drift is defined by

$\displaystyle
F^{G}_{u,s}(\BFx)=\BFR^{\text{uni}}_{u,s}\BFx_{u}+\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFx_{v}dv.$
(6)

This equation is identical to the original drift equation (2) with two
modifications: First, the rates do not depend on $N$. Second, the discrete
variables $k,l\in[N]$ are replaced by continuous variables $u,v\in(0,1]$,
which notably replaces the sum over $l$ by an integral over $v$.

Based on the drift function $\BFF^{G}$ and the initial condition $\BFx_{0}$,
we call the graphon mean field approximation the solution of the differential
equation

$\displaystyle\BFx^{G}(t,\BFx_{0})=\BFx_{0}+\int_{0}^{t}\BFF^{G}(\BFx^{G}(\tau,\BFx_{0}))d\tau.$
(7)

###### Lemma 4.1

Let $\BFF^{G}$ denote the deterministic drift defined in Equation (6) and
$\BFx^{G}(t,\BFx_{0})$ the solution of the graphon mean field approximation at
time $t$ with initial condition $\BFx_{0}$ as in Equation (7). It holds that
$\BFx^{G}$ is well-defined, has a unique solution and is differentiable with
respect to its initial condition. Furthermore, this derivative is Lipschitz
continuous.

The proof is postponed to Section 6.1.

###  4.3 Accuracy of the Approximation

The following result provides a bound on the distance between the stochastic
system and the graphon mean field approximation. This bound is stated as the
$L_{2}$ distance between $\mathbb{E}[\BFX^{G^{N}}]$ and $\BFx^{G}$. Recall
that by definition, $\mathbb{E}[X^{G^{N}}_{k,s}(t)]=\mathbb{P}(S_{k}(t)=s)$ is
the probability for an item $k$ to be in state $s$ at time $t$. Hence, our
theorem shows that the graphon mean field model is an accurate approximation
of the state distribution of the particles. The statement of the theorem
should be interpreted as saying that the distribution of the particles over
the states $\mathcal{S}$ for any time $t\geq 0$ is approximated by the graphon
mean field approximation $\BFx^{G}$ with accuracy depending on the system size
$N$ as well as the distance between $G^{N}$ and graphon $G$. In particular, if
$G^{N}$ converges to $G$ (for the graphon norm), then the graphon mean field
approximation is asymptotically exact.

######  Theorem 4.2 ($L_{2}$ convergence)

Let $\BFX^{G^{N}}(t)=(X^{G^{N}}_{k,s})(t)_{k\in[N],s\in\mathcal{S}}$ be the
stochastic particle system of size $N$ related to a graph instance $G^{N}$.
Let $\BFx^{G}(t)=(x^{G}_{u,s}(t,\BFx_{0}))_{u\in(0,1],s\in\mathcal{S}}$ be the
mean field approximation of the particle system as defined in Equation (7) for
an initial condition $\BFx_{0}=\BFX^{G^{N}}(0)$222Here, $\BFX^{G^{N}}(0)$ is
interpreted as a vector in $L_{2}(0,1]^{\left|\mathcal{S}\right|}$.. Assume
additionally that conditions ($A_{\text{1}}$) \- ($A_{\text{3}}$) are
fulfilled and let $t>0$ be arbitrary but fixed. Then, there exist constants
$C_{A},C_{B}\geq 0$ such that

$\displaystyle\left\|\mathbb{E}[\BFX^{G^{N}}(t)\mid\BFX^{G^{N}}(0),G^{N}]-\BFx^{G}(t,\BFX^{G^{N}}(0))\right\|_{L_{2}}\leq\frac{C_{A}}{N}+C_{B}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|.$
(8)

The proof of Theorem 4.2 is postponed to Section 6.2. The constants
$C_{A},C_{B}$ of Equation (8) depend on the uniform bound of the rates, the
time $t$ and the Lipschitz constants of the graphon. In the subsequent
Corollaries 4.3 and 4.4 we will see, that if the $G^{N}$ is generated by one
of the methods illustrated in Section 3.2, precise bounds on the distance
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|$
can be obtained. To illustrate our results and underline that the constants
are small in practice, we provide examples in Section 5. We point out that
Theorem 4.2 is applicable for a wide range of construction methods for
$G^{N}$. The subsequent corollaries illustrate cases of stochastic and
deterministic sampling methods. We emphasize, however, that any method that
allows the construction of densely connected graphs, for which bounds of
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|$
are attainable, is viable. At last, we want to point out that by using the
same framework, it seems feasible to extend our results to interactions of
triplets or higher order interactions seems feasible. Yet, due to our generic
choice of transition rates, this would be linked to increasingly heavy
notations for the dynamics while only giving little more insight into the
accuracy of the graphon mean field method.

###  4.4 Case Specific Bounds for
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}$

The subsequent corollaries, give specific bounds for the case that $G^{N}$ was
generated as described in Section 3.2. In the first case, if $G^{N}$ is
obtained though discretization of $G$, Corollary 4.3 shows that the accuracy
is of order $\mathcal{O}(1/N)$. Our second Corollary 4.4 specifies
big-$\mathcal{O}$ convergence rates if the graph $G^{N}$ is sampled
stochastically from the graphon. In this case, the accuracy is with high
probability of order $\mathcal{O}(\sqrt{\frac{\log{N}}{N}})$.

####  4.4.1 Case 1: Graphon Discretization

The corollary give accuracy bounds for the approximation in the case that the
graph $G^{N}$ is obtained as a discretized version of $G$. The result gives
bounds on the difference between the distributions of the particles in the
stochastic system and the approximate values obtained through the graphon mean
field approximation.

######  Corollary 4.3 (Dense Heterogeneous System)

Assume ($A_{\text{1}}$) \- ($A_{\text{3}}$) and let $\BFx^{G}$ and
$\BFX^{G^{N}}$ be defined as in Theorem 4.2. Let $k\in[N],s\in\mathcal{S}$ and
$t\geq 0$ be arbitrary but fixed. If $G^{N}$ is generated by the deterministic
sampling method of 3.2, i.e., a discretization of the graphon $G$, it holds
that

$\displaystyle\left\|\mathbb{E}[\BFX^{G^{N}}(t)\mid\BFX^{G^{N}}(0),G^{N}]-\BFx^{G}(t,\BFX^{G^{N}}(0))\right\|_{L_{2}}$
$\displaystyle\leq\frac{C_{A}+C_{B}C_{G^{N}}}{N}$ (9)

The proof is postponed to Section 6.3. The constants $C_{A},C_{B}$ are as in
Theorem 4.2. The additional constant $C_{G^{N}}$ relates to the discretization
error of the deterministic sampling method. For this case, it is noteworthy
that the accuracy of the approximation aligns with the results one obtains for
the homogeneous setting as described in Gast (2017), Gast and Van Houdt (2017)
while allowing for heterogeneous connectivity and rates among the population.

####  4.4.2 Case 2: Random Graph

Our second corollary provides accuracy bounds for interacting particle systems
on dense random graphs. By the definition of the graph sampling methods, see
Section 3.2, with high probability, the number of connected neighbors for each
particle is of order $N$. This ensures that for large enough systems that the
neighborhood of each node serves as a local representation of the overall
system state. This leads to the following result:

######  Corollary 4.4 (Graphon System Approximation)

Assume the conditions ($A_{\text{1}}$) \- ($A_{\text{3}}$) as in Theorem 4.2
for a symmetric graphon $G$. Let
$\BFX^{G^{N}}(t)=(X^{G^{N}}_{k,s})(t)_{k\in[N],s\in\mathcal{S}}$ be a
stochastic particle system of size $N$ with $G^{N}$ obtained through the
stochastic sampling method as defined in 3.2. Let $t\geq 0$ and
$k\in[N],s\in\mathcal{S}$ be arbitrary but fixed. Then, with probability at
least $1-2/N$ and for ‘large enough’ $N$, as defined in Definition 3.3, the
graph $G^{N}$ is sampled such that

$\displaystyle\left\|\mathbb{E}[\BFX^{G^{N}}(t)\mid\BFX^{G^{N}}(0),G^{N}]-\BFx^{G}(t,\BFX^{G^{N}}(0))\right\|_{L_{2}}\leq\frac{C_{A}}{N}+C_{B}\psi_{G}(N)$
(10)

where
$\psi_{G}(N):=\sqrt{\frac{8\log{(N)}}{N}}+2\sqrt{\frac{(L_{G}^{2}-K_{G}^{2})}{N^{2}}+\frac{K_{G}}{N}}$
with $L_{G}$ being the Lipschitz constant of the graphon $G$ and $K_{G}$ the
size of the partition as defined in 3.1.

The proof of the corollary is postponed to Section 6.3.

##  5 Numerical Experiments

In this section, we present two examples which support the statements of our
theoretical results and illustrate the applicability of the framework. For the
first example, we look at a load balancing model with communication
restrictions of the servers imposed by a graph. In this example, the focus is
on the stochastically sampled graph, which imposes heterogeneous rates due to
connectedness of the servers. For the dynamics of the system, we see each node
as a dispatcher / server pair applying the JSQ($2$) policy with respect to
itself and connected servers whenever a job arrives. For the second example,
we consider a heterogeneous bike-sharing system. Here, the focus lies on the
heterogeneity of the popularity of the stations, which affects the flow of
bikes through the system.

###  5.1 Load Balancing Example

####  5.1.1 Model

The considered load balancing model is a jump process in the Markovian
setting. We consider a system where the servers-dispatcher pairs are connected
through a graph structure where each server-dispatcher is represented by a
node. All servers have a finite maximal queue length $K_{L}\in\mathbb{N}$. The
connections between the server-dispatcher pairs are denoted by $G^{N}_{kl}$
and are sampled according to Section 3.2 using the stochastic sampling method.
The graphon $G$ used to sample the edges is the same as Figure 1(a). Based on
the described connectivity structure, jobs arrive to a server-dispatch pair
following a Poisson arrival process with rate $\lambda_{L}>0$. Arriving jobs
are distributed according to the JSQ($2$) policy idea, i.e., the dispatcher
considers its own server state as well as another randomly sampled but
connected server and forwards the job to the server with the lesser load. In
the case that both servers have the same queue length, the job is assigned
randomly among the two. In case both servers have a full queue, the job is
discarded. The service time of a job is exponentially distributed with mean
$\mu_{L}$ and jobs are handled in first come, first served order. For a system
of size $N$ with graph instance $G^{N}$ and state
$\BFX^{G^{N}}=(X^{G^{N}}_{k,s})_{k\in[N],s=0..K_{L}}$ the transitions of the
system are

$\displaystyle\BFX^{G^{N}}\to\BFX^{G^{N}}+\BFe_{k,s+1}^{N}-\BFe_{k,s}^{N}$
(11) $\displaystyle\quad\text{ at rate
}\lambda_{L}X^{G^{N}}_{k,s}\sum_{l\in[N]}(\frac{G^{N}_{kl}}{d^{(N)}(k)}+\frac{G^{N}_{kl}}{d^{(N)}(l)})\Bigl{(}\frac{1}{2}X_{l,s}+\sum_{s_{l}\geq
s+1}X_{l,s_{l}}\Bigr{)}\mathds{1}_{s<K_{L}}$ (12)
$\displaystyle\BFX^{G^{N}}\to\BFX^{G^{N}}+\BFe_{k,s-1}^{N}-\BFe_{k,s}^{N}\quad\text{
at rate }X^{G^{N}}_{k,s}\mu_{L}\mathds{1}_{s>0}$ (13)

where $\BFe_{k,s}$ is a $N\times\left|\mathcal{S}\right|$ matrix having a one
at the $(k,s)$ entry and zero values everywhere else and $d^{(N)}(k)$ is the
degree of node $k$. Equation (13) corresponds to the completion of a job by
server $k$ when the queue is of size $1\leq s\leq K_{L}$. The second type of
transitions, equation (12), corresponds to adding a job to $k$ having $0\leq
s\leq K_{L}-1$ jobs in the buffer. In this case, the queue size is increased
by one from $s$ to $s+1$. To explain the transition rate we see that the
servers $k$ can be selected in two ways. By selecting $k$ or another connected
server $l$ first and the other second. In the case that both queues have equal
length, the chance that the job to arrives at server $k$ is $1/2$. If both
buffers are full, the job is discarded. In the case that a node associated to
server $l$ is isolated, i.e., has no edges, we define
$\frac{G^{N}_{kl}}{d^{(N)}(l)}$ to be zero.

####  5.1.2 Drift and Graphon Mean Field Approximation

For the deterministic drif, we replace the values of
$\frac{G^{N}}{d^{(N)}(k)}$ by the ones of the graphon $\frac{G}{d(v)}$ with
$d(v):=\int_{0}^{1}G_{v,\nu}d\nu$ and the sums over particles by an integral
over $(0,1]$. For a given state $\BFx=(x_{s})_{s=0,\dots,K_{L}}$ with
$x_{s}\in L_{2}(0,1]$, the drift evaluated at
$(u,s)\in(0,1]\times\\{0,\dots,K_{L}\\}$ is defined by

$\displaystyle
F^{G}_{u,s}(\BFx)=x_{u,s}\mathds{1}_{s<K_{L}}\int_{0}^{1}\lambda_{L}\left(\frac{G_{u,v}}{d(u)}+\frac{G_{u,v}}{d(v)}\right)\bigl{(}\frac{1}{2}x_{v,s}+\sum_{s_{v}\geq
s+1}x_{v,s_{v}}\bigr{)}dv\ -x_{u,s}\mu_{L}\mathds{1}_{s>0}.$

The graphon mean field approximation is then defined as in Equation (7).

####  5.1.3 Derivation of Accuracy Bounds

While edges between the server-dispatcher pair are sampled according to the
stochastic sampling method of Section 3.2, the dependence on of the rates on
the node degree, prevents a direct application of the results of Corollary
4.4. In particular, to apply Corrollary 4.4 it is assumed that the graph edges
are scaled by a factor of $1/N$ instead of the node degree. To resolve this
issue and obtain accuracy bounds similar to the one given by the Corollary, we
start by defining
$\tilde{G}^{N}_{kl}:=(\frac{G^{N}_{kl}}{d^{(N)}(k)}+\frac{G^{N}_{kl}}{d^{(N)}(l)})\frac{N}{16}$,
$r^{N,\text{pair}}_{k,l,s\to
s+1,s_{l}}=\lambda_{L}(\frac{1}{2}\mathds{1}_{s_{l}=s}+\mathds{1}_{s_{l}>s})\mathds{1}_{s<K_{L}}$,
and $r^{N,\text{uni}}_{k,s\to s-1}=\mu_{L}\mathds{1}_{s>0}$ with
$r^{N,\text{pair}}$ being the pairwise transitions and $r^{N,\text{uni}}$ the
uniateral transitions for $s,s^{\prime},s_{l}\in\\{0,...,K_{L}\\}$. This
allows to obtain the drift of the stochastic system similar to the one
introduced in Equation (2) namely

$\displaystyle F^{G^{N}}_{k,s}(\BFX)=X_{k,s}\sum_{s_{l}\geq
s}\sum_{l\in[N]}16\frac{\tilde{G}^{N}_{kl}}{N}r^{N,\text{pair}}_{k,l,s\to
s+1,s_{l}}X_{l,s^{\prime}}-X_{k,s}r^{N,\text{uni}}_{k,s\to s-1}.$

We proceed similarly for the drift of the graphon mean field approximation by
defining
$\tilde{G}_{uv}:=(\frac{G_{u,v}}{d(u)}+\frac{G_{u,v}}{d(v)})\frac{1}{16}$,
$r^{\text{pair}}_{k,l,s\to
s+1,s_{l}}:=\lambda_{L}(\frac{1}{2}\mathds{1}_{s_{l}=s}+\mathds{1}_{s_{l}>s})\mathds{1}_{s<K_{L}}$
and $r^{N,\text{uni}}_{k,s\to s-1}:=\mu_{L}\mathds{1}_{s>0}$ to rewrite

$\displaystyle F^{G}_{u,s}(\BFx)=x_{u,s}\int_{0}^{1}\sum_{s_{v}\geq
s}16\tilde{G}_{uv}r^{\text{pair}}_{k,l,s\to s+1,s_{v}}x_{v,s_{v}}dv-
x_{u,s}r^{N,\text{uni}}_{k,s\to s-1}.$

Note that for the above reformulations, the result of our main Theorem 4.2 is
still applicable and yields constants $C_{A},C_{B}\geq 0$ such that
$\left\|\mathbb{E}[\BFX^{G^{N}}(t)\mid
G^{N}]-\BFx^{G}(t)\right\|_{L_{2}}\leq\frac{C_{A}}{N}+C_{B}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\tilde{G}^{N}-\tilde{G}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}$.
Using the definition of the $\tilde{G}^{N},\tilde{G}$, we obtain
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\tilde{G}^{N}-\tilde{G}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}\leq\frac{1}{16}\Bigl{(}\left\|G^{N}-G\right\|_{L_{2}}\frac{N}{\min_{k}d^{N}(k)}+\sup_{u\in(0,1]}\left|\sum_{k=1}^{N}\frac{N}{d^{N}(k)}\mathbf{1}_{u\in(k-1/N,k/N]}-\frac{1}{d(u)}\right|\Bigr{)}$.
For
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}$,
we can use the bound as in Corollary 4.4 which is of order
$\mathcal{O}(\sqrt{\frac{\log N}{N}})$ with probability at least
$1-\frac{2}{N}$ and large enough $N$ (as in Definition 3.3). To obtain similar
bounds for the node degree, the multiplicative Chernoff bound can be used,
i.e.,
$\mathbb{P}\Bigl{(}\frac{N}{d^{N}(k)}\geq\frac{N}{(1-\gamma)\mathbb{E}[d^{N}(k)]}\Bigr{)}=\mathbb{P}\Bigl{(}d^{N}(k)\leq(1-\gamma)\mathbb{E}[d^{N}(k)]\Bigr{)}\leq\exp(-\frac{\mathbb{E}[d^{n}(k)]\gamma^{2}}{3})$.
Taking $\gamma=\sqrt{\frac{3\log{N}}{\mathbb{E}[d^{N}(k)]}}$ yields
$\mathbb{P}\Bigl{(}\frac{N}{d^{N}(k)}\geq\frac{N}{(1-\gamma)\mathbb{E}[d^{N}(k)]}\Bigr{)}\leq\frac{2}{N}$.
By the triangle inequality we have
$\left|\sum_{k=1}^{N}\frac{N}{d^{N}(k)}\mathbf{1}_{u\in(k-1/N,k/N]}-\frac{1}{d(u)}\right|\leq\left|\sum_{k=1}^{N}(\frac{N}{d^{N}(k)}-\mathbb{E}[\frac{N}{d^{N}(k)}])\mathbf{1}_{u\in(k-1/N,k/N]}\right|+\left|\sum_{k=1}^{N}\mathbb{E}[\frac{N}{d^{N}(k)}]\mathbf{1}_{u\in(k-1/N,k/N]}-\frac{1}{d(u)}\right|$.
We apply Chernoff bound the first summand to show that the difference
concentrates around zero, i.e.,
$\mathbb{P}\Bigl{(}\frac{N}{d^{N}(k)}\leq\frac{N}{(1-\gamma)\mathbb{E}[d^{N}(k)]}\Bigr{)}\geq
1-2/N$ and therefore $\frac{N}{d^{N}(k)}-\frac{N}{\mathbb{E}[d^{N}(k)]}\leq
N(\frac{1}{(1-\gamma)\mathbb{E}[d^{N}(k)]}-\frac{1}{\mathbb{E}[d^{N}(k)]})=\mathcal{O}(\sqrt{\frac{\log{N}}{N}})$
with probability at least $1-2/N$. The case
$\frac{N}{d^{N}(k)}\leq\frac{N}{(1-\gamma)\mathbb{E}[d^{N}(k)]}$ follows using
the same arguments. For the second summand we use the statement of Lemma A.7
to obtain an error bound of order $\mathcal{O}(1/N)$. This shows that for
large enough $N$,
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\tilde{G}^{N}-\tilde{G}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|_{L_{2}}$
is of order $\mathcal{O}(\sqrt{\frac{\log N}{N}})$.

####  5.1.4 Implementation

For the simulation, we set the parameters of the load balancing system as
follows: The arrival rate is set to $\lambda_{L}=1$, the server rate to
$\mu_{L}=1.1$ for all servers, and the maximal server capacity is $K_{L}=10$.
For each system size $N$, we sample a graph that remains fixed for all
simulations and obtain the sample mean by averaging over $2000$ sampled
trajectories for all system sizes. To compute the solution of the graphon mean
field approximation, we discretize the drift by a simple step discretization.
To be more precise, let $\gamma\in\mathbb{N}$ be the discretization parameter
and $U^{\gamma}=\\{U_{i}^{\gamma}\\}_{i=1..\gamma}$ be the even partition of
the interval $(0,1]$ with $U_{i}^{\gamma}:=(i-1/\gamma,i/\gamma]$. We define
by
$F^{\gamma}_{B}(\BFx)_{u,s}:=\sum_{i=1}^{\gamma}\mathds{1}_{U_{i}^{\gamma}}(u)F_{B}(\BFx)_{i/\gamma,s}$
the values of the discretized version of $F_{B}$. Throughout the numerical
experiments, we set $\gamma=100$. We choose this simple discretization method
as it provides good approximation results and low computation times, i.e., for
a rudimentary implementation using a NumPy ode solver, we are able to solve
the discretization in a few seconds.333For intricate intensity and graphon
functions it might be necessary to refine and tune the discretization
approach. The Figures 2 and 3 show the results of our numerical experiments.
The first Figure 2 shows the approximation accuracy for a single server and
different system sizes. We see that already for $N=40$ the sample mean is very
close to the approximation value, which supports the statement of Corollary
4.4. In the second Figure 3, we compare the values for the fraction of servers
with at least $s$ jobs for $s=1,...,4$ as described in the figure caption.
Remarkably here, even for small system sizes, the values obtained by the
sample mean are close to the ones obtained by the approximation. As we observe
the by Corollary 4.4 suggested behavior for single particles inb Figure 2, the
increased accuracy for $N=20,30$ is likely caused by an averaging effect and a
well-connected graph instance $G^{N}$.

Figure 2: The figure shows the evolution of
$\mathbb{P}(X^{G^{N}}_{k,s}(t)=1)=\mathbb{E}[X^{G^{N}}_{k,s}(t)]$, the
probability for a server-dispatcher pair to have $s=0,...,3$ jobs for
$t\in[0,2]$. In each plot the results for one system size $N=20,30,40$ are
displayed. Throughout, the sample mean values are plotted against the values
of the graphon mean field approximation.

Simulation and Approximation for the Load Balancing Example

Figure 3: The figure compares the values for the fraction of servers with at
least $s$ jobs obtained by the sample mean and the approximation. The
quantities are calculated as
$\mathbb{E}[Q_{s}(t)]=\mathbb{E}[\sum_{k\in[N]}\frac{1}{N}\sum_{s^{\prime}\geq
s}X^{G^{N}}_{k,s^{\prime}}]$ for the sample mean and
$q_{s}(t)=\int_{0}^{1}\sum_{s^{\prime}\geq s}x^{G}_{u,s^{\prime}}du$ for the
approximation.

###  5.2 Heterogeneous Bike-Sharing System

####  5.2.1 Model

We consider a bike-sharing model following the setup used in Fricker and Gast
(2016), Fricker et al. (2012). The model consists of $N\in\mathbb{N}$ bike
stations, representing the particles in the system. Each station has finite
capacity $K_{B}\in\mathbb{N}$. The system has a fleet of bikes of size
$M:=\lfloor\alpha N\rfloor$ which, at time $t=0$, is evenly distributed
amongst the stations, i.e., $\alpha$ bikes per station. The evolution of the
system depends on the movement of bikes. We differ between the two cases,
bikes in the system can either be stationary or in transit between stations.
To align with the notations used in the theory part of the paper, we denote by
$S_{k}^{N}(t),i\in[N]$ the number of bikes at station $i$ for time $t\geq 0$.
For the system, heterogeneity comes from the varying popularity of stations.
Hence, we denote by $p_{B}:(0,1]\mapsto\mathbb{R}_{\geq 0}$ the popularity
function. Based on the popularity function $p_{B}$, we define the graphon of
the bike sharing system by $G_{B}(u,v):=p_{B}(v)/\int_{0}^{1}p_{B}(\nu)d\nu$.
For the system of size $N$, connectivity for the stations is obtained by
deterministic sampling as described in Section 3.2, i.e., we discretize the
graphon based on the system size $N$. For two stations $k,l\in[N]$
connectivity is therefore defined by
$G^{N,B}_{kl}:=G_{B}(k/N,l/N)=p_{B}(l/N)/\int_{0}^{1}p_{B}(\nu)d\nu$. Based on
the connectivity between stations, it remains to define the dynamics of the
system. For a given state $(S_{1},...,S_{N})\in[K_{B}]^{N}$ the bikes move
between stations in the following way:

  * •

Customers arrive at a stations $k\in[N]$ with rate $\lambda_{B}$ leading to
the transitions

$\displaystyle S\to S-\BFe_{k}^{N}$ at rate
$\displaystyle\lambda_{B}\mathds{1}_{S_{k}>0}.\phantom{\frac{p_{B}(k/N)}{\int_{0}^{1}p_{B}(\nu)d\nu}(M-\sum_{k}S_{k}^{N})}$

  * •

The travel time of bikes between stations is exponentially distributed with
rate $\mu_{B}>0$. Hence, the for station $k\in[N]$, the arrival rate of bikes
is the product of the scaled popularity of station times the traveling bikes
$(M-\sum_{k}S_{k}^{N})$ weighted by the travel time, i.e.,
$\frac{p_{B}(k/N)}{\int_{0}^{1}p_{B}(\nu)d\nu}\bigl{(}M-\sum_{k}S_{k}^{N}\bigr{)}\mu_{B}$.
This implies the transition

$\displaystyle S\to S+\BFe_{k}^{N}$ at rate
$\displaystyle\frac{p_{B}(k/N)}{\int_{0}^{1}p_{B}(\nu)d\nu}(M-\sum_{k}S_{k}^{N})\mu_{B}\mathds{1}_{S_{k}<K_{B}}.$

For the above, the notation $\BFe_{k}^{N}$ refers to the unit vector of size
$N$ having a one at entry $k$ and zeros everywhere else.

####  5.2.2 Drift & Limiting System

To define the limiting system, we derive the drift implied by the above
transitions and the defined graphon $G_{B}$. We start by considering the
indicator state representation as outlined in Section 2, i.e., the drift is a
vector of size $K+1$ of $L_{2}(0,1]$ functions. In contrast to the stochastic
rates, for the drift, the sum over particles is replaced by an integral over
and the values of $G^{N}$ replaced by the graphon values. For a state
$\BFx=(x_{s})_{s=0,\dots,K}$ with $x_{s}\in L_{2}(0,1]$, the drift of the
system at $(u,s)\in(0,1]\times\mathcal{S}=\\{0,\dots,K_{B}\\}$ is defined by

$\displaystyle
F^{G}_{u,s}(\BFx):=x_{u,s}\left(\frac{p_{B}(u)}{\int_{0}^{1}p_{B}(\nu)d\nu}\bigl{(}M-\int_{0}^{1}\sum_{s=0}^{K_{B}}x_{v,s}dv\bigr{)}\mu_{B}\mathds{1}_{s<K_{B}}-\lambda_{B}\mathds{1}_{s>0}\right).$

By definition of the system and particularly the graphon, it is immediate that
the assumptions of Theorem 4.2 and Corollary 4.3 hold, therefore guaranteeing
the accuracy of the approximation.

####  5.2.3 Implementation & Results

For our simulations, we set the popularity function of the bike sharing system
to $p_{B}(v):=1-0.5v$. The travel rate and customer arrival rate are
$\mu_{B}=1$ and $\lambda_{B}=1$ respectively. For the plots of Figure 4 and
Figure 5 we calculate the sample mean by averaging over $7500$ simulations for
each system size. For each plot in the first Figure 4, the mean field
trajectory for single item and state is plotted against the sample mean. Here
the horizontal axis represents time $t$ and the vertical axis the probability
for the item to be in the state. In the second Figure 5, a comparison of the
sample mean against the values of the approximation for fixed time and state
is shown. As used for the theoretical results, the state of the stochastic
system is represented as a step function with constant value
$\mathbb{E}[X^{G^{N}}_{k,s}]$ on the intervals $(k-1/N,k/N]$ for $k=1..N$. In
accordance with our theoretical results, the plot shows that with increasing
system size the graphon mean field approximation becomes more accurate and
captures the state distribution of the particles well.

Simulation and Approximation for the Load Balancing Example

Figure 4: The plot shows the sample mean of $\mathbb{E}[X^{G^{N}}_{k,s}(t)]$
for $N=20,50$ and the value of the graphon mean field approximation
$\BFx^{G}_{k/N,s}(t)$ for $t\in[0,3]$. The plots are generated for the states
$s=3,6,10$ and $k=14,35$. Figure 5: For a fixed $t=2.5$ and state $s=3$, the
plots compare the graphon mean field approximation $\BFx^{G}$ against sample
means for systems of size $N=20,50$. As for the main theorem of our results,
we represent the values of the stochastic system as step function with
constant values on the intervals $(k-1/N,k/N],k\in[N]$. The figure shows that
for increasing system size, the values of the stochastic system indeed
approach the trajectory of the deterministic system. In the upper plots, the
approximation is plotted against values of the sample mean for one system size
$N=20,50$ each. In the lower plot, both functions are overlain for better
comparison.

##  6 Proofs

This section provides the proof of the main statements of this paper.

###  6.1 Proof of Lemma 4.1

###### Proof 6.1

Proof. We prove the uniqueness and continuous differentiability of the
differential equation (7) as well as the Lipschitz properties of the drift.
The proof is the consequence of two lemmas:

  * •

We show in Lemma 6.2 below that the drift is locally Lipschitz continuous in
$L_{2}$. By Driver (2003), this implies the existence of a local solution and
the uniqueness of this solution.

  * •

In Lemma 6.4 we show that the directional derivatives of $F^{G}$ are well-
defined and Lipschitz continuous. This property ensures that the $\BFx^{G}$ is
also continuously differentiable (for a proof of this property in general
Banach spaces see for example Driver (2003)).

######  Lemma 6.2 (Local Lipschitz Continuity of the Drift)

Let
$L_{\BFF^{G}}=2(C_{R}\left|\mathcal{S}\right|+2C_{R}\left|\mathcal{S}\right|^{2}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|)$
where $C_{R}$ is the bound on the rate functions,
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|$
the operator norm of the graphon $G$ and $\left|\mathcal{S}\right|$ the size
of the finite state space. Then for all
$\BFx_{.}=(x_{.,s})_{s\in\mathcal{S}},\BFy_{.}=(y_{.,s})_{s\in\mathcal{S}}$
with $x_{.,s},y_{.,s}\in L_{2}(0,1]$ and
$\left\|x\right\|,\left\|y\right\|\leq 1$ with
$\left|x_{s}(u)\right|,\left|y_{u,s}\right|\leq 1$ almost everywhere444Here,
almost everywhere means that the property holds except on a subset of $(0,1]$
with measure zero., we have:

$\displaystyle\left\|\BFF^{G}(\BFx)-\BFF^{G}(\BFy)\right\|_{L_{2}}\leq
L_{\BFF^{G}}\left\|\BFx-\BFy\right\|_{L^{2}}.$

###### Proof 6.3

Proof. Define
$\hat{\boldsymbol{F}}^{G}_{u,s}(\BFx,\BFy):=\BFR^{\text{uni}}_{u,s}\BFx_{u}+\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFy_{v}\
dv$, i.e., we replace $\BFx_{v}$ by $\BFy_{v}$ under the integral. Applying
the triangle inequality to the $L_{2}$ norm gives

$\displaystyle\left\|\BFF^{G}(\BFx)-\BFF^{G}(\BFy)\right\|_{L_{2}}$
$\displaystyle=\left\|\BFF^{G}(\BFx)-\hat{\boldsymbol{F}}^{G}(\BFx,\BFy)+\hat{\boldsymbol{F}}^{G}(\BFx,\BFy)-F(\BFy)\right\|_{L_{2}}$
$\displaystyle\leq\left\|\BFF^{G}(\BFx)-\hat{\boldsymbol{F}}^{G}(\BFx,\BFy)\right\|_{L_{2}}+\left\|\hat{\boldsymbol{F}}^{G}(\BFx,\BFy)-\BFF^{G}(\BFy)\right\|_{L_{2}}.$

Writing out the definitions and using the bounds of the rate vectors /
matrices and the bound of the graphon one immediately obtains

$\displaystyle\left\|\BFF^{G}(\BFx)-\hat{\boldsymbol{F}}^{G}(\BFx,\BFy)\right\|_{L_{2}}$
$\displaystyle=\sqrt{\sum_{s\in\mathcal{S}}\int_{0}^{1}\left(\BFR^{\text{uni}}_{u,s}\BFx_{u}+\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFx_{v}dv-\BFR^{\text{uni}}_{u,s}\BFx_{u}-\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFy_{v}\
dv\right)^{2}du}$
$\displaystyle=\sqrt{\sum_{s\in\mathcal{S}}\int_{0}^{1}\left(\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}(\BFx_{v}-\BFy_{v})dv\right)^{2}du}\leq
2\left|\mathcal{S}\right|^{3}C_{R}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|\left\|\BFx-\BFy\right\|_{L_{2}}$

and

$\displaystyle\left\|\hat{\boldsymbol{F}}^{G}(\BFx,\BFy)-\BFF^{G}(\BFy)\right\|_{L_{2}}$
$\displaystyle=\sqrt{\sum_{s\in\mathcal{S}}\int_{0}^{1}\left(\BFR^{\text{uni}}_{u,s}(\BFx_{u}-\BFy_{u})+(\BFx_{u}-\BFy_{u})^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFy_{v}\
dv\right)^{2}du}$
$\displaystyle\leq\bigl{(}2C_{R}\left|\mathcal{S}\right|^{2}+2\left|\mathcal{S}\right|^{3}C_{R}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|\bigr{)}\left\|\BFx-\BFy\right\|_{L_{2}}.$

Last, we conclude that
$\left\|F^{G}(\BFx)-F^{G}(\BFy)\right\|_{L_{2}}\leq(2(C_{R}\left|\mathcal{S}\right|^{2}+2C_{R}\left|\mathcal{S}\right|^{3}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|))\left\|\BFx-\BFy\right\|_{L_{2}}$,
where the $\left|\mathcal{S}\right|^{2}$ and $\left|\mathcal{S}\right|^{3}$
terms come from the vector and matrix notation used for the unilateral and
pairwise rates as well as the sum over $\mathcal{S}$.

######  Lemma 6.4 (Lipschitz Continuous Directional Derivative of $\BFF^{G}$)

The directional derivative of the drift $F^{G}$ for
$\BFx=(x_{s})_{s\in\mathcal{S}}$ with $x_{s}\in L_{2}(0,1]$ is given by

$\displaystyle
D_{x}F^{G}_{u,s}(\BFx)(\BFh):=\BFR^{\text{uni}}_{u,s^{\prime}}\BFh_{u}+\BFh_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFx_{v}dv+\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFh_{v}dv$
(14)

where $\BFR^{\text{uni}}_{u,s},\BFR^{\text{pair}}_{u,v,s}$ and $G$ are
Lebesgue integrable functions. Furthermore $D_{x}F^{G}_{u,s}(\BFx)(h)$ is
Lipschitz-continuous in $x$.

###### Proof 6.5

Proof. Recall the definition of
$\BFF^{G}(\BFx)=\BFR^{\text{uni}}_{u,s^{\prime}}\BFx_{u}+\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s^{\prime}}G_{u,v}\BFx_{v}\
dv$. We define $D_{x}F^{G}_{u,s}(\BFx)(\BFh)$ for
$(u,s)\in(0,1]\times\mathcal{S}$ as in Equation (14). To see that
$D_{x}\BFF^{G}(x)(h)$ is a directional derivative of $\BFF^{G}$ in $\BFx\in
L_{2}(0,1]^{\mathcal{S}}$ in direction $h\in L_{2}$, we show that it fulfills

$\displaystyle\lim_{\epsilon\downarrow
0}\frac{1}{\epsilon}\left\|\BFF^{G}(\BFx+\epsilon\BFh)-\BFF^{G}(\BFx)-D_{x}\BFF^{G}(\BFx)(\epsilon\BFh)\right\|_{L_{2}}=0.$

By definition

$\displaystyle\left\|\BFF^{G}(\BFx+\epsilon\BFh)-\BFF^{G}(x)-D_{x}\BFF^{G}(\BFx)(\BFh)\right\|^{2}_{L_{2}}$
$\displaystyle=\int_{0}^{1}\sum_{s}\left(\BFR^{\text{uni}}_{u,s}(\BFx_{u}+\epsilon\BFh_{u})+(\BFx_{u}+\epsilon\BFh_{u})^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}(\BFx_{v}+\epsilon\BFh_{v})\
dv\right.$
$\displaystyle\qquad-\BFR^{\text{uni}}_{u,s^{\prime}}\BFx_{u}-\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s^{\prime}}G_{u,v}\BFx_{v}\
dv$
$\displaystyle\qquad\left.-\epsilon\BFR^{\text{uni}}\BFh_{u}-\epsilon\BFh_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFx_{v}dv-\epsilon\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFh_{v}dv\right)^{2}du$
$\displaystyle=\int_{0}^{1}\sum_{s}\left(\epsilon\BFR^{\text{uni}}_{u,s}\BFh_{u}+\epsilon\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFh_{v}dv+\epsilon\BFh_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFx_{v}\
dv\right.$
$\displaystyle\qquad+\epsilon^{2}\BFh_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFh_{v}\
dv$
$\displaystyle\qquad\left.-\epsilon(\BFR^{\text{uni}}\BFh_{u}-\BFh_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFx_{v}dv-\BFx_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFh_{v}dv)\right)^{2}du$
$\displaystyle=\int_{0}^{1}\sum_{s}\left(\epsilon^{2}\BFh_{u}^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s}G_{u,v}\BFh_{v}dv\right)^{2}du.$

Taking dividing by $\epsilon$ and taking $\epsilon\to 0$ shows that
$D_{x}\BFF^{G}(\BFx)(\BFh)$ is indeed a directional derivative of $\BFF^{G}$
in $\BFx$. By definition, $D_{x}\BFF^{G}_{u,s}(\BFx)(\BFh)$ is linear in
$\BFx$ and as $\BFR^{\text{uni}}_{u,s},\BFR^{\text{pair}}_{u,v,s}$ and
$G_{u,v}$ are bounded functions, the Lipschitz continuity in $L_{2}$ follows.

###  6.2 Proof of Theorem 4.2

In the following, we present the proof of Theorem 4.2.

###### Proof 6.6

Proof. Let $\BFX^{G^{N}}_{t}$ be the $L_{2}(0,1]^{\left|\mathcal{S}\right|}$
representation of the stochastic system as defined in 2.4. For a fixed $t\geq
0$, define

$\displaystyle\boldsymbol{\nu}^{N}(\tau)=\mathbb{E}[\BFx^{G}\bigl{(}t-\tau,\BFX^{G^{N}}(\tau)\bigr{)}]$
(15)

for which we suppress the dependence of the expectation on the graph $G^{N}$
and initial state $\BFX^{G^{N}}(0)$ from here on. Using the definition of
$\boldsymbol{\nu}^{N}$ we rewrite

$\displaystyle\left\|\mathbb{E}[\BFX^{G^{N}}(t)]-x^{G}(t,\BFX^{G^{N}}(0))\right\|_{L_{2}}$
$\displaystyle=\sqrt{\sum_{s\in\mathcal{S}}\int_{0}^{1}\left(\mathbb{E}[X^{G^{N}}_{u,s}(t)]-x^{G}_{u,s}(t,\BFX^{G^{N}}(0))\right)^{2}du}$
$\displaystyle=\sqrt{\sum_{s\in\mathcal{S}}\int_{0}^{1}\left(\nu_{u,s}^{N}(t)-\nu_{u,s}^{N}(0)\right)^{2}du}.$
(16)

To compare the drift of the ODE against the drift of the stochastic process
and ultimately bound them, we first want to rewrite
$\nu_{u,s}^{N}(t)-\nu_{u,s}^{N}(0)=\int_{0}^{t}\frac{d}{d\tau}\nu^{N}_{u,s}(\tau)$,
with $\frac{d}{d\tau}\nu^{N}(\tau)$ being the quantity that fulfills
$\boldsymbol{\nu}^{N}(b)-\boldsymbol{\nu}^{N}(a)=\int_{a}^{b}\frac{d}{d\tau}\boldsymbol{\nu}^{N}(\theta)d\theta$
for arbitrary $a,b\in[0,t]$ with $a<b$. In Lemma A.1 we show that
$\frac{d}{d\tau}\boldsymbol{\nu}(\theta)$ exists almost everywhere for
$\theta\in(0,t)$ and is almost everywhere equal to

$\displaystyle\mathbb{E}\left[\sum_{s^{\prime}}\int_{0}^{1}\left[D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\left(\BFF^{G}(\BFX^{G^{N}}(\tau))-\BFF^{G^{N}}(\BFX^{G^{N}}(\tau))\right)\right]_{v,s^{\prime}}dv\right]$
(17)
$\displaystyle\qquad+\mathbb{E}\left[\tilde{R}_{1}(\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)))\right]d\tau.$
(18)

In the above sum,
$D_{x}\BFx^{G}(\tau,\BFX^{G^{N}}(\tau))\left(\BFF^{G}(\BFX^{G^{N}}(\tau))-\BFF^{G^{N}}(\BFX^{G^{N}}(\tau))\right)$
is the directional derivative of $\BFx^{G}$ in its initial condition in
direction $\BFF^{G}(\BFX^{G^{N}}(\tau))-\BFF^{G^{N}}(\BFX^{G^{N}}(\tau))$. We
aim to bound the sum by using the properties of the derivative of $\BFx^{G}$
and the differences between the drift of the deterministic and stochastic
system with respect to the $L_{2}$ norm. The technical details to bound the
remainder term and the difference between the drifts are moved to the Lemmas
A.3 and A.5. From the application of the latter lemmas, one obtains the bounds

$\displaystyle
C_{D_{x(0)}}\left(\frac{2L_{\BFR^{\text{pair}}}+16C^{2}_{\BFR^{\text{pair}}}K_{\BFR^{\text{pair}}}}{N}-\frac{16C^{2}_{\BFR^{\text{pair}}}K^{2}_{\BFR^{\text{pair}}}}{N^{2}}+C_{\BFR^{\text{pair}}}^{2}\left|\mathcal{S}\right|^{2}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G-G^{N}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|\right)$

for Equation (17) and $C_{\tilde{R}}/N$ for Equation (18). It follows by
applying the bounds and rearranging terms, that

$\displaystyle\left\|\mathbb{E}[\BFX^{G^{N}}(t)]-x^{G}(t,\BFX^{G^{N}}(0))\right\|_{L_{2}}\leq\frac{C_{A}}{N}+C_{B}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|$

for some finite constants $C_{A},C_{B}>0$ additionally to the previous bound
also depend on $t$. This concludes the proof.

###  6.3 Proof of Corollary 4.4 and Corollary 4.3

In this section, we prove the corollaries of Theorem 4.2. In each case, we use
an additional lemma to obtain a bound on
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|$.

###### Proof 6.7

Proof.[Proof of Corollary 4.4] By the application of Theorem 4.2 we have

$\displaystyle\left\|\mathbb{E}[\BFX^{G^{N}}(t)]-x^{G}(t,\BFX^{G^{N}}(0))\right\|_{L_{2}}$
$\displaystyle=\frac{C_{A}}{N}+C_{B}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|.$

It therefore remains to bound the $L_{2}$ distance
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|$
between the graphon $G$ and graph $G^{N}$. By application of Lemma 3.2 with
probability at least $1-\delta$ and for ‘large enough’ $N$

$\displaystyle\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|\leq\sqrt{\frac{4\log{(2N/\delta)}}{N}}+2\sqrt{\frac{(L_{G}^{2}-K_{G}^{2})}{N^{2}}+\frac{K_{G}}{N}}=:\psi_{\delta,G}(N).$

Defining and substitution $\delta=2/N$ into the above equation concludes the
proof.

###### Proof 6.8

Proof.[Proof of Corollary 4.3] The corollary is a direct implication of
Theorem 4.4 and Lemma A.7: due to $G^{N}$ representing a discretized version
of $G$ it directly follows that
$\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G-G^{N}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|$
is of order $\mathcal{O}(1/N)$, which concludes the proof.

##  7 Conclusion

In this paper, we study an approximation for a system of particles interacting
on a graph. We show that when the interacting graph converges to a graphon,
the underlying behavior of the stochastic system converges to a deterministic
limit, which we call the graphon mean field approximation. While this result
is similar to other results in the literature – that show that graphon mean
field approximations are asymptotically exact in some settings –, our main
contribution is to provide precise bounds on the accuracy of this
approximation. We showed indeed that the distance between the original
finite-$N$ system and the graphon mean field approximation can be bounded by a
term $\mathcal{O}(1/N)$ that depends on the number of particles of the system
plus a term
$\mathcal{O}(\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G^{N}-G\right|\kern-1.07639pt\right|\kern-1.07639pt\right|)$
that depends on the distance between the original graph $G^{N}$ and the
graphon $G$, when measured as a $L_{2}$ operator.

This paper aims to be methodological. It shows that it is possible to obtain
bounds for a system with a graph structure, and not mere asymptotic
convergence results. To keep the presentation reasonable, we intentionally
considered a relatively simple model, for instance by restricting our
attention to pairwise interactions between nodes or unilateral jumps, and by
considering that the rate functions $r^{N}$ are discretized versions of the
limiting rates $r$. We believe that by doing so, the proof is easier to follow
and could then be adapted to more general cases.

The focus of this paper is to study the case of dense graphs that converge to
graphons. This implies that the total number of edges per node is of order
$\mathcal{O}(N)$. We believe that our methodology could be applicable in the
not-so-dense case when the number of edges per node goes to infinity at a sub-
linear rate. This would probably give bounds that converge more slowly to
zero. Complementary to this paper, another interesting question is the case of
sparse graphs, where the number of neighbors per node remains bounded when the
number of nodes $N$ goes to infinity. Yet, this would require a fundamentally
different approach: studying such a problem is out of the scope of our tools
since the graphon mean field approximation is not asymptotically exact for
sparse graphs.

## References

  * Abbe (2018) Abbe E (2018) Community Detection and Stochastic Block Models. _Foundations and Trends® in Communications and Information Theory_ 14(1-2):1–162, ISSN 1567-2190, 1567-2328, URL http://dx.doi.org/10.1561/0100000067.
  * Allmeier and Gast (2022) Allmeier S, Gast N (2022) Mean Field and Refined Mean Field Approximations for Heterogeneous Systems: It Works! _Proceedings of the ACM on Measurement and Analysis of Computing Systems_ 6(1):13:1–13:43, URL http://dx.doi.org/10.1145/3508033.
  * Aurell et al. (2022) Aurell A, Carmona R, Dayanıklı G, Laurière M (2022) Finite State Graphon Games with Applications to Epidemics. _Dynamic Games and Applications_ 12(1):49–81, ISSN 2153-0793, URL http://dx.doi.org/10.1007/s13235-021-00410-2.
  * Avella-Medina et al. (2018) Avella-Medina M, Parise F, Schaub MT, Segarra S (2018) Centrality measures for graphons: Accounting for uncertainty in networks. URL http://dx.doi.org/10.1109/TNSE.2018.2884235.
  * Bayraktar et al. (2023) Bayraktar E, Chakraborty S, Wu R (2023) Graphon mean field systems. _The Annals of Applied Probability_ 33(5):3587–3619, ISSN 1050-5164, 2168-8737, URL http://dx.doi.org/10.1214/22-AAP1901.
  * Bayraktar and Wu (2021) Bayraktar E, Wu R (2021) Mean field interaction on random graphs with dynamically changing multi-color edges. _Stochastic Processes and their Applications_ 141:197–244, ISSN 0304-4149, URL http://dx.doi.org/10.1016/j.spa.2021.07.005.
  * Bet et al. (2024) Bet G, Coppini F, Nardi FR (2024) Weakly interacting oscillators on dense random graphs. _Journal of Applied Probability_ 61(1):255–278, ISSN 0021-9002, 1475-6072, URL http://dx.doi.org/10.1017/jpr.2023.34.
  * Bhamidi et al. (2019) Bhamidi S, Budhiraja A, Wu R (2019) Weakly interacting particle systems on inhomogeneous random graphs. _Stochastic Processes and their Applications_ 129(6):2174–2206, ISSN 0304-4149, URL http://dx.doi.org/10.1016/j.spa.2018.06.014.
  * Blanchard and Brüning (2015) Blanchard P, Brüning E (2015) _Mathematical Methods in Physics: Distributions, Hilbert Space Operators, Variational Methods, and Applications in Quantum Physics_ , volume 69 of _Progress in Mathematical Physics_ (Cham: Springer International Publishing), ISBN 978-3-319-14044-5, URL http://dx.doi.org/10.1007/978-3-319-14045-2.
  * Braverman and Dai (2017) Braverman A, Dai JG (2017) Stein’s method for steady-state diffusion approximations of M/Ph/n+M systems. _The Annals of Applied Probability_ 27(1), ISSN 1050-5164, URL http://dx.doi.org/10.1214/16-AAP1211.
  * Braverman et al. (2017) Braverman A, Dai JG, Feng J (2017) Stein’s method for steady-state diffusion approximations: An introduction through the Erlang-A and Erlang-C models. _Stochastic Systems_ 6(2):301–366.
  * Budhiraja et al. (2019) Budhiraja A, Mukherjee D, Wu R (2019) Supermarket model on graphs. _The Annals of Applied Probability_ 29(3), ISSN 1050-5164, URL http://dx.doi.org/10.1214/18-AAP1437.
  * Caines and Huang (2021) Caines PE, Huang M (2021) Graphon Mean Field Games and Their Equations. _SIAM Journal on Control and Optimization_ 59(6):4373–4399, ISSN 0363-0129, 1095-7138, URL http://dx.doi.org/10.1137/20M136373X.
  * Decreusefond et al. (2012) Decreusefond L, Dhersin JS, Moyal P, Tran VC (2012) Large graph limit for an SIR process in random network with heterogeneous connectivity. _The Annals of Applied Probability_ 22(2):541–575, ISSN 1050-5164, 2168-8737, URL http://dx.doi.org/10.1214/11-AAP773.
  * Delmas et al. (2023) Delmas JF, Frasca P, Garin F, Tran VC, Velleret A, Zitt PA (2023) Individual based SIS models on (not so) dense large random networks.
  * Driver (2003) Driver BK (2003) Analysis tools with applications. _Lecture notes_ .
  * Fricker and Gast (2016) Fricker C, Gast N (2016) Incentives and redistribution in homogeneous bike-sharing systems with stations of finite capacity. _EURO Journal on Transportation and Logistics_ 5(3):261–291, ISSN 2192-4384, URL http://dx.doi.org/10.1007/s13676-014-0053-5.
  * Fricker et al. (2012) Fricker C, Gast N, Mohamed H (2012) Mean field analysis for inhomogeneous bike sharing systems. _AofA_ , volume DMTCS Proceedings vol. AQ, 23rd Intern. Meeting on Probabilistic, Combinatorial, and Asymptotic Methods for the Analysis of Algorithms (AofA’12) (Montreal, Canada: DMTCS), URL http://dx.doi.org/10.46298/dmtcs.3006.
  * Ganguly (2022) Ganguly A (2022) _Non-Markovian Interacting Particle Systems on Large Sparse Graphs: Hydrodynamic Limits and Marginal Characterizations_. Ph.D. thesis, Brown Univerity.
  * Ganguly and Ramanan (2022) Ganguly A, Ramanan K (2022) Hydrodynamic Limits of non-Markovian Interacting Particle Systems on Sparse Graphs.
  * Gast (2017) Gast N (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):1–26, ISSN 24761249, URL http://dx.doi.org/10.1145/3084454.
  * Gast and Van Houdt (2017) Gast N, Van Houdt B (2017) A Refined Mean Field Approximation. _Proceedings of the ACM on Measurement and Analysis of Computing Systems_ 1(2):33:1–33:28, URL http://dx.doi.org/10.1145/3154491.
  * Keliger (2023) Keliger D (2023) Universality of Sis Epidemics Starting from Small Initial Conditions. URL http://dx.doi.org/10.2139/ssrn.4543174.
  * Kurtz (1970) Kurtz TG (1970) Solutions of ordinary differential equations as limits of pure jump markov processes. _Journal of Applied Probability_ 7(1):49–58, ISSN 0021-9002, 1475-6072, URL http://dx.doi.org/10.2307/3212147.
  * Kurtz (1971) Kurtz TG (1971) Limit theorems for sequences of jump Markov processes approximating ordinary differential processes. _Journal of Applied Probability_ 8(2):344–356, ISSN 0021-9002, 1475-6072, URL http://dx.doi.org/10.2307/3211904.
  * Le Boudec et al. (2007) Le Boudec JY, McDonald D, Mundinger J (2007) A Generic Mean Field Convergence Result for Systems of Interacting Objects. _Fourth International Conference on the Quantitative Evaluation of Systems (QEST 2007)_ , 3–18 (Edinburgh, Scotland, UK: IEEE), ISBN 978-0-7695-2883-0, URL http://dx.doi.org/10.1109/QEST.2007.8.
  * Lovász (2012) Lovász L (2012) _Large Networks and Graph Limits_ , volume 60 of _Colloquium Publications_ (Providence, Rhode Island: American Mathematical Society), ISBN 978-0-8218-9085-1, URL http://dx.doi.org/10.1090/coll/060.
  * McKean (1967) McKean HP (1967) Propagation of chaos for a class of non-linear parabolic equations. _Stochastic Differential Equations (Lecture Series in Differential Equations, Session 7, Catholic Univ., 1967)_ 41–57.
  * Mitzenmacher (2001) Mitzenmacher M (2001) The power of two choices in randomized load balancing. _IEEE Transactions on Parallel and Distributed Systems_ 12(10):1094–1104, ISSN 10459219, URL http://dx.doi.org/10.1109/71.963420.
  * Norman (1972) Norman MF (1972) _Markov Processes and Learning Models_ , volume 84 (Academic Press New York), ISBN 0-12-521450-2.
  * Ramanan (2022) Ramanan K (2022) Beyond mean-field limits for the analysis of large-scale networks. _Queueing Systems_ 100(3-4):345–347, ISSN 0257-0130, 1572-9443, URL http://dx.doi.org/10.1007/s11134-022-09845-9.
  * Roy et al. (2023) Roy A, Singh C, Narahari Y (2023) Recent advances in modeling and control of epidemics using a mean field approach. _Sādhanā_ 48(4):207, ISSN 0973-7677, URL http://dx.doi.org/10.1007/s12046-023-02268-z.
  * Rutten and Mukherjee (2023) Rutten D, Mukherjee D (2023) Mean-field Analysis for Load Balancing on Spatial Graphs. _Abstract Proceedings of the 2023 ACM SIGMETRICS International Conference on Measurement and Modeling of Computer Systems_ , 27–28, SIGMETRICS ’23 (New York, NY, USA: Association for Computing Machinery), ISBN 9798400700743, URL http://dx.doi.org/10.1145/3578338.3593552.
  * Stein (1986) Stein C (1986) Approximate computation of expectations. _Lecture Notes-Monograph Series_ 7:i–164.
  * Van Der Hofstad (2017) Van Der Hofstad R (2017) _Random Graphs and Complex Networks_ (Cambridge: Cambridge University Press), ISBN 978-1-316-77942-2, URL http://dx.doi.org/10.1017/9781316779422.
  * Vizuete et al. (2020) Vizuete R, Frasca P, Garin F (2020) Graphon-Based Sensitivity Analysis of SIS Epidemics. _IEEE Control Systems Letters_ 4(3):542–547, ISSN 2475-1456, URL http://dx.doi.org/10.1109/LCSYS.2020.2971021.
  * Zhao and Mukherjee (2024) Zhao Z, Mukherjee D (2024) Optimal Rate-Matrix Pruning For Heterogeneous Systems. _ACM SIGMETRICS Performance Evaluation Review_ 51(4):26–27, ISSN 0163-5999, URL http://dx.doi.org/10.1145/3649477.3649492.
  * Zhao et al. (2024) Zhao Z, Mukherjee D, Wu R (2024) Exploiting Data Locality to Improve Performance of Heterogeneous Server Clusters. _Stochastic Systems_ ISSN 1946-5238, URL http://dx.doi.org/10.1287/stsy.2022.0040.

##  Appendix A Lemmas

This section contains the most technical lemmas of the paper.

###  A.1 Derivative of $\nu$ (Lemma A.1)

The following Lemma derives a representation of the derivative of $\nu$, as
defined in the proof of Theorem 4.2, which depends on the derivative of the
graphon mean field approximation with respect to the initial condition, the
transitions of the stochastic system and the deterministic drift.

###### Lemma A.1

For $\nu^{N}(\tau)$ and $\tau\in(0,t)$ we define

$\displaystyle\frac{d}{d\tau}\nu^{N}(\tau):=\lim_{h\downarrow
0}\frac{1}{h}\int_{\tau}^{\tau+h}\lim_{r\downarrow
0}\frac{1}{r}\left(\mathbb{E}[\BFx^{G}(t-(\theta+r),\BFX^{G^{N}}(\theta+r))]-\mathbb{E}[\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta))]\right)d\theta.$

for which, by construction and continuity properties of $\BFx^{G}$ and
$\BFX^{G^{N}}$, it holds that for arbitrary $a,b\in[0,T],a<b$,
$\int_{a}^{b}\frac{d}{d\tau}\nu^{N}(\tau)=\nu^{N}(b)-\nu^{N}(a)$. Furthermore,
the right-hand side is almost everywhere equal to

$\displaystyle\mathbb{E}[\sum_{s^{\prime}}\int_{0}^{1}\left[D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\left(\BFF^{G}(\BFX^{G^{N}}(\tau))-\BFF^{G^{N}}(\BFX^{G^{N}}(\tau))\right)\right]_{v,s}du]$
$\displaystyle\quad+\mathbb{E}[\tilde{R}_{1}(\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)))].$

which is well defined and bounded due to the properties of the drifts
$\BFF^{G^{N}},\BFF^{G}$ and the approximation $\BFx^{G}$. In this context,
$\tilde{R}_{1}$ refers to the summed Taylor remainder term defined by

$\displaystyle\tilde{R}_{1}\bigl{(}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\bigr{)}:=$
$\displaystyle\sum_{k,s_{k},s_{k}^{\prime}}R_{1}(\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)),\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N})$
$\displaystyle\qquad\times\left[X^{G^{N}}_{k,s_{k}}r^{N,\text{uni}}_{k,s_{k}\to
s_{k}^{\prime}}+X_{k,s_{k}}\sum_{l\in[N]}\sum_{s_{l}}r^{N,\text{pair}}_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}\frac{G^{N}_{kl}}{N}X_{l,s_{l}}\right].$

###### Proof A.2

Proof. Recall, that $\BFX^{G^{N}}$ is in the following handled as a $L_{2}$
function as described in Section 2.4. We start the proof by rewriting

$\displaystyle\frac{d}{d\tau}\nu(\tau)$ $\displaystyle:=\lim_{h\downarrow
0}\frac{1}{h}\int_{\tau}^{\tau+h}\lim_{r\downarrow
0}\frac{1}{r}\left(\mathbb{E}[\BFx^{G}(t-(\theta+r),\BFX^{G^{N}}(\theta+r))]-\mathbb{E}[\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta))]\right)d\theta$
$\displaystyle=\lim_{h\downarrow
0}\frac{1}{h}\int_{\tau}^{\tau+h}\lim_{r\downarrow
0}\frac{1}{r}\left(\mathbb{E}[\BFx^{G}(t-(\theta+r),\BFX^{G^{N}}(\theta+r))]-\mathbb{E}[\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta+r))]\right)d\theta$
$\displaystyle\quad+\lim_{h\downarrow
0}\frac{1}{h}\int_{\tau}^{\tau+h}\lim_{r\downarrow
0}\frac{1}{r}\left(\mathbb{E}[\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta+r))]-\mathbb{E}[\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta))]\right)d\theta$
$\displaystyle=\lim_{h\downarrow
0}\frac{1}{h}\int_{\tau}^{\tau+h}\lim_{r\downarrow
0}\frac{1}{r}\left(\mathbb{E}\left[\BFx^{G}(t-\theta,\BFx^{G}(-r,\BFX^{G^{N}}(\theta+r)))-\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta+r))]\right]\right)d\theta$
$\displaystyle\quad+\lim_{h\downarrow
0}\frac{1}{h}\int_{\tau}^{\tau+h}\lim_{r\downarrow
0}\frac{1}{r}\left(\mathbb{E}\left[\mathbb{E}[\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta+r))-\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta))\mid\BFX^{G^{N}}(\tau)]\right]\right)d\theta$
$\displaystyle=\lim_{h\downarrow
0}\frac{1}{h}\int_{\tau}^{\tau+h}\mathbb{E}\left[\lim_{r\downarrow
0}\frac{1}{r}\left(\BFx^{G}(t-\theta,\BFx^{G}(-r,\BFX^{G^{N}}(\theta+r)))-\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta+r))\right)\right]d\theta$
$\displaystyle\quad+\lim_{h\downarrow
0}\frac{1}{h}\int_{\tau}^{\tau+h}\mathbb{E}\left[\lim_{r\downarrow
0}\frac{1}{r}\left(\mathbb{E}[\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta+r))-\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta))\mid\BFX^{G^{N}}(\tau)]\right)\right]d\theta$

For the first equality, an artificial zero is added by adding and subtracting
the term $\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta+r))$ and using dominated
convergence in combination with the boundedness of the solution of the
differential equation and the stochastic system almost everywhere. For the
second equality, we use the tower property to write the difference in the
conditional expectation with respect to $\BFX^{G^{N}}(\tau)$. The last
equality, follows by dominated convergence which allows to take the limit
inside the outer expectation and by definition of $\BFx^{G}$ as well as the
drift $\BFF^{G}$ which allows to rewrite
$\BFx^{G}(t-(\theta+r),\BFX^{G^{N}}(\theta+r))=\BFx^{G}(t-\theta,\BFx^{G}(-r,\BFX^{G^{N}}(\theta+r)))$.
By definition, $\BFx^{G}$ is continuously differentiable in $L_{2}$ with
respect to its initial condition. The property follows by a general version of
Grönwalls inequality and the Lipschitz properties of the drift, also see
Driver (2003) for a proof of this property. Using the differentiability in
combination with the Lebesgue differentiation theorem for $L_{p}$ spaces and
application of the chain rule we can rewrite

$\displaystyle\lim_{h\downarrow
0}\frac{1}{h}\int_{\tau}^{\tau+h}\mathbb{E}\left[\lim_{r\downarrow
0}\frac{1}{r}\left(\BFx^{G}(t-\theta,\BFx^{G}(-r,\BFX^{G^{N}}(\theta+r)))-\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta+r))\right)\right]d\theta$
$\displaystyle=\mathbb{E}\left[D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\BFF^{G}(\BFX^{G^{N}}(\tau))\right]$
(19)

where equality holds almost everywhere and
$D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\BFF^{G}(\BFX^{G^{N}}(\tau))$ is the
directional derivative of $\BFx^{G}$ in its initial condition with direction
$\BFF^{G}(\BFX^{G^{N}}(\tau))$. For the second integral, we apply Lebesgue
differentiation theorem in combination with the definition of the transitions
for the stochastic system. In particular, we use that the probability for a
small $d\tau>0$ the state of the stochastic system changes to
$\BFX^{G^{N}}(\tau+d\tau)=\BFX^{G^{N}}(\tau)+\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N}$
is given by $d\tau X^{G^{N}}_{k,s_{k}}(\tau)r^{N,\text{uni}}_{k,s_{k}\to
s_{k}^{\prime}}+d\tau
X^{G^{N}}_{k,s_{k}}(\tau)\sum_{l\in[N]}\sum_{s_{l}}r^{N,\text{pair}}_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}\frac{G^{N}_{kl}}{N}X^{G^{N}}_{l,s_{l}}(\tau)+d\tau)$.
This shows

$\displaystyle\lim_{h\downarrow
0}\int_{\tau}^{\tau+h}\mathbb{E}\left[\lim_{r\downarrow
0}\frac{1}{r}\left(\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta+r))-\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta))\right)\right]d\theta$
$\displaystyle=\mathbb{E}\left[\lim_{d\tau\downarrow
0}\frac{1}{d\tau}\sum_{k,s_{k},s_{k}^{\prime}}\left(\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)+\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N})-\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\right)\right.$
$\displaystyle\qquad\times\left.\left(d\tau
X^{G^{N}}_{k,s_{k}}(\tau)r^{N,\text{uni}}_{k,s_{k}\to s_{k}^{\prime}}+d\tau
X^{G^{N}}_{k,s_{k}}(\tau)\sum_{l\in[N]}\sum_{s_{l}}r^{N,\text{pair}}_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}\frac{G^{N}_{kl}}{N}X^{G^{N}}_{l,s_{l}}(\tau)+o(d\tau)\right)\right]$

where equality holds almost everywhere. Lastly, by continuous
differentiability of $\BFx^{G}$ with respect to its initial condition, we
apply the $L_{2}$ version of the Taylor expansion, see Lemma A.9, for
$\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)+\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N})$
around $\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))$. Note again, that the equality
for the expansion hold almost everywhere. This allows to rewrite

$\displaystyle\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)+\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N})=$
$\displaystyle\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))+\int_{0}^{1}\sum_{s^{\prime}}\left[D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\left(\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N}\right)\right]_{v,s^{\prime}}$
$\displaystyle+R_{1}(\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)),\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N}).$

Here,
$D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\left(\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N}\right)$
is the directional derivative of $\BFx^{G}$ in direction
$\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N}$ which should be interpreted
as $L_{2}$ similar as for the $L_{2}$ representation of $\BFX^{G^{N}}$. By
$R_{1}(\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)),\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N})$
we denote the Taylor remainder term of first order defined as

$\displaystyle
R_{1}(\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)),\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N})$
$\displaystyle=\int_{0}^{1}(1-\delta)\int_{0}^{1}\sum_{s^{\prime}}\left[\left(D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)+\delta(e_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N}))-D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\right)\right.$
$\displaystyle\qquad\left.\times(e_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N})\right]_{v,s^{\prime}}\
dv\ d\delta.$

Adding this into the previous right-hand side and rearranging the terms yields
by using the linearity of the directional derivative

$\displaystyle\mathbb{E}\left[\sum_{k,s_{k},s_{k}^{\prime}}\left(\int_{0}^{1}\sum_{s^{\prime}}\left[D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))(\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N})\right]_{v,s^{\prime}}dv+R_{1}(\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)),\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N})\right)\right.$
$\displaystyle\qquad\times\left.\left(X^{G^{N}}_{k,s_{k}}(\tau)r^{N,\text{uni}}_{k,s_{k}\to
s_{k}^{\prime}}+X^{G^{N}}_{k,s_{k}}(\tau)\sum_{l\in[N]}\sum_{s_{l}}r^{N,\text{pair}}_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}\frac{G^{N}_{kl}}{N}X^{G^{N}}_{l,s_{l}}(\tau)\right)\right]$
$\displaystyle=\mathbb{E}\left[\int_{0}^{1}\sum_{s^{\prime}}\left[D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\BFF^{G^{N}}(\BFX^{G^{N}}(\tau))\right]_{v,s^{\prime}}dv\right]+\mathbb{E}\left[\tilde{R}_{1}\bigl{(}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\bigr{)}\right].$

with

$\displaystyle\tilde{R}_{1}\bigl{(}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\bigr{)}:=$
$\displaystyle\sum_{k,s_{k},s_{k}^{\prime}}R_{1}(\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)),\BFe_{k,s_{k}^{\prime}}^{N}-\BFe_{k,s_{k}}^{N})$
(20) $\displaystyle\times\left(X^{G^{N}}_{k,s_{k}}r^{N,\text{uni}}_{k,s_{k}\to
s_{k}^{\prime}}+X^{G^{N}}_{k,s_{k}}\sum_{l\in[N]}\sum_{s_{l}}r^{N,\text{pair}}_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}\frac{G^{N}_{kl}}{N}X^{G^{N}}_{l,s_{l}}\right).$

Ultimately, we rewrite the initial integral by using the formulations of
Equations (19) and (20) to

$\displaystyle\lim_{h\downarrow
0}\frac{1}{h}\int_{\tau}^{\tau+h}\lim_{r\downarrow
0}\frac{1}{r}\left(\mathbb{E}[\BFx^{G}(t-(\theta+r),\BFX^{G^{N}}(\theta+r))]-\mathbb{E}[\BFx^{G}(t-\theta,\BFX^{G^{N}}(\theta))]\right)d\theta$
$\displaystyle=\mathbb{E}[\sum_{s^{\prime}}\int_{0}^{1}\left[D_{x}\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau))\left(F^{G}(\BFX^{G^{N}}(\tau))-\BFF^{G^{N}}(\BFX^{G^{N}}(\tau))\right)\right]_{v,s^{\prime}}dv\mid
G^{N},\BFX^{G^{N}}(0)]$
$\displaystyle\qquad+\mathbb{E}[\tilde{R}_{1}(\BFx^{G}(t-\tau,\BFX^{G^{N}}(\tau)))\mid
G^{N},\BFX^{G^{N}}(0)]$

where equality holds almost everywhere. This concludes the proof.

###  A.2 Bound on the Taylor Remainder Term (Lemma A.3)

The subsequent Lemma states an upper bound for the remainder term arising in
Lemma A.1 which ultimately is used to obtain the bound for Theorem 4.2.

###### Lemma A.3

For any attainable state $\BFX\in\mathcal{X}^{N}$ of the stochastic process
and $\tau\geq 0$, the expectation of the Taylor remainder term
$\tilde{R}_{1}\bigl{(}x^{G}_{u,s}(t,\BFX)\bigr{)}$555In order for
$\BFx^{G}_{t}$ to be properly defined, $\BFX$ should be understood in this
context as a $L_{2}(0,1]^{\left|\mathcal{S}\right|}$ function. as defined in
Equation (20) is of order $\mathcal{O}(1/N)$.

###### Proof A.4

Proof. For a state $\BFX\in\mathcal{X}^{N}$, recall the definition of Equation
(20)

$\displaystyle\tilde{R}_{1}\bigl{(}x^{G}_{u,s}(\tau,\BFX)\bigr{)}=\sum_{k\in[N],s_{k},s_{k}^{\prime}\in\mathcal{S}}$
$\displaystyle
R_{1}(x^{G}_{u,s}(\tau,.),\BFX,e^{(k,s_{k}^{\prime})}-e^{(k,s_{k})})$
$\displaystyle\quad\times\bigl{(}X^{G^{N}}_{k/N,s_{k}}r^{N,\text{uni}}_{k,s_{k}\to
s_{k}^{\prime}}+X^{G^{N}}_{k/N,s_{k}}\sum_{s_{l}\in\mathcal{S}}\int_{0}^{1}r^{\text{pair}}_{k/N,v,s_{k}\to
s_{k}^{\prime},s_{l}}G^{N}_{k/N,v}X^{G^{N}}_{v,s_{l}}dv\bigr{)}$

where we use that $r^{N,\text{pair}}_{k,l,s_{k}\to
s_{k}^{\prime},s_{l}}=r^{\text{pair}}_{k/N,v,s_{k}\to s_{k}^{\prime},s_{l}}$
for $v\in(l-1/N,l]$ and $R_{1}$ given by

$\displaystyle
R_{1}(x^{G}_{u,s}(\tau,.),\BFX,\BFe^{N}_{(k,s_{k})}-\BFe^{N}_{(k,s_{k})})=$
$\displaystyle\int_{0}^{1}(1-\delta)\sum_{s^{\prime}}\int_{0}^{1}\left[\left(D_{x}\BFx^{G}\bigl{(}\tau-\delta,\BFX+\delta(\BFe^{N}_{(k,s_{k})}-\BFe^{N}_{(k,s_{k})})\bigr{)}-D_{x}\BFx^{G}\bigl{(}\tau,\BFX\bigr{)}\right)(\BFe^{N}_{(k,s_{k}^{\prime})}-\BFe^{N}_{(k,s_{k})})\right]_{v,s^{\prime}}dv\
d\delta.$

To be clear, as in the proof of Lemma A.1,
$\left[D_{x}\BFx^{G}(t,\BFx)(\BFy)\right]_{v,s^{\prime}}$ refers to the
directional derivative of $\BFx^{G}$ in its initial condition in direction
$\BFy$ evaluated at $(v,s^{\prime})$. By the properties of $\BFx^{G}$, as
stated in Lemma 4.1, $D_{x}x^{G}$ is Lipschitz continuous in with respect to
the initial condition. Thus, we can bound

$\displaystyle\sum_{s^{\prime}}||\left(D_{x}x^{G}\bigl{(}\tau,\BFX+\delta(\BFe^{N}_{(k,s_{k})}-\BFe^{N}_{(k,s_{k})}))-D_{x}x^{G}_{u,s}\bigl{(}\tau,\BFX\bigr{)}\right)(\BFe^{N}_{(k,s_{k}^{\prime})}-\BFe^{N}_{(k,s_{k})})||_{L_{2}}$
$\displaystyle\leq\sum_{s^{\prime}}L_{D_{x(0)}x^{G}}\left\|(\BFe^{N}_{(k,s_{k})}-\BFe^{N}_{(k,s_{k})})\right\|_{L_{2}}$

Application of Cauchy-Schwarz then yields

$\displaystyle
R_{1}(x^{G}_{u,s}(t,.),\BFX,\BFe^{N}_{(k,s_{k}^{\prime})}-\BFe^{N}_{(k,s_{k})})$
$\displaystyle\leq\sum_{k,s_{k},s_{k}^{\prime}}L_{D_{x}x^{G}}\left\|(\BFe^{N}_{(k,s_{k})}-\BFe^{N}_{(k,s_{k})}\right\|^{2}_{L_{2}}\left|\bigl{(}X^{G^{N}}_{k/N,s_{k}}r^{N,\text{uni}}_{k,s_{k}\to
s_{k}^{\prime}}+X^{G^{N}}_{k/N,s_{k}}\sum_{s_{l}\in\mathcal{S}}\int_{0}^{1}r^{\text{pair}}_{k/N,v,s_{k}\to
s_{k}^{\prime},s_{l}}G^{N}_{k/N,v}X^{G^{N}}_{v,s_{l}}dv\right|.$

By the boundedness of the rates and the finite state space, it is immediate to
see that terms of the form

$\displaystyle X^{G^{N}}_{k/N,s_{k}}r^{\text{uni}}_{k,s_{k}\to
s_{k}^{\prime}}+X^{G^{N}}_{k/N,s_{k}}\sum_{s_{l}\in\mathcal{S}}\int_{0}^{1}r^{\text{pair}}_{k/N,v,s_{k}\to
s_{k}^{\prime},s_{l}}G^{N}_{k/N,v}X^{G^{N}}_{v,s_{l}}dv$

are bounded. Additionally, by definition of $e^{(k,s_{k})}$ which are non-zero
only on a sub-interval of $(0,1]$ of size $\frac{1}{N}$,
$\left\|e^{N}_{(k,s_{k})}-\BFe^{N}_{(k,s_{k})}\right\|^{2}_{L_{2}}\leq 2/N$.
Lastly, application of the mentioned properties and bounds show that there
exists a constant $C_{\tilde{R}}>0$ such that
$\tilde{R}_{1}\bigl{(}x^{G}_{u,s}(t,\BFX)\bigr{)}=C_{\tilde{R}}/N=\mathcal{O}(1/N)$.

###  A.3 Bound on the Difference of Stochastic and Deterministic Drift (Lemma
A.5)

Lemma A.5 is concerned with bounding the difference between the drift of the
stochastic and deterministic system. As stated in the Lemma, the bound depends
mainly on the $L_{2}$ difference between rates and between the graphon and
graph of the stochastic system.

###### Lemma A.5

Let $D_{x}x^{G}(t,\BFX)(\BFy)$ be the directional derivative of $\BFx^{G}$
with respect to the initial condition in direction $\BFy$ in $L_{2}$. Then,

$\displaystyle\int_{0}^{1}\sum_{s^{\prime}}\left[D_{x}x^{G}(t,\BFX)\left(\BFF^{G}\bigl{(}\BFx^{G}(t,\BFX)\bigr{)}-\BFF^{G^{N}}(\BFx^{G}(t,\BFX))\right)\right]_{v,s^{\prime}}dv$

is bounded by
$6C_{D_{x(0)}}\left(\frac{2L_{\BFR^{\text{pair}}}+16C^{2}_{\BFR^{\text{pair}}}K_{\BFR^{\text{pair}}}}{N}-\frac{16C^{2}_{\BFR^{\text{pair}}}K^{2}_{\BFR^{\text{pair}}}}{N^{2}}+C_{\BFR^{\text{pair}}}^{2}\left|\mathcal{S}\right|^{2}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G-G^{N}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|\right)$.

###### Proof A.6

Proof. Implied by the finite values of $\BFX\in\mathcal{X}^{N}$

$\displaystyle\int_{0}^{1}\sum_{s^{\prime}}\left[D_{x}x^{G}(t,\BFX)\left(\BFF^{G}\bigl{(}\BFx^{G}(t,\BFX)\bigr{)}-\BFF^{G^{N}}(\BFx^{G}(t,\BFX))\right)\right]_{v,s^{\prime}}dv$
$\displaystyle\leq
C_{N,D_{x}x^{G}}\int_{0}^{1}\sum_{s^{\prime}}\left(\BFF^{G}_{v,s^{\prime}}\bigl{(}\BFx^{G}(t,\BFX)\bigr{)}-\BFF^{G^{N}}_{v,s^{\prime}}(\BFx^{G}(t,\BFX))\right)^{2}dv.$

where $C_{N,D_{x}x^{G}}$ is essentially obtained by taking the maximum over
all $\BFX\in\mathcal{X}^{N}$. To conclude it remains to be shown that
$\int_{0}^{1}\left(F^{G}_{v,s^{\prime}}\bigl{(}\BFx^{G}(t,\BFX)\bigr{)}-F^{G^{N}}_{v,s^{\prime}}(\BFx^{G}(t,\BFX))\right)^{2}dv$
is small. By definition we are looking at

$\displaystyle\int_{0}^{1}\sum_{s^{\prime}}\left(F^{G}_{u,s^{\prime}}\bigl{(}\BFx^{G}(t,\BFX)\bigr{)}-F^{G^{N}}_{u,s^{\prime}}(\BFx^{G}(t,\BFX))\right)^{2}du$
$\displaystyle=\int_{0}^{1}\sum_{s^{\prime}}\left(\BFR^{\text{uni}}_{u,s^{\prime}}\BFx^{G}_{u}(t,\BFX)+\BFx^{G}_{u}(t,\BFX)^{T}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s^{\prime}}G_{u,v}\BFx^{G}_{v}(t,\BFX)\
dv\right.$
$\displaystyle\quad\left.-\BFR^{N,\text{uni}}_{u,s^{\prime}}\BFx^{G}_{u}(t,\BFX)-\BFx^{G}_{u}(t,\BFX)^{T}\int_{0}^{1}\BFR^{N,\text{pair}}_{u,v,s^{\prime}}G^{N}_{u,v}\BFx^{G}_{v}(t,\BFX)dv\right)^{2}du$
$\displaystyle\leq
2\int_{0}^{1}\sum_{s^{\prime}}\left(\bigl{(}\BFR^{\text{uni}}_{u,s^{\prime}}-\BFR^{N,\text{uni}}_{u,s^{\prime}}\bigr{)}\BFx^{G}_{u}(t,\BFX)\right)^{2}du$
(21)
$\displaystyle\quad+2\int_{0}^{1}\sum_{s^{\prime}}\left(\BFx^{G}_{u}(t,\BFX)^{T}\bigl{(}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s^{\prime}}G_{u,v}\BFx^{G}_{v}(t,\BFX)dv-\int_{0}^{1}\BFR^{N,\text{pair}}_{u,v,s^{\prime}}G^{N}_{u,v}\BFx^{G}_{v}(t,\BFX)dv\bigr{)}\right)^{2}\
du$ (22)

where the inequality follows as for $f,g\in L_{2}(0,1]$ it holds that
$\left\|f+g\right\|^{2}\leq 2(\left\|f\right\|^{2}+\left\|g\right\|^{2})$. We
separately bound the two summands. The application of Lemma A.7 directly gives
a bound for the first

$\displaystyle
2\int_{0}^{1}\sum_{s^{\prime}}\left(\bigl{(}\BFR^{\text{uni}}_{u,s^{\prime}}-\BFR^{N,\text{uni}}_{u,s^{\prime}}\bigr{)}\BFx^{G}_{u}(t,\BFX)\right)^{2}du\leq
2(\frac{2L_{\BFR^{\text{uni}}}+16C^{2}_{\BFR^{\text{uni}}}K_{\BFR^{\text{uni}}}}{N}-\frac{16C^{2}_{\BFR^{\text{uni}}}K^{2}_{\BFR^{\text{uni}}}}{N^{2}}).$

For the second term, first, we add and subtract the term
$\BFx^{G}_{u}(t,\BFX)^{T}\int_{0}^{1}\BFR^{N,\text{pair}}_{u,v,s^{\prime}}G_{u,v}\BFx^{G}_{v}(t,\BFX)dv$
and use again the previously applied inequality to get the two terms

$\displaystyle
4\int_{0}^{1}\sum_{s^{\prime}}\left(\BFx^{G}_{u}(t,\BFX)^{T}\bigl{(}\int_{0}^{1}\BFR^{\text{pair}}_{u,v,s^{\prime}}-\BFR^{N,\text{pair}}_{u,v,s^{\prime}}\bigr{)}G_{u,v}\BFx^{G}_{v}(t,\BFX)dv\right)^{2}\text{
and }$ (i) $\displaystyle
4\int_{0}^{1}\sum_{s^{\prime}}\left(\BFx^{G}_{u}(t,\BFX)^{T}\int_{0}^{1}\BFR^{N,\text{pair}}_{u,v,s^{\prime}}\bigl{(}G_{u,v}-G^{N}_{u,v}\bigr{)}\BFx^{G}_{v}(t,\BFX)\bigr{)}dv\right)^{2}.$
(ii)

At its core, it remains to bound the difference of
$\BFR^{\text{pair}}-\BFR^{N,\text{pair}}$ rate matrices and the difference
$G-G^{N}$ between the graphon and the graph. Similar to the bound for the
unilateral rates, we use Lemma A.7 to bound Equation (i), i.e.,

$\displaystyle\eqref{eq:drift_difference_pair_rates}\leq
4(\frac{2L_{\BFR^{\text{pair}}}+16C^{2}_{\BFR^{\text{pair}}}K_{\BFR^{\text{pair}}}}{N}-\frac{16C^{2}_{\BFR^{\text{pair}}}K^{2}_{\BFR^{\text{pair}}}}{N^{2}}).$

The second Equation (ii) can be bounded by
$4C_{\BFR^{\text{pair}}}^{2}\left|\mathcal{S}\right|^{2}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G-G^{N}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|$.
To obtain this bound, we use the bounds $C_{\BFR^{\text{pair}}}$ for
$\BFR^{N,\text{pair}},\BFR^{\text{pair}}$, the fact that $\BFR^{\text{pair}}$
is a $\mathcal{S}\times\mathcal{S}$ matrix, and lastly that
$\left\|x^{G}\right\|_{L_{2}}=1$. Application of the bounds yields

$\displaystyle\sum_{s^{\prime}}\langle
D_{x(0)}x^{G}_{u,s}(t,\BFX)_{s^{\prime}},F^{G}_{s^{\prime}}\bigl{(}\BFx^{G}(t,\BFX)\bigr{)}-F^{G^{N}}_{s^{\prime}}(\BFx^{G}(t,\BFX))\rangle_{L_{2}}$
$\displaystyle\leq
6C_{D_{x(0)}}\left(\frac{2L_{\BFR^{\text{pair}}}+16C^{2}_{\BFR^{\text{pair}}}K_{\BFR^{\text{pair}}}}{N}-\frac{16C^{2}_{\BFR^{\text{pair}}}K^{2}_{\BFR^{\text{pair}}}}{N^{2}}+C_{\BFR^{\text{pair}}}^{2}\left|\mathcal{S}\right|^{2}\left|\kern-1.07639pt\left|\kern-1.07639pt\left|G-G^{N}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|\right)$

which concludes the proof.

###  A.4 Bound on the Approximation of a Piecewise Lipschitz Function (Lemma
A.7)

In several part of the paper, the difference between a piecewise Lipschitz
function and a discretized version needs to be quantified. The next lemma
states an upper bound depending on the coarseness of the discretization as
well as properties induces by the definition of piecewise Lipschitz
continuity.

###### Lemma A.7

Let $B$ be a bounded piecewise Lipschitz function with $K\in\mathbb{N}$
blocks, Lipschitz constant $L_{B}<\infty$ as defined in 3.1 and bound
$C_{B}\geq 0$ such that $\left|B(u,v)\right|\leq C_{B}$. Furthermore, let
$B^{N}$ be its discretized version defined by
$B^{N}(u,v)=\sum_{i,j=1}^{N}B(i/N,j/N)\mathds{1}_{I^{N}_{i}}(u)\mathds{1}_{I^{N}_{j}}(v)$.
Assuming that $N$ is large enough as in Definition 3.3, it holds for $f\in
L_{2}(0,1]$ with $\left\|f\right\|_{L_{2}}\leq 1$ that

$\displaystyle\int_{0}^{1}\left|\int_{0}^{1}\bigl{(}B(u,v)-B^{N}(u,v)\bigr{)}f(v)dv\right|du\leq\frac{2L_{B}+16C_{B}^{2}K}{N}-\frac{16C_{B}^{2}K^{2}}{N^{2}}.$

Note, in the case of $B$ being Lipschitz continuous, i.e., $K=0$, the right-
hand side reduces to the anticipated bound of $2L_{B}/N$.

###### Proof A.8

Proof. By Cauchy-Schwarz inequality on $L_{2}(0,1]$ it follows that

$\displaystyle\int_{0}^{1}\bigl{(}B(u,v)-B^{N}(u,v)\bigr{)}f(v)dv$
$\displaystyle\leq\int_{0}^{1}\bigl{(}B(u,v)-B^{N}(u,v)\bigr{)}^{2}dv\left\|f\right\|^{2}_{L_{2}}$
$\displaystyle\leq\int_{0}^{1}\bigl{(}B(u,v)-B^{N}(u,v)\bigr{)}^{2}dv.$

For the rest of the proof we refer to Avella-Medina et al. (2018) as the
reasoning is identical as for their proof of Theorem 1, staring at Equation
(31). The sole adaptation made is to include the case that $B$ can take values
outside of $(0,1]$, leading to $2C_{B}$ as the bound of the difference
$\left|B(u,v)-B^{N}(u,v)\right|$.

###  A.5 Taylor Expansion with Remainder for Banach Spaces (Lemma A.9)

For completeness, we reformulate the result from Blanchard and Brüning (2015)
pp. 524-525.

###### Lemma A.9

Suppose $E,F$ are real Banach spaces, $U\subset E$ an open and nonempty subset
and $f\in\mathcal{C}^{n}(U,F)$ ($n$-times continuously differentiable). Given
$x_{0}\in U$ choose $r>0$ such that $x_{0}+B_{r}\subset U$, where $B_{r}$ is
the open ball in $E$ with center $0$ and radius $r$. Then for all $\ell\in
B_{r}$ we have, using the abbreviation $(\ell)^{k}=(\ell,...,\ell)$,

$\displaystyle
f(x_{0}+\ell)=\sum_{k=0}^{n}\frac{1}{k!}f^{(k)}(x_{0})(\ell)^{k}+R_{n}(f,x_{0},\ell),$
(23)

where the remainder $R_{n}$ has the form

$\displaystyle
R_{n}\left(f,x_{0},\ell\right)=\frac{1}{(n-1)!}\int_{0}^{1}(1-\theta)^{n-1}\left[f^{(n)}\left(x_{0}+\theta\ell\right)-f^{(n)}\left(x_{0}\right)\right](\ell)^{n}d\theta.$
(24)