Accuracy of the Graphon Mean Field Approximation for Interacting Particle Systems

1. Introduction

Mean field approximations are ubiquitously used in the study of large-scale stochastic systems. Its mathematical foundations date back to the 1960s and 1970s with the seminal papers of Kurtz (1970, 1971), McKean (1967), and Norman (1972). Whereas 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, for example, 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 classic 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 classic mean field methods complicated if not infeasible.

To overcome these difficulties, mean field models that use graphons as a description of diversity of individuals and the variation of their interaction have been proposed (Budhiraja et al. 2019, Vizuete et al. 2020, Bayraktar et al. 2023, Keliger 2023, Bet et al. 2024). One example for when this approach finds widespread application are epidemic models for which graphons characterize the epidemic propagation (Decreusefond et al. 2012, Aurell et al. 2022, Delmas et al. 2023, 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. Whereas 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 that 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 ${\mathbf{x}}^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_k^{G^N}(t)=s)$ this probability, our main result shows that

where $⦀G^N-G{⦀}_{L_2}$ is the $L_2$ distance of the step graphon representation of $G^N$ and the graphon G. The $L_2$ norm that we use is essentially equivalent to the distance implied by graphon cut norm (see the discussion in appendix B of Avella-Medina et al. 2018).

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_{kl}^N≔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 a system composed of heterogeneous particles. In the second case (which we call stochastic sampling), $G^N$ is a random graph generated from G, where an edge is present between two nodes k and $\mathcal{\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 nonuniform 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

where, 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 Corollaries 1 and 2. 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, at which each server similarly acts as a dispatcher and keeps or forwards the job according to the join-the-shortest-of-two-queues (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, whereas the numerical complexity of the approximation only slightly increases compared with 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. In Section 5, we present the two numerical examples. Finally, we conclude in Section 7 in which we discuss limitations and future work.

1.3. Related Work

1.3.1. 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 (e.g., Abbe 2018, Bhamidi et al. 2019, Bayraktar and Wu 2021, Bet et al. 2024, 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) and 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) and Aurell et al. (2022). In the not-so-dense setting, available results are more limited with some newer references being Bayraktar et al. (2023) and 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 (2024), Ramanan (2022), and Ganguly (2022).

1.3.2. 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), and Budhiraja et al. (2019). Here, the authors study a variety of load-balancing systems with dynamics based on compatibility or locality constraints, and these give rise to intricate connectivity between dispatchers and servers. Whereas 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.

1.3.3. Generator and Stein’s Method.

For our proofs, we adapt techniques used in Gast and Van Houdt (2017), Allmeier and Gast (2022), and Gast (2017), which, in turn, rely on the use of Stein’s (1986) method. The method is used to estimate and bound the distance between two random variables through their respective generators. Since the publication of works Braverman et al. (2017) and Braverman and Dai (2017), Stein’s (1986) 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_{kl}^N$ 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 kth object at time $t\ge 0$ is denoted by $S_k^{G^N}(t)$. As we see later, $G^N$ can correspond to a random graph for which $G_{kl}^N\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 ${\mathbf{S}}^{G^N}(t)≔(S_1^{G^N},\dots ,S_N^{G^N})(t)\in {\mathcal{S}}^N$. We assume that the process ${\mathbf{S}}^{G^N}≔{({\mathbf{S}}^{G^N}(t))}_{t\ge 0}$ is a continuous time Markov chain 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 two ways:

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

• (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}_{kl}^N/N$.

Note that the rate functions are assumed to be heterogeneous; that is, 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. We show in Example 5.1 that our result can also be used if the scaling factor depends not on the system size N but on the average degrees. We further want to point out that we restrict our framework to the interaction of two particles, and these can be utilized to model many relevant interacting particle systems on graphs such as, for example, 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 between more than two particles by modeling, for instance, the interaction between more than two particles as hypergraphs. 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 use a binary representation of the state based on indicator functions. We denote the new representation by ${\mathbf{X}}^{G^N}={(X^{G^N}{}_{k,s}(t))}_{\stackrel{k\in [N]\hspace{0.17em}s\in \mathcal{S}}{t\ge 0}}$, where

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

Whereas 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 ${\mathbf{e}}_{k,s}^N$ a matrix of size $N\times |\mathcal{S}|$ whose $(k,s)$ component is equal to one, 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}+{\mathbf{e}}_{k,s_k}^N-{\mathbf{e}}_{k,s_k^{\prime }}^N$ at rate

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 ${\mathbf{X}}^{G^N}$, we define what we call the drift of the system of N particles as the expected change for ${\mathbf{X}}^{G^N}$ in state $\mathbf{X}\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 ${\mathbf{F}}^{G^N}(\mathbf{X})$. By using Equation (2), it is equal to

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

In the following, it is 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

where ${\mathbf{R}}_{k,s}^{N,\text{uni}}$ is a row vector whose $s^{\prime }$ component is $r_{k,s^{\prime }\to s}^{N,\text{uni}}$ if $s^{\prime }\ne s$ and $-{\sum }_{\tilde{s}}{r}_{k,s\to \tilde{s}}^{N,\text{uni}}$ for $s=s^{\prime }$ and ${\mathbf{R}}_{k,l,s}^{N,\text{pair}}$ 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 }\ne 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 ${\mathbf{X}}^{G^N}$, $G^N$, and ${\mathbf{F}}^{G^N}$ as Functions from $(0,1]$

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

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_{uv}^N=G_{kl}^N$ for $u\in ((k-1)/N,k/N]$ and $v\in ((l-1)/N,l/N]$, $F_{u,s}^{G^N}=F_{k,s}^{G^N}$, and ${\mathbf{R}}_{k,s}^{N,\text{uni}}={\mathbf{R}}_{u,s}^{N,\text{uni}}$. By using this notation, the sum of $k\in [N]$—for instance in (3)—can be replaced by an integral, that is, $F_{k,s}^{G^N}(\mathbf{X})=F_{u,s}^{G^N}(\mathbf{X})={\mathbf{R}}_{u,s}^{N, \text{uni}}{\mathbf{X}}_u+{\mathbf{X}}_u^T{\int }_0^1{\mathbf{R}}_{u,v,s}^{N, \text{pair}}{\mathbf{X}}_v{G}_{uv}^{N}dv$ for $u\in ((k-1)/N,k/N]$.

2.5. Notations

Throughout the paper, matrices and vectors are written in bold letters, that is, $\mathbf{X},\mathbf{x},\dots$, and regular letters are used to denote scalars such as $X_{u,s},r_{u,s\to s^{\prime }}^{\text{uni}}$. The indices $s,s_k,s_k^{\prime },s_l,\dots$ 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\le \frac{C}{N}$. Such a constant C might depend on the model considered (the graphon or the rate functions). Most of our results are 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 deal with vectors of $L_2(0,1]$ functions, that is, $f,g\in L_2{(0,1]}^{\mathcal{S}}$. For two vectors g, f we denote the scalar product by ${〈f,g〉}_{L_2(0,1]}={\sum }_{s\in \mathcal{S}}{\int }_0^1{f}_{u,s}{g}_{u,s}\hspace{0.17em}du$ and the induced norm ${\Vert f\Vert }_{L_2}=\sqrt{{\sum }_{s\in \mathcal{S}}{\int }_0^1{f}_{u,s}^{2}du}$. For a function $G:{(0,1]}^2\to \mathbb{R}$ denote by Gf the function defined as ${(Gf)}_{u,s}={\int }_0^1{G}_{uv}{f}_{v,s}\hspace{0.17em}dv$. We denote the $L_2$ norm-operator of G as

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 function² $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_{kl}^N$ 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 a 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 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 nonoverlapping intervals ${\mathcal{A}}_k=(a_{k-1},a_k]$ with $0=a_0<a_1<\cdots <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}$,

$$|G_{uv}-G_{u^{\prime }{v}^{\prime }}|\le L_G(|u-u^{\prime }|+|v-v^{\prime }|).$$

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.

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, that is, $G_{kl}^N=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, that is, $G_{kl}^N=\text{Bernoulli}(G(k/N,l/N))$, and $G_{lk}^N=G_{kl}^N$ 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 a finite graph $G^N$ converges to G), we use the $L_2$ operator norm as defined in Section 2.5. More precisely, for a measurable function $f:(0,1]\to \mathbb{R}$, we denote by ${\Vert f\Vert }_{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 $⦀G⦀$ the operator norm of G in $L_2$, that is

The distance between two graphons (or finite graphs) G and $G^{\prime }$ is measured as $⦀G-G^{\prime }⦀$.

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 2) 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 1

(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 Section 3.2. Then, for large enough N (as in Definition 2) 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

$$⦀G-G^N{⦀}_{L_2}\le \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).$$(4)

Definition 2

(Large Enough N (Avella-Medina et al. 2018). Given a piecewise Lipschitz graphon G with partition of $(0,1]$ of nonoverlapping intervals ${\mathcal{A}}_k=(a_{k-1},a_k]$ as in Definition 1 and $\delta <e^{-1}$. The quantity N is called large enough if

$$\frac{2}{N}<\underset{k\in \{0..K_G\}}{\min }(a_k-a_{k-1}),$$(5a)

$$\frac{1}{N}\log \left(\frac{2N}{\delta }\right)+\frac{2K_G+3L_G}{N}<\underset{u\in (0,1]}{\sup }{\displaystyle {\int }_0^1{G}_{u,v}}\hspace{0.17em}dv\hspace{1em}\text{and}$$(5b)

$$N{e}^{-N/5}<\delta .$$(5c)

Implied by Theorem 1, 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 $⦀G-G^N{⦀}_{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 now state the assumptions that we use to state the accuracy of the graphon mean field approximation.

• (A₁) The state space $\mathcal{S}$ is finite.

• (A₂) The graphon G is piecewise Lipschitz continuous.

• (A₃) 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:

  $$r_{k/N,s_k\to s_k^{\prime }}^{\text{uni}}=r_{k,s_k\to s_k^{\prime }}^{N,\text{uni}}\hspace{1em} \text{and} \hspace{1em}{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}}\hspace{1em} \text{for} k,l\in [N].$$(6)

The first assumption on the finiteness of the state-space is not necessary per se but simplifies the exposition of the results. Our results can be generalized to countable state-space $\mathcal{S}$ but at the price of extra assumptions that would clutter the readability of the paper. Imposing a finite state-space simplifies the definition of the $L_2$ space. To allow for nonfinite $\mathcal{S}$, one needs to add not just a uniform bound on the transition rates (such as $(A_3)$) but a bound on their norm in this infinite-dimensional space: in our proof, some bounds depend on $|\mathcal{S}|$ (see, for instance, Online Lemma 7) because $(A_3)$ only provides a bound on the infinite-norm of the rate functions. When $\mathcal{S}$ is finite, all such norms are equivalent. The second assumption is a regularity assumption that we use to obtain bounds on the rate at which a sequence of graphs converges to graphons. As shown in Lovász (2012, lemma 11.33), some regularity assumptions on G are needed to show the convergence. In practice, our bounds 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_3$) can be generalized by replacing the equality in (6) with bounds for the distance of $r^{N,·}$ and $r^{·}$, therefore modifying the bounds of the theorem to include terms of the form $\Vert r^{N,·}-r^{·}\Vert$.

4.2. The Graphon Mean Field Approximation

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

where ${\mathbf{R}}^{\text{uni}}$ and ${\mathbf{R}}^{\text{pair}}$ are the functions corresponding to $r^{\text{uni}}$ and $r^{\text{pair}}$ as we did for ${\mathbf{R}}^{N,\text{uni}}$ and ${\mathbf{R}}^{N,\text{pair}}$ in Section 2.4. This equation is identical to the original drift Equation (3) except that we replace the sum over discrete variables by integrals over $u,v\in (0,1]$.

Based on the drift function ${\mathbf{F}}^G$ and the initial condition ${\mathbf{x}}_0$, we call the graphon mean field approximation the solution of the differential equation

Lemma 2.

Let ${\mathbf{F}}^G$ denote the deterministic drift defined in Equation (7) and ${\mathbf{x}}^G(t,{\mathbf{x}}_0)$ the value of the solution of the graphon mean field approximation at time t with initial condition ${\mathbf{x}}_0$ as in Equation (8). It holds that, for all t and all initial conditions ${\mathbf{x}}_0$, ${\mathbf{x}}^G(t,{\mathbf{x}}_0)$ is well-defined and unique (in $L_2$) 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}[{\mathbf{X}}^{G^N}]$ and ${\mathbf{x}}^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\ge 0$ is approximated by the graphon mean field approximation ${\mathbf{x}}^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 1

($L_2$ Convergence). Let ${\mathbf{X}}^{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 ${\mathbf{x}}^G(t)={(x^G{}_{u,s}(t,{\mathbf{x}}_0))}_{u\in (0,1],s\in \mathcal{S}}$ be the mean field approximation of the particle system as defined in Equation (8) for an initial condition ${\mathbf{x}}_0={\mathbf{X}}^{G^N}(0)$.³ Assume additionally that conditions ($A_1$)–($A_3$) are fulfilled and let $t>0$ be arbitrary but fixed. Then, there exist constants $C_A,C_B\ge 0$ such that

$${\Vert \mathbb{E}[{\mathbf{X}}^{G^N}(t)|{\mathbf{X}}^{G^N}(0),G^N]-{\mathbf{x}}^G(t)\Vert }_{L_2}\le \frac{C_A}{N}+C_B⦀G^N-G⦀.$$(9)

The proof of Theorem 1 is postponed to Section 6.2. The constants $C_A,C_B$ of Equation (9) depend on the uniform bound of the rates, the time t, and the Lipschitz constants of the graphon. In the subsequent Corollaries 1 and 2, we see that, if the $G^N$ is generated by one of the methods illustrated in Section 3.2, precise bounds on the distance $⦀G^N-G⦀$ 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 1 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 $⦀G^N-G⦀$ 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. Yet, because of our generic choice of transition rates, this would be linked to increasingly heavy notations for the dynamics, only giving little more insight into the accuracy of the graphon mean field method.

A nice interpretation of Theorem 1 is in terms of total variation distance. Indeed, by definition of $X_s^{G^N}(t)$ as an indicator function (1), we have $\mathbb{E}[X_{s,u}^{G^N}(t)]=\mathbb{P}(S_{\lfloor Nu\rfloor }(t)=s)$. The total variation distance between two distributions $\mu$ and $\nu$ is equal to $\frac{1}{2}{\sum }_{s\in \mathcal{S}}|{\mu }_s-{\nu }_s|=\frac{1}{2}{\Vert \mu -\nu \Vert }_1$. Hence, by the Cauchy–Schwarz inequality, one has

The above equation shows that, on average, the law of a particle k converges to $x_{k/N,·}^G(t)$. Note that this does not mean that the law of each particle k converges to $x_{k/N,·}^G(t)$. For instance, changing all edges of a single particle only changes $G^N$ by a factor $O(1/N)$.

4.4. Case-Specific Bounds for $⦀G^N-G{⦀}_{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 1 shows that the accuracy is of order $\mathcal{O}(1/N)$. Our second Corollary 2 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}\left(\sqrt{\frac{\log N}{N}}\right)$.

4.4.1. Case 1: Graphon Discretization.

The corollary provides a bound on the accuracy of 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 1

(Dense Heterogeneous System). Assume ($A_1$)–($A_3$) and let ${\mathbf{x}}^G$ and ${\mathbf{X}}^{G^N}$ be defined as in Theorem 1. Let $k\in [N],s\in \mathcal{S}$ and $t\ge 0$ be arbitrary but fixed. If $G^N$ is generated by the deterministic sampling method of Section 3.2, that is, a discretization of the graphon G, it holds that

$${\Vert \mathbb{E}[{\mathbf{X}}^{G^N}(t)|{\mathbf{X}}^{G^N}(0),G^N]-{\mathbf{x}}^G(t,{\mathbf{X}}^{G^N}(0))\Vert }_{L_2}\le \frac{C_A+C_B{C}_{G^N}}{N}.$$(10)

The proof is postponed to Section 6.3. The constants $C_A,C_B$ are as in Theorem 1. 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) and Gast and Van Houdt (2017), 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 2

(Graphon System Approximation). Assume the conditions ($A_1$)–($A_3$) as in Theorem 1 for a symmetric graphon G. Let ${\mathbf{X}}^{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 Section 3.2. Let $t\ge 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 2, the graph $G^N$ is sampled such that

$${\Vert \mathbb{E}[{\mathbf{X}}^{G^N}(t)|{\mathbf{X}}^{G^N}(0),G^N]-{\mathbf{x}}^G(t,{\mathbf{X}}^{G^N}(0))\Vert }_{L_2}\le \frac{C_A}{N}+C_B{\psi }_G(N),$$(11)

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 (1).

In the above theorems, the constants $L_G$ and $K_G$ are a measure of heterogeneity of the nodes of the graphon. If G represents a stochastic bloc model, then $L_G=0$ and $K_G$ corresponds to the number of blocks. If the graphon represents a continuous popularity (as in Figure 1 or in our bike-sharing example), then $K_G=0$. A graphon with a large $K_G$ or a large $L_G$ would represent a graph in which the nodes have a very heterogeneous behavior.

The proof of the corollary is postponed to Section 6.3.

5. Numerical Experiments

In this section, we present two examples that 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 because of the 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 in which the server–dispatcher pairs are connected through a graph structure in which 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_{kl}^N$ 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–dispatcher pair following a Poisson arrival process with rate ${\lambda }_L>0$. Arriving jobs are distributed according to the JSQ(2) policy idea; that is, 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 ${\mathbf{X}}^{G^N}={(X^{G^N}{}_{k,s})}_{k\in [N],s=0..K_L}$ the transitions of the system are

where ${\mathbf{e}}_{k,s}$ is a $N\times |\mathcal{S}|$ 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 (14) corresponds to the completion of a job by server k when the queue is of size $1\le s\le K_L$. The second type of transition, Equation (13), corresponds to adding a job to k having $0\le s\le 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 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, that is, has no edges, we define $\frac{G_{kl}^N}{d^{(N)}(l)}$ to be zero.

5.1.2. Drift and Graphon Mean Field Approximation.

For the deterministic drift, 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 $\mathbf{x}={(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

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

5.1.3. Derivation of Accuracy Bounds.

Whereas edges between the server–dispatcher pair are sampled according to the stochastic sampling method of Section 3.2, the dependence of the rates on the node degree prevents a direct application of the results of Corollary 2. In particular, to apply Corollary 2, it is assumed that the graph edges are scaled by a factor of $1/N$ instead of the node degree. Yet we can obtain accuracy bounds similar to the one given by the corollary by changing the definition of the graph $G^N$.

Corollary 3.

In the load-balancing example, there exists $C>0$ such that, with probability at least $1-2/N$ and for large enough N, the graph $G^N$ is sampled such that

$${\Vert \mathbb{E}[{\mathbf{X}}^{G^N}(t)|{\mathbf{X}}^{G^N}(0),G^N]-{\mathbf{x}}^G(t,{\mathbf{X}}^{G^N}(0))\Vert }_{L_2}\le C\sqrt{\frac{\log (N)}{N}}.$$(15)

Proof.

For a sketch of the proof, we start by defining ${\tilde{G}}_{kl}^N≔(\frac{G_{kl}^N}{d^{(N)}(k)}+\frac{G_{kl}^N}{d^{(N)}(l)})\frac{N}{16}$, $r_{k,l,s\to s+1,s_l}^{N,\text{pair}}={\lambda }_L(\frac{1}{2}{\mathrm{𝟙}}_{s_l=s}+{\mathrm{𝟙}}_{s_l>s}){\mathrm{𝟙}}_{s<K_L}$, and $r_{k,s\to s-1}^{N,\text{uni}}={\mu }_L{\mathrm{𝟙}}_{s>0}$ with $r^{N,\text{pair}}$ being the pairwise transitions and $r^{N,\text{uni}}$ the unilateral transitions for $s,s^{\prime },s_l\in \{0,\dots ,K_L\}$. This allows us to obtain the drift of the stochastic system similar to the one introduced in Equation (3), namely

$$F_{k,s}^{G^N}(\mathbf{X})=X_{k,s}{\displaystyle \sum _{s_l\ge s}{\displaystyle \sum _{l\in [N]}1}}6\frac{{\tilde{G}}_{kl}^N}{N}{r}_{k,l,s\to s+1,s_l}^{N,\text{pair}}{X}_{l,s^{\prime }}-X_{k,s}{r}_{k,s\to s-1}^{N,\text{uni}}.$$

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_{k,l,s\to s+1,s_l}^{\text{pair}}≔{\lambda }_L(\frac{1}{2}{\mathrm{𝟙}}_{s_l=s}+{\mathrm{𝟙}}_{s_l>s}){\mathrm{𝟙}}_{s<K_L}$ and $r_{k,s\to s-1}^{N,\text{uni}}≔{\mu }_L{\mathrm{𝟙}}_{s>0}$ to rewrite

$$F_{u,s}^G(\mathbf{x})=x_{u,s}{\displaystyle {\int }_0^1{\displaystyle \sum _{s_v\ge s}1}}6{\tilde{G}}_{uv}{r}_{k,l,s\to s+1,s_v}^{\text{pair}}{x}_{v,s_v}dv-x_{u,s}{r}_{k,s\to s-1}^{N,\text{uni}}.$$

Note that, for the above reformulations, the result of our main Theorem 1 is still applicable and yields constants $C_A,C_B\ge 0$ such that ${\Vert \mathbb{E}[{\mathbf{X}}^{G^N}(t)|G^N]-{\mathbf{x}}^G(t)\Vert }_{L_2}\le \frac{C_A}{N}+C_B⦀{\tilde{G}}^N-\tilde{G}{⦀}_{L_2}$. Using the definition of the ${\tilde{G}}^N,\tilde{G}$, we obtain $⦀{\tilde{G}}^N-\tilde{G}{⦀}_{L_2}\le \frac{1}{16}\left({\Vert G^N-G\Vert }_{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|\right)$. For $⦀G^N-G{⦀}_{L_2}$, we can use the bound as in Corollary 2, 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 2). To obtain similar bounds for the node degree, the multiplicative Chernoff bound can be used, that is, $\mathbb{P}\left(\frac{N}{d^N(k)}\ge \frac{N}{(1-\gamma )\mathbb{E}[d^N(k)]}\right)=\mathbb{P}(d^N(k)\le (1-\gamma )\mathbb{E}[d^N(k)])\le \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}\left(\frac{N}{d^N(k)}\ge \frac{N}{(1-\gamma )\mathbb{E}[d^N(k)]}\right)\le \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|\le \left|{\sum }_{k=1}^N(\frac{N}{d^N(k)}-\mathbb{E}[\frac{N}{d^N(k)}\left]\right){\mathbf{1}}_{u\in (k-1/N,k/N]}\right|+|{\sum }_{k=1}^N\mathbb{E}\left[\frac{N}{d^N(k)}\right]$ ${\mathbf{1}}_{u\in (k-1/N,k/N]}-\frac{1}{d(u)}|$. We apply the Chernoff bound to the first summand to show that the difference concentrates around zero, that is, $\mathbb{P}\left(\frac{N}{d^N(k)}\le \frac{N}{(1-\gamma )\mathbb{E}[d^N(k)]}\right)\ge 1-2/N$, and therefore, $\frac{N}{d^N(k)}-\frac{N}{\mathbb{E}[d^N(k)]}\le 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)}\le \frac{N}{(1-\gamma )\mathbb{E}[d^N(k)]}$ follows using the same arguments. For the second summand, we use the statement of Online Lemma 8 to obtain an error bound of order $\mathcal{O}(1/N)$. This shows that, for large enough N, $⦀{\tilde{G}}^N-\tilde{G}{⦀}_{L_2}$ is of order $\mathcal{O}\left(\sqrt{\frac{\log N}{N}}\right)$. Combining the two terms, we obtain a bound of order $O(1/N)+O(\sqrt{\log (N)/N})=O(\sqrt{\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 more than 2,000 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_B^{\gamma }{(\mathbf{x})}_{u,s}≔{\sum }_{i=1}^{\gamma }{\mathrm{𝟙}}_{U_i^{\gamma }}(u)F_B{(\mathbf{x})}_{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, that is, for a rudimentary implementation using a NumPy ordinary differential equation (ODE) solver, we are able to solve the discretization in a few seconds.⁴ Figures 2 and 3 show the results of our numerical experiments. 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 2. In Figure 3, we compare the values for the fraction of servers with at least s jobs for $s=1,\dots ,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 by Corollary 2 suggested behavior for single particles in Figure 2, the increased accuracy for $N=20,30$ is likely caused by an averaging effect and a well-connected graph instance $G^N$.

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) and 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 among the stations, that is, $\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\ge 0$. For the system, heterogeneity comes from the varying popularity of stations. Hence, we denote by $p_B:(0,1]\mapsto {\mathbb{R}}_{\ge 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; that is, we discretize the graphon based on the system size N. For two stations $k,l\in [N]$, connectivity is, therefore, defined by $G_{kl}^{N,B}≔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,\dots ,S_N)\in {[K_B]}^N$, the bikes move between stations in the following way:

• Customers arrive at a station $k\in [N]$ with rate ${\lambda }_B$ leading to the transitions

  $$S\to S-{\mathbf{e}}_k^N \text{at} \text{rate} {\lambda }_B{\mathrm{𝟙}}_{S_k>0}.$$

• The travel time of bikes between stations is exponentially distributed with rate ${\mu }_B>0$. Hence, 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, that is, $\frac{p_B(k/N)}{{\int }_0^1{p}_B(\nu )d\nu }\left(M-{\sum }_k{S}_k^N\right){\mu }_B$. This implies the transition

  $$S\to S+{\mathbf{e}}_k^N\hspace{1em}\hspace{1em}\hspace{1em}\mathit{at} \mathit{\text{rate}}\hspace{1em}\hspace{1em}\hspace{1em}\frac{p_B(k/N)}{{\int }_0^1{p}_B(\nu )d\nu }(M-{\displaystyle \sum _k{S}_k^N}){\mu }_B{1}_{S_k<K_B}.$$

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

5.2.2. Drift and 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, that is, 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 $\mathbf{x}={(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

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

5.2.3. Implementation and 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 Figures 4 and 5, we calculate the sample mean by averaging more than 7,500 simulations for each system size. For each plot in Figure 4, the mean field trajectory for a 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 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.

6. Proofs

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

6.1. Proof of Lemma 2

Proof.

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

• We show in Lemma 3 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 4, we show that the directional derivatives of $F^G$ are well-defined and Lipschitz continuous. This property ensures that the ${\mathbf{x}}^G$ is also continuously differentiable (for a proof of this property in general Banach spaces, see, for example, Driver 2003). □

Lemma 3

(Local Lipschitz Continuity of the Drift). Let $L_{{\mathbf{F}}^G}=2(C_R|\mathcal{S}|+2C_R{|\mathcal{S}|}^2⦀G⦀)$, where $C_R$ is the bound on the rate functions, ⦀G⦀ the operator norm of the graphon G and $|\mathcal{S}|$ the size of the finite state space. Then, for all ${\mathbf{x}}_{.}={(x_{.,s})}_{s\in \mathcal{S}},{\mathbf{y}}_{.}={(y_{.,s})}_{s\in \mathcal{S}}$ with $x_{.,s},y_{.,s}\in L_2(0,1]$, and ${\Vert x\Vert }_{L^2},{\Vert y\Vert }_{L^2}\le 1$ with $|x_{u,s}|,|y_{u,s}|\le 1$ almost everywhere,⁵ we have

$${\Vert {\mathbf{F}}^G(\mathbf{x})-{\mathbf{F}}^G(\mathbf{y})\Vert }_{L_2}\le L_{{\mathbf{F}}^G}{\Vert \mathbf{x}-\mathbf{y}\Vert }_{L^2}.$$

Proof.

Define ${\widehat{F}}_{u,s}^{\mathit{G}}(\mathbf{x},\mathbf{y})≔{\mathbf{R}}_{u,s}^{\text{uni}}{\mathbf{x}}_u+{\mathbf{x}}_u^T{\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}{G}_{u,v}{\mathbf{y}}_v\hspace{0.17em}dv$, that is, we replace ${\mathbf{x}}_v$ by ${\mathbf{y}}_v$ under the integral. Applying the triangle inequality to the $L_2$ norm gives

$$\begin{array}{ll}{\Vert {\mathbf{F}}^G(\mathbf{x})-{\mathbf{F}}^G(\mathbf{y})\Vert }_{L_2}\hfill & ={\Vert {\mathbf{F}}^G(\mathbf{x})-{\widehat{\mathit{F}}}^{\mathit{G}}(\mathbf{x},\mathbf{y})+{\widehat{\mathit{F}}}^{\mathit{G}}(\mathbf{x},\mathbf{y})-F(\mathbf{y})\Vert }_{L_2}\hfill \\ \hfill & \le {\Vert {\mathbf{F}}^G(\mathbf{x})-{\widehat{\mathit{F}}}^{\mathit{G}}(\mathbf{x},\mathbf{y})\Vert }_{L_2}+{\Vert {\widehat{\mathit{F}}}^{\mathit{G}}(\mathbf{x},\mathbf{y})-{\mathbf{F}}^G(\mathbf{y})\Vert }_{L_2}.\hfill \end{array}$$

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

$$\begin{array}{l}{\Vert {\mathbf{F}}^G(\mathbf{x})-{\widehat{\mathit{F}}}^{\mathit{G}}(\mathbf{x},\mathbf{y})\Vert }_{L_2}=\sqrt{{\displaystyle \sum _{s\in \mathcal{S}}{\displaystyle {\int }_0^1{\left({\mathbf{R}}_{u,s}^{\text{uni}}{\mathbf{x}}_u+{\mathbf{x}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{x}}_{v}dv-{\mathbf{R}}_{u,s}^{\text{uni}}{\mathbf{x}}_u-{\mathbf{x}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{y}}_v\hspace{0.17em}dv\right)}^2}}du}\\ =\sqrt{{\displaystyle \sum _{s\in \mathcal{S}}{\displaystyle {\int }_0^1{\left({\mathbf{x}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}({\mathbf{x}}_v-{\mathbf{y}}_v)dv\right)}^2}}du}\le 2{|\mathcal{S}|}^3{C}_R⦀G⦀{\Vert \mathbf{x}-\mathbf{y}\Vert }_{L_2}\end{array},$$

and

$$\begin{array}{l}{\Vert {\widehat{\mathit{F}}}^{\mathit{G}}(\mathbf{x},\mathbf{y})-{\mathbf{F}}^G(\mathbf{y})\Vert }_{L_2}=\sqrt{{\displaystyle \sum _{s\in \mathcal{S}}{\displaystyle {\int }_0^1{\left({\mathbf{R}}_{u,s}^{\text{uni}}({\mathbf{x}}_u-{\mathbf{y}}_u)+{({\mathbf{x}}_u-{\mathbf{y}}_u)}^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{y}}_v\hspace{0.17em}dv\right)}^2}}du}\hfill \\ \le \left(2C_R{|\mathcal{S}|}^2+2{|\mathcal{S}|}^3{C}_R⦀G⦀\right){\Vert \mathbf{x}-\mathbf{y}\Vert }_{L_2}.\hfill \end{array}$$

Last, we conclude that ${\Vert F^G(\mathbf{x})-F^G(\mathbf{y})\Vert }_{L_2}\le (2(C_R{|\mathcal{S}|}^2+2C_R{|\mathcal{S}|}^3⦀G⦀)){\Vert \mathbf{x}-\mathbf{y}\Vert }_{L_2}$, where the ${|\mathcal{S}|}^2$ and ${|\mathcal{S}|}^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 4

(Lipschitz Continuous Directional Derivative of ${\mathbf{F}}^G$). The directional derivative of the drift $F^G$ for $\mathbf{x}={(x_s)}_{s\in \mathcal{S}}$ with $x_s\in L_2(0,1]$ is given by

$$D_{\mathbf{x}}{F}_{u,s}^G(\mathbf{x})(\mathbf{h})≔{\mathbf{R}}_{u,s}^{\text{uni}}{\mathbf{h}}_u+{\mathbf{h}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{x}}_{v}dv+{\mathbf{x}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{h}}_{v}dv,$$(16)

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

Proof.

Recall the definition of $F_{u,s}(\mathbf{x})={\mathbf{R}}_{u,s}^{\text{uni}}{\mathbf{x}}_u+{\mathbf{x}}_u^T{\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}{G}_{u,v}{\mathbf{x}}_v\hspace{0.17em}dv$. We define $D_{\mathbf{x}}{F}_{u,s}^G(\mathbf{x})(\mathbf{h})$ for $(u,s)\in (0,1]\times \mathcal{S}$ as in Equation (16). To see that $D_{\mathbf{x}}{\mathbf{F}}^G(x)(h)$ is a directional derivative of ${\mathbf{F}}^G$ in $\mathbf{x}\in L_2{(0,1]}^{\mathcal{S}}$ in direction $h\in L_2$, we show that it fulfills

$$\underset{ϵ\downarrow 0}{\lim }\frac{1}{ϵ}{\Vert {\mathbf{F}}^G(\mathbf{x}+ϵ\mathbf{h})-{\mathbf{F}}^G(\mathbf{x})-D_{\mathbf{x}}{\mathbf{F}}^G(\mathbf{x})(ϵ\mathbf{h})\Vert }_{L_2}=0.$$

Fix $ϵ>0$. By definition,

$$\begin{array}{l}{\Vert {\mathbf{F}}^G(\mathbf{x}+ϵ\mathbf{h})-{\mathbf{F}}^G(\mathbf{x})-D_{\mathbf{x}}{\mathbf{F}}^G(\mathbf{x})(\mathbf{h})\Vert }_{L_2}^2\hfill \\ ={\displaystyle {\int }_0^1{\displaystyle \sum _s({\mathbf{R}}_{u,s}^{\text{uni}}({\mathbf{x}}_u+ϵ{\mathbf{h}}_u)+{({\mathbf{x}}_u+ϵ{\mathbf{h}}_u)}^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}({\mathbf{x}}_v+ϵ{\mathbf{h}}_v)\hspace{0.17em}dv}}-{\mathbf{R}}_{u,s^{\prime }}^{\text{uni}}{\mathbf{x}}_u-{\mathbf{x}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s^{\prime }}^{\text{pair}}}{G}_{u,v}{\mathbf{x}}_v\hspace{0.17em}dv\hfill \\ \hspace{1em}\hspace{1em}{-ϵ{\mathbf{R}}^{\text{uni}}{\mathbf{h}}_u-ϵ{\mathbf{h}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{x}}_{v}dv-ϵ{\mathbf{x}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{h}}_{v}dv)}^{2}du\hfill \\ ={\displaystyle {\int }_0^1{\displaystyle \sum _s(ϵ{\mathbf{R}}_{u,s}^{\text{uni}}{\mathbf{h}}_u+ϵ{\mathbf{x}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{h}}_{v}dv+ϵ{\mathbf{h}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{x}}_v\hspace{0.17em}dv}}+{ϵ}^2{\mathbf{h}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{h}}_v\hspace{0.17em}dv\hfill \\ \hspace{1em}\hspace{1em}{-ϵ({\mathbf{R}}^{\text{uni}}{\mathbf{h}}_u-{\mathbf{h}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{x}}_{v}dv-{\mathbf{x}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{h}}_{v}dv))}^{2}du\hfill \\ ={\displaystyle {\int }_0^1{\displaystyle \sum _s{\left({ϵ}^2{\mathbf{h}}_u^T{\displaystyle {\int }_0^1{\mathbf{R}}_{u,v,s}^{\text{pair}}}{G}_{u,v}{\mathbf{h}}_{v}dv\right)}^2}}du.\hfill \end{array}$$

Dividing by $ϵ$ and taking $ϵ\to 0$ shows that $D_{\mathbf{x}}{\mathbf{F}}^G(\mathbf{x})(\mathbf{h})$ is indeed a directional derivative of ${\mathbf{F}}^G$ in $\mathbf{x}$. By definition, $D_{\mathbf{x}}{\mathbf{F}}_{u,s}^G(\mathbf{x})(\mathbf{h})$ is linear in $\mathbf{x}$ and as ${\mathbf{R}}_{u,s}^{\text{uni}},{\mathbf{R}}_{u,v,s}^{\text{pair}}$, and $G_{u,v}$ are bounded functions, the Lipschitz continuity in $L_2$ follows. □

6.2. Proof of Theorem 1

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

Proof.

Let ${\mathbf{X}}_t^{G^N}$ be the $L_2{(0,1]}^{|\mathcal{S}|}$ representation of the stochastic system as defined in Section 2.4. For a fixed $t\ge 0$, define

$${\mathit{\nu }}^N(\tau )=\mathbb{E}[{\mathbf{x}}^G(t-\tau ,{\mathbf{X}}^{G^N}(\tau ))]$$(17)

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

$$\begin{array}{l}{\Vert \mathbb{E}[{\mathbf{X}}^{G^N}(t)]-x^G(t,{\mathbf{X}}^{G^N}(0))\Vert }_{L_2}=\sqrt{{\displaystyle \sum _{s\in \mathcal{S}}{\displaystyle {\int }_0^1{(\mathbb{E}[X^{G^N}{}_{u,s}(t)]-x^G{}_{u,s}(t,{\mathbf{X}}^{G^N}(0)))}^2}}du}\hfill \\ =\sqrt{{\displaystyle \sum _{s\in \mathcal{S}}{\displaystyle {\int }_0^1{({\nu }_{u,s}^N(t)-{\nu }_{u,s}^N(0))}^2}}du}.\hfill \end{array}$$(18)

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 }_{u,s}^N(\tau )$, with $\frac{d}{d\tau }{\nu }^N(\tau )$ being the quantity that fulfills ${\mathit{\nu }}^N(b)-{\mathit{\nu }}^N(a)={\int }_a^b\frac{d}{d\tau }{\mathit{\nu }}^N(\tau )d\tau$ for arbitrary $a,b\in [0,t]$ with $a<b$. In Online Lemma 5, we show that $\frac{d}{d\tau }\mathit{\nu }(\tau )$ exists almost everywhere for $\tau \in (0,t)$ and is almost everywhere equal to

$$\mathbb{E}\left[{\displaystyle \sum _{s^{\prime }}{\displaystyle {\int }_0^1{\left[D_x{\mathbf{x}}^G(t-\tau ,{\mathbf{X}}^{G^N}(\tau ))({\mathbf{F}}^G({\mathbf{X}}^{G^N}(\tau ))-{\mathbf{F}}^{G^N}({\mathbf{X}}^{G^N}(\tau )))\right]}_{v,s^{\prime }}}}dv\right]$$(19)

$$+\mathbb{E}\left[{\tilde{R}}_1({\mathbf{x}}^G(t-\tau ,{\mathbf{X}}^{G^N}(\tau )))\right]d\tau ,$$(20)

where the term ${\tilde{R}}_1$ is defined in Equation (22) and refers to a Taylor remainder term.

In the above sum, $D_x{\mathbf{x}}^G(\tau ,{\mathbf{X}}^{G^N}(\tau ))({\mathbf{F}}^G({\mathbf{X}}^{G^N}(\tau ))-{\mathbf{F}}^{G^N}({\mathbf{X}}^{G^N}(\tau )))$ is the directional derivative of ${\mathbf{x}}^G$ in its initial condition in direction ${\mathbf{F}}^G({\mathbf{X}}^{G^N}(\tau ))-{\mathbf{F}}^{G^N}({\mathbf{X}}^{G^N}(\tau ))$. We aim to bound the sum by using the properties of the derivative of ${\mathbf{x}}^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 Online Lemmas 6 and 7. From the application of the latter lemmas, one obtains the bounds

$$C_{D_{x(0)}}\left(\frac{2L_{{\mathbf{R}}^{\text{pair}}}+16C_{{\mathbf{R}}^{\text{pair}}}^2{K}_{{\mathbf{R}}^{\text{pair}}}}{N}-\frac{16C_{{\mathbf{R}}^{\text{pair}}}^2{K}_{{\mathbf{R}}^{\text{pair}}}^2}{N^2}+C_{{\mathbf{R}}^{\text{pair}}}^2{|\mathcal{S}|}^2⦀G-G^N⦀\right)$$

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

$${\Vert \mathbb{E}[{\mathbf{X}}^{G^N}(t)]-x^G(t,{\mathbf{X}}^{G^N}(0))\Vert }_{L_2}\le \frac{C_A}{N}+C_B⦀G^N-G⦀$$

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 2 and Corollary 1

In this section, we prove the corollaries of Theorem 1. In each case, we use an additional lemma to obtain a bound on $⦀G^N-G⦀$.

Proof (Proof of Corollary 2).

By the application of Theorem 1, we have

$${\Vert \mathbb{E}[{\mathbf{X}}^{G^N}(t)]-x^G(t,{\mathbf{X}}^{G^N}(0))\Vert }_{L_2}=\frac{C_A}{N}+C_B⦀G^N-G⦀.$$

It, therefore, remains to bound the $L_2$ distance $⦀G^N-G⦀$ between the graphon G and graph $G^N$. By application of Lemma 1 with probability at least $1-\delta$ and for large enough N

$$⦀G^N-G⦀\le \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 (Proof of Corollary 1).

The corollary is a direct implication of Theorem 1 and Online Lemma 8: because of $G^N$ representing a discretized version of G, it directly follows that $⦀G-G^N⦀$ is of order $\mathcal{O}(1/N)$, which concludes the proof. □

7. Conclusion and Discussion

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. Whereas this result is similar to other results in the literature—showing 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 show 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}(⦀G^N-G⦀)$ that depends on the distance between the original graph $G^N$ and the graphon G when measured as an $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 consider 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. In particular, there are two direct follow-ups of this work.

First, the focus of this paper is to study the case of dense graphs that converge to graphons. This implies that the average degree of a node is linear in the number of edges. We believe that our methodology could be adapted in the not-so-dense case in which the average degree, say $d(N)$, goes to infinity at a sublinear rate. To do so, the scaling model should be changed: our model assumes that $G^N$ converges to G and that rates of transitions are of order $O(1/N)$. One needs to consider graphs such that $N{G}^N/d(N)$ converges to the graphon G and to have interaction rates of order $1/d(N)$. Once this is done, the part of our methodology that uses generators (Section 6.2 and Online Lemma 6) would still be applicable but would give us a rate of $O(1/\text{d}(\text{N}))$ instead of $O(1/N)$. We believe that most of the rest of the proof (that concerns the Lipschitz continuity of the graphon mean field approximation) would continue to hold except for Lemma 1 (on stochastic sampling) that requires the average degree to be relatively large. For the case of sparse graphs, the number of neighbors per node remains bounded when the number of nodes N goes to infinity. This requires a fundamentally different approach: studying such a problem is out of the scope of our tools because the graphon mean field approximation is not asymptotically exact for sparse graphs.

Second, in this paper, we restrict our attention to finite-time convergence. We believe that, if the graphon dynamics has a fixed point that attracts all trajectories and one could bound the rate at which the convergence to this fixed point occur, then a bound on the steady-state convergence can be obtained, for instance, by adapting the methodology of Stein presented in Gast (2017) and Ying (2016). The main technical difficulty would be to adapt the perturbation theory from these papers to the case of the solution of the graphon mean field approximation.