1

arXiv:1510.00761v1 [cs.PF] 3 Oct 2015

On the Rate of Convergence of Mean-Field
Models: Stein’s Method Meets the
Perturbation Theory
Lei Ying
School of Electrical, Computer and Energy Engineering
Arizona State University, Tempe, AZ 85287
lei.ying.2@asu.edu
Abstract
This paper studies the rate of convergence of a family of continuous-time Markov chains (CTMC)
to a mean-field model. When the mean-field model is a finite-dimensional dynamical system with a
unique equilibrium point, an analysis based on Stein’s method and the perturbation theory shows
that under some mild conditions, the stationary distributions of CTMCs converge (in the mean-square
sense) to the equilibrium point of the mean-field model if the mean-field model is globally asymptotically
stable and locally exponentially stable. In particular, the mean square difference between the Mth CTMC
in the steady state and the equilibrium point of the mean-field system is O(1/M ), where M is the size
of the Mth CTMC. This approach based on Stein’s method provides a new framework for studying
the convergence of CTMCs to their mean-field limit by mainly looking into the stability of the meanfield model, which is a deterministic system and is often easier to analyze than the CTMCs. More
importantly, this approach quantifies the rate of convergence, which reveals the approximation error
of using mean-field models for approximating finite-size systems.

I. Introduction
The mean-field method is to study large-scale and complex stochastic systems using simple
deterministic models. The idea of the mean-field method is to assume the states of nodes
in a large-scale system are independently and identically distributed (i.i.d.). Based on this
i.i.d. assumption, in a large-scale system, the interaction of a node to the rest of the system
can be replaced with an “average” interaction, and the evolution of the system can then be
modeled as a deterministic dynamical system, called a mean-field model. Then the macroscopic
behaviors of the stochastic system can be approximated using the mean-field model, e.g., the
stationary distribution of the stochastic system may be approximated using the equilibrium
point of the mean-field model. The mean-field method has important applications in various
areas including statistical physics, epidemiology, queueing theory, and game theory. Here are
just a few examples of these applications [1]–[6].
This paper focuses on the use of mean-field models for stationary distribution approximation
of CTMCs, which has been used in queueing networks and recently in cloud computing
systems for quantifying the performance of large-scale communication and computing systems such as data centers. Besides solving the equilibrium point of the mean-field model, a
critical component of the mean-field method in this application is to prove that the stationary
distributions of a family of CTMCs indeed converge to the equilibrium point of the mean-field
model, which justifies the use of the mean-field model.

2

In this paper, we consider a family of CTMCs, where the Mth CTMC is an M-dimensional
continuous-time Markov chain W( M) ∈ U M , where the superscript M is the number of nodes
(or called particles) in the system and U M ⊆ Rn is the state space of the CTMC. We assume U is
a finite state space and the CTMC is irreducible. Without loss of generality, let U = {1, · · · , n}.
We further define
1 M
( M)
xi ( t ) =
1 ( M)
Wm (t)=i
M m∑
=1
( M)

(t) is n
the fraction ofonodes in state i at time t.
This paper focuses on the case such that x( M) = x( M) (t), t ≥ 0 is also an (n-dimensional)
CTMC, i.e., the CTMC is a population process [7], [8]. We remark that many applications of the
mean-field method such as those in queueing networks and epidemiology are for population
processes.
Now let x( M) (∞) denote the stationary distribution of the Mth CTMC. Furthermore, let x(t)
denote the solution of an associated mean-field model and x∗ denote its equilibrium point.
Existing approaches for proving the convergence of x( M) (∞) to x∗ often involve the following
three components.
(1) The first component is to show the convergence of CTMCs to the trajectory of the meanfield model for any finite time interval [0, t], i.e.,

where 1 is the indicator function, so xi

lim sup d(x( M) (s), x(s)) = 0,

M → ∞ 0≤ s ≤ t

(1)

where d(·, ·) is some measure of distance. This can be proved using different techniques
including Kurtz’s theorem [5], [7]–[9], propagation of chaos [10]–[12], or the convergence
of the transition semigroup of CTMCs [4], [13].
(2) The second component is to prove the asymptotic stability of the mean-field model, i.e.,
lim x(t) = x∗ .

t→∞

Lyapunov theorem or LaSalle invariance principle can often be used for proving the
stability.
(3) After establishing the previous two results, we obtained
lim lim x( M) (t) = lim x(t) = x∗ .

t→∞ M→∞

t→∞

The convergence of the stationary distributions can then be concluded if we can prove the
interchange of the limits, i.e.,
lim x( M) (∞) = lim lim x( M) (t) =(a) lim lim x( M) (t) = x∗ ,

M→∞

M→∞ t→∞

t→∞ M→∞

where step (a) is called the interchange of the limits.
Since these approaches are all based on the interchange of the limits and use the finite-time
convergence (equality (1)) as the stepping stone, they are indirect methods of proving
lim x( M) (∞) = x∗ .

M→∞

3

Because of this reason, it is difficult to use these approaches to establish the rate of convergence
of stationary distributions, and to provide bounds on the approximation error kx( M) (∞) − x∗ k
for the finite-size system (i.e., for any fixed M).
This paper directly studies the convergence of the stationary distributions of CTMCs using
Stein’s method [14], [15], which is a method to bound the distance of two probability distributions. Our use of Stein’s method for the rate of convergence was inspired by the work by
Braverman and Dai [16], in which they developed a modular framework with three components
for steady-state diffusion approximations and established the rate of convergence to diffusion
models for M/Ph/n + M queuing systems. The results in this paper also share similar spirit
with the work by Gurvich [17], which establishes the rate of convergence of diffusion models
for steady-state approximations for exponentially ergodic Markovian queues.
This paper focuses on mean-field models instead of diffusion models. Based on Stein’s
method, the paper identifies a fundamental connection between the perturbation theory for
nonlinear systems and the convergence of mean-field models. The perturbation theory shows
that for a stable nonlinear system with exponentially stable equilibrium point, the error of
the first-order approximation of the nonlinear system is at the order of O(ǫ2 ), where ǫ is
the scaling factor of the perturbation. It turns out the mean-square difference between the
stationary distribution of the Mth CTMC and the equilibrium point of the mean-field model
is related to the cumulative error (integrated over infinite time horizon) of using the first-order
approximation of the mean-field model. After quantifying the cumulative error, the following
results are established for finite-dimensional mean-field models.
• If the mean-field model is perfect (see definition in Section II), globally asymptotically
stable and locally exponentially stable, then the stationary distributions of the CTMCs
converges in the mean-square sense to the equilibrium point of the mean-field model
with rate O(1/M) (Theorem 1), specifically,
h
i
1
E k x ( M ) ( ∞ ) − x ∗ k 2 = O ( ).
(2)
M
• If the mean-field model is not perfect, sufficient conditions that guarantee the convergence
of the stationary distributions have been obtained in Corollary 1.
We remark that these results are different from the celebrated law of large numbers for
Markov chains established by Kurtz [7], [8], where the convergence is established for sample
paths of the CTMCs over a finite time interval, not for the stationary distributions of the
CTMCs. The contributions of these results are two-fold: First, it provides a direct method
of studying the convergence of stationary distributions to its mean-field limit. The method
connects the convergence of CTMCs with the stability of the mean-field model. Note that
the mean-field model is a deterministic system, so it is often easier to analyze than the
CTMCs. Second, the method quantifies the rate of convergence, and provides bounds on the
approximation error when using the mean-field limit for approximating finite-size systems.
Finally, recall that
1 M
( M)
xi ( t ) =
1 ( M)
Wm (t)=i
M m∑
=1
which is the average of M Bernoulli random variables. The mean-square error in (2) is at the
same order as the variance of the average of M i.i.d. Beroulli random variables (according
to law of large numbers). While it is obvious that 1 ( M)
are not independent in the Mth
Wm

(t)=i

4

CTMC, the bound on the mean-square difference, however, provides an intuitive support to
approximating a large-scale CTMC based on the i.i.d. assumption and mean-field models.
II. Mean-Field Models
Consider an M-dimensional continuous-time Markov chain W( M) ∈ U M , where the superscript M is the number of nodes (or called particles) in the system and U M ⊆ R M is the state
space of the CTMC. We assume U is a finite state space and the CTMC is irreducible. Without
loss of generality, we assume U = {1, · · · , n}. We further define
( M)
Xi ( t )
( M)

where 1 is the indicator function, so Xi
further define

M

=

∑
m=1

1

( M)

Wm

(t)=i

(t) is the number of nodes in state i at time t. We

x( M ) (t) =

X( M ) (t)
,
M

( M)

(nt) ∈ [0, 1] represents
the fraction of nodes in state i at time t. In this paper, we assume
o
x( M) = x( M) (t), t ≥ 0 is an (n-dimensional) CTMC. We use x( M) (∞) to denote its stationary
distribution.
Furthermore, we have a mean-field model described by the following autonomous dynamical
system:
d
ẋ , x (t) = f (x(t)) x(0) = x and x(t) ∈ D ⊆ [0, 1]n ,
(3)
dt
where D is a compact set. Here, we abuse the notation and use x to denote the initial condition,
which simplifies the notation in the analysis later without causing too much confusion. Assume
the system has a unique equilibrium point and let x∗ denote the equilibrium point. The key
idea of the mean-field analysis is to use the solution of this deterministic dynamical system to
approximate the behavior of the CTMC when M is large; for example, use x∗ to approximate
x ( M ) ( ∞ ).
Let Qx( M),y( M) denote the transition rate of the CTMC from state x M to state y M . The sequence
of CTMCs is called a density-dependent family of CTMCs if the normalized transition rate
1
qx( M) ,y( M) = M
Qx( M) ,y( M) only depends on x( M) and y( M) but is independent of M (see a detailed
definition in [5]). For a density-dependent family of CTMCs, the mean-field model can often
be obtained by choosing
f (x) = ∑ qx,y (y − x)
so xi

y:y 6= x

because qx,y is the transition rate from x to y and y − x is the change of system state when
such a transition occurs.
We next illustrate the idea using an SIS (susceptible-infected-susceptible) model with external
infection source, which is a variation of the SIS model.
( M)
( M)
Example: Let Wm denote the state of an individual such that Wm = 0 if the individual
( M)
( M)
is susceptible and Wm = 1 if the individual is infected. So x0 is the fraction of susceptible
( M)
is the fraction of infected individuals. We assume the recover time of
individuals and x1
an individual follows an exponential distribution with mean 1. Each infected node randomly

5

selects an susceptible node after waiting for a random time that is exponentially distributed
with mean 1/β and infects it. Each susceptible node, after it becomes susceptible, gets infected
by an external infection source after a random time period that is exponentially distributed
with mean 1/α. Therefore, W ( M) , X(m) and x( M) are CTMCs. Specifically, x( M) has the following
transition rates:



−
1
( M) ( M)
( M)

1
(
M
)
(
M
)

if y
=x
+M
Mαx0 + Mβx0 x1 ,



 1 


1
( M)
1
if y( M) = x( M) + M
Mx1 ,
Qx( M),y( M) =
−
1



(
M
)
(
M
)
(
M
)
(
M
)
 − Mαx

− Mβx0 x1 − Mx1 , if y( M) = x( M)

0

0
otherwise.

Note for a given M, computing the stationary distribution of x( M) is not easy because it has
a large state space
2

1 2
, ,··· ,1
M M
and the transition rates are nonlinear functions of the state.
The SIS considered above is a density-dependent family, so we consider the following meanfield model






ẋ0
1
−1
(4)
= f (x) = ∑ qx,y (y − x) = (αx0 + βx0 x1 )
+ x1
−1
1
ẋ1
y:y 6= x

To solve the mean-field model above, we notice that x0 + x1 = 1 always holds, so we only need
to consider
ẋ0 = −αx0 − βx0 (1 − x0 ) + (1 − x0 ).
The equilibrium point can then be obtained by solving
0 = −αx0 − βx0 (1 − x0 ) + (1 − x0 ).
For example, if α = β = 0.5, then
x0∗ = 2 −

√

2 and

x1∗ =

√

2 − 1,

which can be used to approximate the fractions of susceptible and infected populations when
M is large, i.e., the stationary distribution of x( M) . The simulation results of the fraction of
infected population with M = 100, 1, 000 and 100, 000 are shown in Figure 1, from which the
√
( M)
convergence of x0 to 2 − 2 can be observed.


6

Normalized Infected Population (x0)

0.8
100 nodes
1,000 nodes
10,000 nodes

0.75
0.7
0.65
0.6
0.55
0.5
0.45
0.4
2000

2100

2200

2300

2400

2500

2600

2700

2800

2900

3000

Time

Fig. 1: Simulation results of the normalized infected population with M = 100, 1, 000, and
100, 000. α = β = 0.5 in these simulations. The CTMCs were simulated using the uniformization
method. The time is the discrete-time (scaled with M).
III. Stein’s Method for the Rate of Convergence
In this section, we study the convergence of the CTMCs to a mean-field model using Stein’s
methodqand the perturbation theory. Throughout this paper, k · k denotes the 2-norm, i.e.,
kxk = ∑i xi2 , and | · | denotes the absolute value. For two vectors a, b ∈ Rn , a · b is the dot
product. Furthermore, ▽ g(x) denotes the gradient of g(x), and ▽ xi (t, x) refers to differentiating
with respect to the location x, and ẋ is the derivative with respect to time.
Recall the mean-field model defined in (3):
ẋ = f (x(t))

x(0) = x and x(t) ∈ D ⊆ [0, 1]n .

The mean-field model is said to be globally asymptotically stable if given any initial condition
x(0) ∈ D and any ǫ > 0, there exists t(x(0), ǫ) such that

k x ( t ) − x ∗ k ≤ ǫ ∀ t ≥ t ( x (0 ), ǫ ).

The mean-field model is said to be locally exponentially stable if there exist positive constants ǫ,
α and κ such that starting from any initial condition kx(0)k ≤ ǫ,

kx(t) − x∗ k ≤ κ kx(0)k exp (−αt) .

Let g(x) be the solution to the Poisson equation

▽g(x) · ẋ = ▽g(x) · f (x) =
Then,
g(x) = −

Z ∞
0

n

∑ ( xi − xi∗ )2

(5)

i =1

∑ (xi (t, x) − xi∗ )

2

dt

i

when the integral is finite [18], [19], where x(t, x) is the trajectory of the dynamical system with
x as the initial condition. The integral is finite when the mean-field model is asymptotically
stable and locally exponentially stable. Note that − g(x) can be viewed as the cumulative

7

square-deviation of the system state from the equilibrium point when the initial condition is
x.
Now let Gx( M) denote the generator for the Mth CTMC, then
Gx ( M ) g ( x ) =

∑
y:y 6= x

= M

Qx,y (x) ( g(y) − g(x))

∑

y:y 6= x

qx,y (x) ( g(y) − g(x)) .

Since x( M) is irreducible and has finite state space, x( M) has a stationary distribution. Initializing x( M) (0) according to its stationary distribution, and using Ex( M) [·] throughout to denote
stationary expectation, we have
"
#


Ex( M) Gx( M) g(x) = Ex( M) M ∑ qx,y (x) ( g(y) − g(x)) = 0,
(6)
y:y 6= x

where the subscript in the expectation indicates that the expectation is taken over the stationary
distribution of x( M) . Then by taking expectation of the Poisson equation (5) over the stationary
distribution of x( M) and then adding (6) to the equation, we obtain
"
#
n

Ex ( M )

∑ ( xi − xi∗ )2

"i = 1

= Ex ( M ) ▽ g ( x ) · f ( x ) − M

∑
y:y 6= x

qx,y ( g(y) − g(x))

#
(7)

Now adding and subtracting ▽ g(x) · ∑y:y6=x qx,y M(y − x) yields
"
#
n

Ex ( M )

∑ ( xi − xi∗ )2

"i = 1

= Ex ( M ) ▽ g ( x ) · f ( x ) − ▽ g ( x ) ·
"

= Ex ( M ) ▽ g ( x ) ·

f (x) −

∑
y:y 6= x

∑
y:y 6= x

qx,y M(y − x) − M

qx,y M(y − x)

!

−

∑
y:y 6= x

∑
y:y 6= x

qx,y ( g(y) − g(x) − ▽ g(x) · (y − x))
#

qx,y M ( g(y) − g(x) − ▽ g(x) · (y − x)) . (8)

h
2 i
From the equality above, intuitively, that Ex( M) ∑in=1 xi − xi∗
converges to zero as M → ∞
can be established if the followings are true:
• Bounded gradient of g(x) : k▽ g(x)k is bounded by a constant independent of M.
• Convergence of the generator:
"
#
lim Ex( M)

M→∞

•

f (x) −

∑

y:y 6= x

h

#

qx,y M(y − x)
i

= 0.

Bounded transition-rate of the CTMC: Ex( M) ∑y:y6=x qx,y is bounded.

8

•

Diminishing first-order approximation error:

k g(y) − g(x) − ▽g(x) · (y − x)k = O



1
M2



.

Note that g(x) + ▽ g(x) · (y − x) is the first-order Taylor approximation of g(y).
For many CTMCs and the associated mean-field models, the first three conditions mentioned
above can be easily verified. In the following theorem, we will prove that the last condition
holds when the mean-field model is globally asymptotically stable and locally exponentially
stable (see inequality (13)), and then establish the rate of convergence based on that. The
following theorem presents the main result of this paper.
Theorem 1. The stationary distributions of the CTMCs (x( M) (∞)), defined in Section II, converge to
the equilibrium point (x∗ ) of the mean-field model (3) in the mean-square sense with rate 1/M, i.e.,
"
#
 
n
1
∗ 2
Ex ( M ) ∑ ( x i − x i ) = O
M
i =1
when the following conditions hold:
• Bounded transition-rate condition: There exists a constant c > 0 independent of M such that
#
"
Ex ( M )

•
•

∑

y:y 6= x

Bounded state transition condition: There exists a constant c̃ independent of M such that
kx − yk ≤ Mc̃ for any x and y such that qx,y 6= 0.
Perfect mean-field model condition: The mean-field model (3) is given by
f (x) =

∑
y:y 6= x

•
•

qx,y ≤ c.

qx,y M(y − x)

∀ x.

Partial derivative condition: The function f (x) is twice continuously differentiable.
Stability condition: The mean-field model is globally asymptotically stable and is locally exponentially stable.

Remark. The first four conditions are easy to verify, so only the stability condition requires
nontrivial work. Since a dynamical system has an exponentially stable equilibrium point if
and only if the linearized system (at the equilibrium) is exponentially stable (see Theorem 4.15
in [20]), the local exponential stability can be verified by calculating the eigenvalues of the
state matrix of the mean-field model. When the parameters of the mean-field model are given,
the local exponential stability can be numerically verified. The global asymptotical stability in
general is studied using the Lyapunov theorem.
Remark. It is worth to pointing out that if the mean-field model is unstable but the perfect
mean-field model assumption holds. Kurtz’s theorem indicates that the sample paths of the
CTMCs converge to the trajectory of the mean-field model for any finite time interval, which
implies that the CTMCs are “unstable” as well.
Proof: We first prove the theorem assuming the mean-field model is globally exponentially
stable, and then extend it to the general case. Under the perfect mean-field model assumption,

9

equation (8) becomes
Ex ( M )

"

n

∑ (xi − xi∗ )2

"i =1

= Ex ( M ) −

∑

y:y 6= x

#
#

qx,y M ( g(y) − g(x) − ▽ g(x) · (y − x)) ,

where
g(x) = −

Z ∞
0

∑ ( xi (t, x) − xi∗ )

2

dt.

i

We next focus on the following term,

− ( g(y) − g(x) − ▽g(x) · (y − x))
(9)
Z ∞


2
2
=
∑ ( xi (t, y) − xi∗ ) − ( xi (t, x) − xi∗ ) − 2 ( xi (t, x) − xi∗ ) ▽xi (t, x) · (y − x) dt (10)
0

i

Note that we exchanged the order of integration and differentiation for the third term. This is
can be done because
Z ∞
2 ( xi (t, x) − xi∗ ) ▽ xi (t, x) · (y − x) dt
0


is finite, which can be proved using the fact that both xi (t, x) − xi∗ and k▽ xi (t, x)k decay
exponentially fast to zero as t increases (apply inequalities (22) and (33) with z = 1), and the
fact that ky − xk is bounded due to the bounded state transition condition.
We next define
ei (t) = xi (t, y) − xi (t, x) − ▽ xi (t, x) · (y − x),
i.e.,

xi (t, y) = ei (t) + xi (t, x) + ▽ xi (t, x) · (y − x),

so

( xi (t, y) − xi∗ )2 − ( xi (t, x) − xi∗ )2 − 2 ( xi (t, x) − xi∗ ) ▽xi (t, x) · (y − x)

= (ei (t) + xi (t, x) − xi∗ + ▽xi (t, x) · (y − x))2 − ( xi (t, x) − xi∗ )2 − 2 ( xi (t, x) − xi∗ ) ▽xi (t, x) · (y − x)

= e2i (t) + (▽xi (t, x) · (y − x))2 + 2ei (t)▽xi (t, x) · (y − x) + 2ei (t) ( xi (t, x) − xi∗ )

= ei (t) ( ei (t) + 2▽xi (t, x) · (y − x) + 2 ( xi (t, x) − xi∗ )) + (▽xi (t, x) · (y − x))2 .

According to the perturbation theory, in particular, inequality (39), when the system is
exponentially stable, we have that


1
.
|ei (t)| ≤ ke(t)k = O
M2
According to the bounded state transition condition,

k x − yk ≤

c̃
.
M

10

Furthermore, both ▽ xi (t, x) and xi (t, x) are bounded (see inequalities (22) and (33)) by constants
independent of M and t. Therefore, we can choose a constant b and a sufficiently large M̃ such
that for any M ≥ M̃,

|ei (t) + 2▽xi (t, x) · (y − x) + 2 ( xi (t, x) − xi∗ )| ≤ b,

which implies that

| g(y) − g(x) − ▽g(x) · (y − x)| ≤ b

Z ∞
0

√

≤ b n

∑ |ei (t)| dt +

i
Z ∞
0

Z ∞
0

ke(t)k dt +

∑ (▽xi (t, x) · (y − x))2 dt,

Z ∞
0

(11)

i

∑ (▽xi (t, x) · (y − x))

2

dt, (12)

i

where the last
√ inequality is based on the following relation between 1-norm and 2-norm:
∑i |ei (t)| ≤ nke(t)k. In Section IV, we will show that under the exponential stability assumption,
Z
∞

0

ke(t)k dt = O(1/M2 ) see inequality (40).

From the bounded state transition condition, ky − xk2 ≤
Z ∞
0

c̃2
x
(
t,
x
)
·
(
y
−
x
))
dt
≤
(▽
∑ i
M2
i
2

Z ∞
0

c̃2
.
M2

Therefore,

∑ k▽xi (t, x)k

2

dt.

i

Now according to inequality (33), there exist positive constants b1 and b2 , both independent
M, such that
Z ∞
0

2
∑ (▽xi (t, x) · (y − x)) dt ≤
i

c̃2
M2

Z ∞
0

b1 exp(−b2 t) dt ≤

b1 1
.
b2 M2

Therefore, we can conclude that

| g(y) − g(x) − ▽g(x) · (y − x)| = O



1
M2



,

(13)

which implies that
Ex ( M )

"

n

∑ ( xi − xi∗ )2

i =1

#

=O



1
M



Ex ( M )

"

∑

#

qx,y .

y:y 6= x

Finally, using the bounded transition rate condition, we conclude
"
#
 
n
1
∗ 2
.
Ex ( M ) ∑ ( x i − x i ) = O
M
i =1

(14)

(15)

Now consider the case that the mean-field model is not globally exponentially stable, but
is globally asymptotically stable and locally exponentially stable. Recall that D ⊆ [0, 1]n is
compact. According to the definition of global asymptotic stability (Definition 4.4 in [20]),
given any ǫ > 0, there exists a finite time t′ such that

kx(t)k ≤ ǫ

11

for any t ≥ t′ . For any finite t, following a similar analysis as in Section IV (or Section 10.1
in [20]), ke(t, x)k = O(1/M2 ) holds1 Therefore, we can bound the term (10) by separating the
integration into two intervals: from 0 to t′ , and from t′ to ∞, where t′ is chosen such that
x(t) converges exponentially to the equilibrium point after t′ . Since ke(t′ , x)k = O(1/M2 ), the
analysis above applies to the integration over [t′ , ∞). Hence, the result holds.
The theorem above requires a perfect mean-field model and bounded state transitions. Both
conditions can be relaxed, but the rate of convergence will be different.
Corollary 1. Assume partial derivative condition and the stability condition in Theorem 1 hold. The
stationary distributions of the CTMCs converge (in the mean square sense) to equilibrium point of the
mean-field model, i.e.,
"
#
n

∑ ( xi − xi∗ )2

lim Ex( M)

M→∞

when the following conditions also hold:
"

f (x) −

lim Ex( M)

M→∞

lim Ex( M)

M→∞

=0

i =1

∑
y:y 6= x

"

= 0.

(16)

qx,y Mky − xk2 = 0

(17)

max ky − xk = 0

(18)

qx,y M(y − x)

∑
y:y 6= x

lim

#
#

M →∞ x,y:q x,y >0

We say that the mean-field model is asymptotically accurate when condition (16) holds, which replaces
the perfect mean-field model condition. Conditions (17) and (18) replace the bounded state transition
condition.
Proof: First recall that we have
"
#
n

Ex ( M )

∑ ( xi − xi∗ )2

"i = 1

= Ex ( M ) ▽ g ( x ) ·
≤



f (x) −


max k▽ g(x)k Ex( M)
x

Ex ( M )

"

∑
y:y 6= x

"

∑
y:y 6= x

qx,y M(y − x)

f (x) −

∑
y:y 6= x

!

−

∑
y:y 6= x

qx,y M(y − x)
#

qx,y M | g(y) − g(x) − ▽ g(x) · (y − x)| .

qx,y M ( g(y) − g(x) − ▽ g(x) · (y − x))
#

+

#

(19)
(20)

By choosing z = 1 in Section IV, it is easy to verify according to inequality (33) that
(maxx k▽g(x)k) is upper bounded by a constant independent of M. Therefore, under condition
(16), (19) → 0 as M → ∞.
1 This holds without exponential stability, but the constant in O (1/M2 ) may be a function of t if the system is not
exponentially stable.

12

A careful examination of inequality (40) shows that
Z ∞


2
′
ke(t)k dt = O ky − xk exp α3 ky − xk .
0

So under condition (18), we have

Z ∞
0



ke(t)k dt = O ky − xk2 .

When condition (18) holds, following the analysis that leads to inequality (12), we can again
show that there exists a constant b̃ independent M such that
Z ∞
√ Z ∞
2
ke(t)k dt +
| g(y) − g(x) − ▽g(x) · (y − x)| ≤ b̃ n
∑ (▽xi (t, x) · (y − x)) dt.
0

0

i

According to inequality (33), we also have
Z ∞
0

Therefore, we have



2
2
.
x
(
t,
x
)
·
(
y
−
x
))
dt
=
O
k
y
−
x
k
(▽
∑ i
i

(20) = O Ex( M)

"

∑
y:y 6= x

qx,y Mky − xk2

#!

,

which converges to zero according to condition (17). Hence, the corollary holds.
Remark. When the mean-field model is asymptotically accurate, the convergence rate depends
on the convergence rates of (16) and (17).
Example: Let us first go back to the SIS model introduced in Section II. A closed-form
solution can be obtained for x0 (t). Again assume α = β = 0.5, then the solution of the ordinary
differential equation is
√
√ 
√ 
√
e− 2t (2 + 2) x0 (0) − 2 + 2 + (2 − 2) x0 (0) − 2
x0 ( t ) =
,
√ 
√ 
√
−
2t
x 0 (0 ) − 2 + 2 + x 0 (0 ) − 2 − 2
−e
√
which converges to 2 − 2 as t → ∞ independent of x0 (0). Therefore, it is easy to verify
that the system is globally, asymptotically stable. Furthermore, the linearized system at the
equilibrium is
ǫ̇0 = − (α + β(1 − x0∗ ) + 1) ǫ0 ,

where ǫ0 = x0 − x0∗ and x0∗ is the equilibrium value, so the equilibrium point is locally
exponentially stable. Furthermore, the mean-field model is perfect in this case and it can be
easily verified that all other conditions in Theorem√1 hold. So in √
the mean square sense,
stationary distributions convergeh to the x0i = 2 − 2 and x1 =
2 − 1 with rate 1/M.
2
Numerical evaluation of MEx( M) ( x0 − x0∗ ) versus M is shown in Figure 2, where M varies
h
2 i
from 100 to 1,000. We can see that MEx( M) ∑in=1 xi − xi∗
varies within the interval [0.21,
0.27] while the size of the system
increases by 10 times √
(from 100 to 1,000). The standard
√
deviation (deviation
from 2 − 2) is 0.02177 (3.72%×(2 − 2))) when M = 100 and is 0.0068
√
(1.16%×(2 − 2)) when M = 1, 000.

M× Mean Square Difference

13

0.28

0.26

0.24

0.22

0.2
100

200

300

400

500
600
System size (M)

700

800

900

1000

i
h
Fig. 2: Numerical evaluation of MEx( M) ( x0 − x0∗ )2 versus M.
IV. The Perturbation Theory
In this section, we summarize the results of the perturbation theory for nonlinear systems.
These results are special cases of those results presented in [20] because we only need to
consider a perturbation to the initial condition. Furthermore, the mean-field model considered
in this paper is an autonomous system, which again is a special case of the nonlinear system
considered in [20]. For these reasons, the analysis of the perturbation results can be simplified.
On the other hand, the perturbation method introduced in [20] only states that the 2-norm of
the following error is is at the order of ky − xk2 independent of t (under certain conditions)
e(t) = x(t, y) − x(t, x) − ▽x(t, x) · (y − x)

Our result on the rate of convergence, however, requires such an upper bound on the cumulative error, i.e., an upper bound on
Z
∞

0

e(t) dt.

Therefore, it is necessary to go through the detailed analysis for the system considered in this
paper to establish the result for the cumulative error. For the completeness of the paper and
the easy reference of the reader, we next introduce the perturbation results in [20] with a more
detailed calculation of ke(t)k, which shows that not only the approximation error is bounded,
but the upper bound decays exponentially to zero as t increases. The analysis closely follows
[20].
Consider the system
ẋ = f (x)
(21)
where f : D ⊆ [0, 1]n → Rn . Without the loss of generality, we assume x∗ = 0. We are interested
in comparing the solution of this nominal system with the system with a perturbation on the
initial condition x(0) = x + ǫz, where z = 1ǫ (y − x) and is an n-dimensional vector. For the
mean-field analysis considered in this paper, ǫ = 1/M. Under the condition of Theorem 1, for
any neighboring states x and y,
1
kzk = k(y − x)k ≤ c̃.
ǫ

14

Let x(t, ǫ) to denote the solution of the dynamical system with initial perturbation ǫ. Note
that the dependence of the solution on y − x is omitted to simplify the notation. The analysis
holds for any y and x. We next first repeat the assumptions on the nominal dynamical system.
Assumption 1. For any i, the function f i (x) is twice continuously differentiable. Therefore, the
∂f
Jacobian matrix of f (x), denoted by ∂x , is Lipschitz. In other words, there exists a constant
L > 0 such that
∂f
∂f
(x) −
( y ) ≤ L k x − y k.
∂x
∂x


Assumption 2. The dynamical system (21) has a unique equilibrium point and is exponentially
stable. In other words, there exist positive constants α and κ such that starting from any initial
condition x(0) ∈ D ,
kx(t)k ≤ κ kx(0)k exp (−αt) .
(22)
Under this assumption, according to Theorem 4.14 in [20], there exist a Lyapunov function V (x)
and positive constants cu , cl , cd , and c p such that for any x ∈ D , the following inequalities hold
c l k x k2 ≤ V ( x ) ≤ c u k x k2
V̇ (x) ≤ −cd kxk2
k▽V (x)k ≤ c p kxk.



We first consider the finite Taylor series for x(t, ǫ) in terms of ǫ :
x(t, ǫ) = x(0) (t) + ǫx(1) (t) + e(t),

(23)

x(0, ǫ) = x + ǫz,

(24)

and
where
x(0) (t) = x(t, 0) and x(1) (t) =

dx
(t, ǫ)
dǫ

.
ǫ =0

Substituting (23) into the dynamical system equation, we get
ẋ (t, ǫ) = ẋ(0) (t) + ǫẋ(1) (t) + ė(t) = f (x(t, ǫ))

= h




(0)

(x

(0)

(t)) + h

(1)

(x

(≤1)

(t))ǫ + Re (t, ǫ),

where x(≤1) = x(0) , x(1) . The zero-order term h(0) is given by




with x(0) (0) = x,
ẋ (0) (t) = h(0) x(0) (t) = f x(0) (t)

(25)
(26)

15

which is the nominal system without the perturbation on the initial condition. The first-order
term is given by


d
(≤1)
(1)
f (x(t, ǫ))
x
(t) =
h
dǫ
ǫ =0
∂f
dx
=
(x(t, ǫ)) (t, ǫ)
∂x
dǫ
ǫ =0
∂ f (0)
(x (t))x(1) (t).
=
∂x
Recall that

∂f
∂x

is the Jacobian matrix. Therefore, we have
ẋ(1) (t) =

∂ f (0)
(x (t))x(1) (t)
∂x

with

x(1) (0) = z.

(27)

We next study e(t) = x(t, ǫ) − x(0) (t) − ǫx(1) (t). Combining the results above, we have


∂f
ė (t) = f (x(t, ǫ)) − f x(0) (t) − ǫ (x(0) (t))x(1) (t)
∂x
e(0) = 0.
Now by defining


ρ x(≤1) (t), e(t), ǫ
= f


γ x(≤1) (t), ǫ
= f

we obtain



 ∂f
e(t) + x(0) (t) + ǫx(1) (t) − f x(0) (t) + ǫx(1) (t) − (x(0) (t))e(t)
∂x




∂
f
x(0) (t) + ǫx(1) (t) − f x(0) (t) − ǫ (x(0) (t))x(1) (t),
∂x






∂ f (0)
(x (t))e(t) + ρ x(≤1) (t), e(t), ǫ + γ x(≤1) (t), ǫ .
∂x
Note that both ρ and γ are n-dimensional vectors. It is easy to see that


(≤1)
ρ x
(t), 0, ǫ = 0.
ė (t) =

(28)

(29)

According to Taylor’s theorem and the mean value theorem, we have




∂2 f l
(1)
(1)
γl x(≤1) (t), ǫ
= ǫ2 x(1) (t)T H( f l )(ξ )x(1) (t) = ǫ2 ∑
( ξ ) xi ( t ) x j ( t )
∂xi ∂x j
i,j
for ξ = x(0) (t) + αǫx(1) (t) for some 0 ≤ α ≤ 1. H( f l ) is the Hessian matrix of function f l (x).
Then we have




∂2 f l
(1)
(1)
( ξ ) xi ( t ) x j ( t ) .
γl x(≤1) (t), ǫ
= ǫ2 ∑
∂xi ∂x j
i,j
Furthermore, we have

∂ fl
∂f
∂ρl  (≤1)
x
( t ), e ( t ), ǫ =
(e(t) + x(0) (t) + ǫx(1) (t)) − l (x(0) (t)).
∂ei
∂xi
∂xi

16

According to the mean-value theorem and (29), we have that


 ∂f
∂ f (0)
(0)
(1)
(≤1)
(ẽ (t) + x (t) + ǫx (t)) − (x (t)) e(t),
ρ x
( t ), e ( t ), ǫ =
∂x
∂x

where ẽ(t) = ae(t) for some 0 ≤ a ≤ 1. According to the Lipschitz condition in Assumption
(1) and the Cauchy - Schwarz inequality, we have




ρ x(≤1) (t), e(t), ǫ
≤ L ke(t)k + ǫkx(1) (t)k ke(t)k.

Now we utilize the assumption that the nominal system (21) converges to the equilibrium
point exponentially fast from any initial condition in the domain. We use the Lyapunov function
in Assumption (2) to bound ke(t)k. We start from

=
=
≤( a)
=

≤

V̇ (e(t))
▽V (e(t)) · ė (t)
▽V (e(t)) · f (e(t)) + ▽V (e(t)) · (ė(t) − f (e(t)))
−cd V (e(t)) + ▽V (e(t)) · (ė (t) − f (e(t)))

∂ f (0)
∂f
∂f
−cd V (e(t)) + ▽V (e(t)) ·
(x (t))e(t) − (0)e(t) + (0)e(t) − f (e(t))
∂x
∂x
∂x




+ρ x(≤1) (t), e(t), ǫ + γ x(≤1) (t), ǫ

∂f
∂f
∂ f (0)
(x (t))e(t) − (0)e(t) +
(0)e(t) − f (e(t))
−cd V (e(t)) + k▽V (e(t))k
∂x
∂x
∂x



 
(≤1)
(≤1)
+ ρ x
( t ), e ( t ), ǫ + γ x
( t ), ǫ

where inequality (a) is due to assumption (2) and the last inequality is a result of the Cauchy
Schwarz inequality. Note that based on Assumption (1) and the mean-value theorem, we have
∂ f (0)
∂f
(x (t))e(t) − (0)e(t)
∂x
∂x
∂f
(0)e(t) − f (e(t))
∂x

≤ Lkx(0) (t)kke(t)k
≤ Lke(t)k2 .

We also know that

where we define



ρ x(≤1) (t), e(t), ǫ


γ x(≤1) (t), ǫ
v
u
u
A(t) = t ∑
l

to simplify the notation.



≤ L ke(t)k + ǫkx(1) (t)k ke(t)k

= ǫ2 A ( t ),



∂2 f l
(1)
(1)
(
ξ
)
x
(
t
)
x
(
t
)
∑ ∂xi ∂x j
i
j
i,j

!2

17

Summarizing the results above, we get
V̇ (e(t))


(0)



(1)

≤ −cd V (e(t)) + L k▽V (e(t))k kx (t)k + ǫkx (t)k + 2ke(t)k ke(t)k + k▽V (e(t))k A(t)ǫ2


(0)
(1)
≤ −cd V (e(t)) + Lc p kx (t)k + ǫkx (t)k + 2ke(t)k ke(t)k2 + c p A(t)ǫ2 ke(t)k
q

c p  (0)
cp
≤ −cd V (e(t)) + L
kx (t)k + ǫkx(1) (t)k + 2ke(t)k V (e(t)) + √ A(t)ǫ2 V (e(t)) .
cl
cl
p
Define W (t) = V (t), then we have

cp
c
L c p  (0)
Ẇ (e(t)) ≤ − d W (e(t)) +
kx (t)k + ǫkx(1) (t)k + 2ke(t)k W (e(t)) + √ A(t)ǫ2 .
2
2 cl
c1
By the comparison lemma in [20], we have

t
cp
W (t) ≤ φ(t, 0)W (0) + √ ǫ2
φ(t, τ ) A(τ ) dτ
c1
0
Z
cp 2 t
φ(t, τ ) A(τ ) dτ
= √ ǫ
c1
0

Z

(30)
(31)

where the transition function φ(t, τ ) is


Z

L c p t  (0)
cd
(1)
kx (γ)k + ǫkx (γ)k + 2ke(γ)k dγ ,
φ(t, τ ) = exp − (t − τ ) +
2
2 cl τ

and the equality holds because e(0) = 0.
The following lemma proves that the first-order system (27) converges exponentially starting
from any initial condition in D .

Lemma 1. The first-order system (27) is exponentially stable for any solution x(0) (t) that starts from
D.

Proof: That the nominal system is exponentially stable implies that the following linear
system
∂f
(0 ) x
ẋ =
∂x
∂f

is (globally) exponentially stable, and ∂x (0) is Hurwitz (Corollary 4.3 in [20]), which further
implies that there exists a positive definite symmetric matrix P such that
V (x) = x T Px
is a Lyapunov function for the linear system such that
V̇ (x) ≤ −kxk2 .

(32)

18

We start from
V̇ (x(1) (t))

=
=

▽V (x(1) (t)) · ẋ(1) (t)


∂ f (0)
∂f
∂f
(1)
(1)
(1)
(1)
(1)
(x (t))x (t) − (0)x (t)
▽V (x (t)) · (0)x (t) + ▽V (x (t)) ·
∂x
∂x
∂x

≤ ( a ) − x (1) ( t )
≤
≤

2

+ 2λmax (P) L x(0) (t)

x (1) ( t )

2



2λmax (P)
L x (0) ( t ) V x (1) ( t )
λmax (P)
λmin (P)
 


1
2λmax (P)
(0)
V x (1) ( t )
−
−
L x (t)
λmax (P)
λmin (P)

−

1

V (x(1) (t)) +

where inequality (a) is based on (32) and the definition of V (x), Assumption (1) and the
mean-value theorem, and λmax (P) is the largest eigenvalue of matrix P.
By the comparison lemma, we have


Z t
1
2λmax (P)
(0)
V (t) ≤ exp −
x (τ ) dτ V (0)
t+
L
λmax (P)
λmin (P)
0




2λmax (P) Lκ kx(0)k
1
≤ exp
exp −
t V (0 ),
αλmin (P)
λmax (P)
where the last inequality holds because the exponential convergence assumption (2) yields
Z t
0

x

(0)

(τ ) dτ ≤

Z ∞
0

x(0) (τ ) dτ ≤

κ kx(0)k
.
α

Recall that x(1) (0) = z, so V (0) = z T Pz = p and




1
2λmax (P) Lκ kx(0)k
p
(1)
2
exp −
exp
t .
kx (t)k ≤
λmin (P)
αλmin (P)
λmax (P)

(33)

From the lemma above and assumption (1), we have that there exists a constant µ such that




µp
1
2λmax (P) Lκ kx(0)k
(1)
2
A(t) ≤ µkx (t)k ≤
exp −
exp
t .
λmin (P)
αλmin (P)
λmax (P)
Consider the set such that ke(t)k ≤

cd cl
4c p L ,

we have

φ(t, τ )
(34)


Z t

c
L cp
kx(0) (γ)k + ǫkx(1) (γ)k dγ
(35)
≤ exp − d (t − τ ) +
2
2 cl τ
!!

√
λmax (P) p
L c p κ kx(0)k
λmax (P) Lκ kx(0)k
cd
,(36)
+ 2ǫ p
exp
≤ exp − (t − τ ) +
4
2 cl
α
αλmin (P)
λmin (P)

where last inequality yields from (22) and (33).
Recall we have inequality (31)
cp
W ( t ) ≤ √ ǫ2
c1

Z t
0

φ(t, τ ) A(τ ) dτ.

(37)

19

Substituting the bounds on φ(t, τ ) and A(τ ), we obtain


Z t
1
cd
2√
τ dτ,
W (t) ≤ ǫ c1 Z(kx(0)k)
exp − (t − τ ) −
4
λmax (P)
0
where
Z(kx(0)k)

!!

√
λ
(
P
)
p
c
λmax (P) Lκ kx(0)k
L p κ kx(0)k
max
×
+ 2ǫ p
exp
2 cl
α
αλmin (P)
λmin (P)


2λmax (P) Lκ kx(0)k
µp
exp
λmin (P)
αλmin (P)



!
c p µp
λmax (P) Lκ kx(0)k
2λmax (P)
√ c p L λmax (P)
κL c p
p
=
exp
+
kx(0)k + ǫ p
.
exp
c1 λmin (P)
α 2cl
λmin (P)
cl
αλmin (P)
λmin (P)
cp
=
exp
c1

In other words, we have

Z t




cd
1
ke(t)k ≤ ǫ Z(kx(0)k)
exp − (t − τ ) −
τ dτ
4
λmax (P)
0





cd 
1
1
 ǫ2 Z(kx(0)k) c
exp − λ (P) t − exp − 4 t ,
1
d
max
=
4 − λmax ( P ) 
 ǫ2 Z(kx(0)k)t exp − cd t ,
4
2

(38)
if

cd
4

6=

1
λmax (P )

(39)
.

otherwise.

Then

Z ∞
0


 ǫ2 Z(kx(0)k) 4λmax (P)
cd
ke(t)k dt ≤
 ǫ2 Z(kx(0)k) 16
,
c2

cd
4

6= λmax1 (P)
.
otherwise.

if

d

= ǫ2 Z(kx(0)k)

4λmax (P)
.
cd

It is easy to see that with properly defined α1 , α2 , α3 and α4 , we have



Z ∞
kx(0)k
√
kx(0)k
2
.
ke(t)k dt ≤ ǫ pα1 exp α2
+ ǫ pα3 exp α4
α
α
0

(40)

We keep the terms kx(0)k and α to show that the cumulative error depends on the initial condition and the convergence rate of the mean-field model. Furthermore, p = z T Pz ≤ λmax (P)kzk2 .
V. Conclusion
This paper studies the convergence of the stationary distributions of a family of CTMCs to
the mean-field limit. When the mean-field model is perfect, the rate of convergence (the meansquare difference) has been proved to be O(1/M). Based on Stein’s method for bounding
the distance of probability distributions and the perturbation theory for nonlinear systems, a
fundamental connection between the convergence to the mean-field limit and the stability of
the mean-field model has been established.

20

Acknowledgement
The author is very grateful to Jim Dai and Anton Braverman. Jim’s seminar on Stein’s method
for the steady-state diffusion approximations inspired this work. The discussions with Jim and
Anton had continuously stimulated the author during the writing of this paper.
References
[1] L. P. Kadanoff, “More is the same; phase transitions and mean field theories,” Journal of Statistical Physics, vol. 137, no.
5-6, pp. 777–797, 2009.
[2] N. T. J. Bailey, The mathematical theory of infectious diseases and its applications. Hafner Press, 1975.
[3] F. Baccelli, F. Karpelevich, M. Y. Kelbert, A. Puhalskii, A. Rybko, and Y. M. Suhov, “A mean-field limit for a class of
queueing networks,” Journal of statistical physics, vol. 66, no. 3-4, pp. 803–825, 1992.
[4] N. D. Vvedenskaya, R. L. Dobrushin, and F. I. Karpelevich, “Queueing system with selection of the shortest of two
queues: An asymptotic approach,” Problemy Peredachi Informatsii, vol. 32, no. 1, pp. 20–34, 1996.
[5] M. Mitzenmacher, “The power of two choices in randomized load balancing,” Ph.D. dissertation, University of California
at Berkeley, 1996.
[6] J.-M. Lasry and P.-L. Lions, “Mean field games,” Japanese Journal of Mathematics, vol. 2, no. 1, pp. 229–260, 2007.
[7] T. G. Kurtz, “Limit theorems for sequences of jump markov processes approximating ordinary differential processes,” J.
Appl. Probab., vol. 8, no. 2, pp. 344–356, 1971.
[8] ——, Approximation of population processes. SIAM, 1981, vol. 36.
[9] L. Ying, R. Srikant, and X. Kang, “The power of slightly more than one sample in randomized load balancing,” in Proc.
IEEE Int. Conf. Computer Communications (INFOCOM), Hong Kong, 2015.
[10] A.-S. Sznitman, “Topics in propagation of chaos,” in Ecole d’été de probabilités de Saint-Flour XIX1989, 1991, pp. 165–251.
[11] V. Anantharam and M. Benchekroun, “A technique for computing sojourn times in large networks of interacting queues,”
Probability in the engineering and informational sciences, vol. 7, no. 04, pp. 441–464, 1993.
[12] M. Bramson, Y. Lu, and B. Prabhakar, “Asymptotic independence of queues under randomized load balancing,” Queueing
Systems, vol. 71, no. 3, pp. 247–292, 2012.
[13] A. Mukhopadhyay and R. R. Mazumdar, “Analysis of load balancing in large heterogeneous processor sharing systems,”
arXiv preprint arXiv:1311.5806, 2013.
[14] C. Stein, “A bound for the error in the normal approximation to the distribution of a sum of dependent random variables,”
in Proc. Sixth Berkeley Symp. Math. Stat. Prob., 1972, pp. 583–602.
[15] ——, “Approximate computation of expectations,” Lecture Notes-Monograph Series, vol. 7, pp. i–164, 1986.
[16] A. Braverman and J. Dai, “Stein’s method for steady-state diffusion approximations of m/ph/n + m systems,” arXiv
preprint arXiv:1503.00774, 2015.
[17] I. Gurvich et al., “Diffusion models and steady-state approximations for exponentially ergodic markovian queues,” Adv.
in Appl. Probab., vol. 24, no. 6, pp. 2527–2559, 2014.
[18] A. D. Barbour, “Stein’s method and Poisson process convergence,” J. Appl. Probab., pp. 175–184, 1988.
[19] F. Gotze, “On the rate of convergence in the multivariate clt,” Ann. Probab., pp. 724–739, 1991.
[20] H. K. Khalil, Nonlinear systems. Prentice Hall, 2001.

