arXiv:1007.5030v1 [math.PR] 28 Jul 2010

Analysis of a Splitting Estimator for Rare Event
Probabilities in Jackson Networks
Jose Blanchet, Kevin Leder and Yixi Shi
February 16, 2018
Abstract
We consider a standard splitting algorithm for the rare-event simulation of overflow
probabilities in any subset of stations in a Jackson network at level n, starting at a fixed
initial position. It was shown in [8] that a subsolution to the Isaacs equation guarantees
that a subexponential number of function evaluations (in n) suffice to estimate such
overflow probabilities
within a given relative accuracy. Our analysis here shows that in

fact O n2β+1 function evaluations suffice to achieve a given relative precision, where
β is the number of bottleneck stations in the network. This is the first rigorous analysis
that allows to favorably compare splitting against directly computing the overflow
probability of interest, which can be evaluated by solving a linear system of equations
with O(nd ) variables.

1

Introduction

The development of rare-event simulation algorithms for overflow probabilities in stable open
Jackson networks has been the subject of a substantial amount of papers in the literature
during the last decades (see Section 2 for the specification of an open Jackson network).
A couple of early references on the subject are [20] and [1]. Subsequent work which has
also been very influential in the development of efficient algorithms for overflows of Jackson
networks include [22], [12], [13], [16], [14], [9], [19], [11] and [8]. The survey papers of [15]
and [7] provide additional references on this topic.
The two most popular approaches that are applied to the construction of efficient rareevent simulation algorithms are importance sampling and splitting (see [3]). Importance
sampling involves simulating the system under consideration (in our case the Jackson network) according to a different set of probabilities in order to induce the occurrence of the
rare event. Then, one attaches a weight to each simulation corresponding to the likelihood
ratio of the observed outcome relative to the nominal / original distribution. In splitting, on
the other hand, there is no attempt to bias the behavior of the system. Instead, the idea is
to decompose the occurrence of the rare event of interest (in our case overflow in a Jackson
network) into a sequence of nested “milestone” events whose subsequent occurrence is not
rare. The rare event occurs when the last of the milestone events occurs. The idea is to
1

2
keep splitting the particles as they reach subsequent milestones. Of course, each particle is
attached a weight corresponding to the total number of times it has split so that the overall
estimation (which is the sum of the weights corresponding to the particles that make it to
the last milestone) provides an unbiased estimator of the probability of interest.
The most popular performance measure for efficiency analysis of rare-event simulation
algorithms for Jackson networks corresponds to that of “asymptotic optimality” or “weak
efficiency”. In order to both explain the computational complexity implied by this notion and
to put in perspective our contributions let us discuss the class of problems we are interested
in. Starting from any fixed state, we consider the problem of computing the probability that
the total number of customers in any fixed set of stations in the network reaches level n prior
to reaching the origin. That is, the probability that the sum of the queue lengths in any
given set of stations reaches level n within a busy period. The number of stations in the
whole network is assumed to be d and the number of bottleneck stations (i.e. stations with
the maximum traffic intensity in equilibrium) is β.
Weak efficiency guarantees that a subexponential number of replications (as a function of
the overflow level, say n) suffices to compute the underlying overflow probability of interest
within a given relative accuracy. In contrast, as we shall explain in Section 2, overflow
probabilities in the setting of Jackson networks can be computed by solving a linear system
of equations
with O(nd )1 unknowns. It is well known that Gaussian elimination then takes

O n3d operations to find an exact solution to such linear system. Moreover, since in our
case the associated linear system
has some sparsity properties the linear equations can be

3d−2
solved in as many as O n
operations (see the discussion in Section 2). Our analysis for
the solution of the associated linear system of equations is not intended to be exhaustive.
Our objective is simply to make the point that naive Monte Carlo (which indeed takes an
exponential number of replication in n to achieve a given relative accuracy) is not the natural
benchmark that one should be using in order to test the performance of an efficient simulation
estimator for overflows in Jackson networks. Rather, a more natural benchmark is the
application of a straightforward method for solving the associated system of linear equations.
It would be interesting to provide a detailed study of various methods for solving linear system
of equations (such as multigrid procedures) that are suitable for our environment and even
combine them with the ideas behind efficient simulation procedures. This, however, would
be the subject of an entire paper and therefore is left as a topic for future research.
Our goal here is to analyze a class of splitting algorithms similar to those introduced in
[22] for the evaluation of overflow probabilities at level n. Further analysis was given in [8],
where the authors provide necessary and sufficient conditions for the design of the “milestone
events” in order to achieve subexponential complexity in n. Our contribution is to show
that if the milestone
 events are properly placed as suggested by [8], the splitting algorithm
2β+1
requires O n
function evaluations to achieve a fixed relative error. Since clearly
the

number of bottleneck stations β is at most d, the complexity of splitting is O n2d+1 , which
is substantially smaller than that of the direct solution of the associated linear system. Our
analysis therefore provides theoretical justification for the superior performance observed
1

Given two non-negative functions f (·) and g (·), we say f (n) = O (g (n)) if there exists c, n0 ∈ (0, ∞)
such that f (n) ≤ cg (n) all n ≥ n0 .

3
when applying splitting algorithms compared to directly solving the associated linear system.
We believe that our results shed light into the type of performance that can be expected
when applying particle algorithms beyond the setting of Jackson networks. This feature
should be emphasized, specially given the fact that a linear time algorithm for computing
overflows in Jackson networks has been developed very recently (see [4]). Contrary to particle
methods, which are versatile and that can in principle be applied in great generality, the
algorithm in [4] takes advantage of certain properties of Jackson networks which are not
shared by all classes of systems.
In addition, our results also provide interesting connections to recent performance analysis
studied in the context of state-dependent importance sampling algorithms for a class of
Jackson networks. These connections might eventually help guide the users of rare event
simulation algorithms decide when to apply importance sampling or splitting. For instance,
consider the overflow at level n of the total population of a tandem network with d stations.
The work of [9] proposes an importance sampling estimator based on the subsolution of an
associated Isaacs equation. In particular, [9] shows that if exponential tiltings are applied
using the gradient of the associated subsolution as the tilting parameter (depending on the
current state), the corresponding algorithm is weakly efficient. Turns out that there are many
subsolutions that can be constructed varying certain so-called “mollification parameters”. A
recent analysis based on Lyapunov inequalities given
 in [6] shows that a natural selection
2(d−β)+1
of mollification parameters guarantees O n
function evaluations to achieve a given
relative error.Our analysis here therefore guarantees that one can achieve a running time of
order O nd+1 if one chooses importance sampling when there are more than d/2 bottleneck
stations in the network and splitting if there are less than d/2 bottleneck stations. Although
our analysis is still not sharp we believe that our results provide a significant step in order
to understand the connections between splitting and importance sampling.
The rest of the paper is organized as follows. A brief discussion on complexity and
efficiency considerations is given in Section 2. Then we discuss the necessary large deviations
asymptotics for Jackson networks required for our analysis in Section 3. The introduction
of the splitting algorithm as well as connections to the theory developed in [8] is given in
Section 4. Our complexity analysis is finally given in Section 5.

2

Complexity and Efficiency

We shall review concepts of efficiency and complexity in rare event simulation. We start
our discussion in the context of a generic class of rare event simulation problems. Consider
a sequence of events {En , n = 1, 2...} with pn , P (En ) → 0 as n → ∞ (Without loss
of generality, we might assume that pn → 0 exponentially fast as n ր 0.) The design of
an efficient rare-event simulation algorithm is typically associated with the construction of
an unbiased estimator, say p̂n , such that pn = E [p̂n ]. A number of m i.i.d. replications
(1)
(m)
{p̂n , ..., p̂n } is produced, the average of which forms an estimate of pn , namely
m

p̂n (m) =

1 X (j)
p̂ .
m j=1 n

4
By virtue of Chebyshev’s inequality we obtain the following property for the relative error,
|p̂n (m) − pn |/pn , of the estimate
Var (p̂n )
.
(1)
mp2n ǫ2
Hence, for a pre-determined upper bound ǫ of relative error, if we choose the number of
replications m such that
m ≥ ǫ−2 δ −1 (cvn )2 ,
(2)
P (|p̂n (m) − pn |/pn > ǫ) ≤

where cvn2 = Var (p̂n ) /p2n is the squared coefficient of variation of pˆn , we can guarantee that
the relative error is no larger than ǫ with probability at least 1 − δ.
Equation (2) stipulates that m needs to grow at least at the same rate as cvn2 does in
order to keep the relative error within a desirable threshold. If cvn2 grows at a subexponential rate (i.e. if log (cvn )2 = o (log pn ) , as n ր ∞) the estimator is said to be asymptotic
optimality, logarithmical efficient or weakly efficient. In this case, the number of replication
needs to increase subexponentially in n to achieve a prescribed level of relative accuracy.
The name “asymptotic optimality” is derived from the fact that weak efficiency implies that
the exponential rate of decay to zero of the Ep̂2n coincides with that of p2n and therefore is
maximal (by virtue of Jensen’s inequality).
Obviously, one has to keep in mind that weak efficiency measures the optimality of the
estimator for a given level of computational budget. Coming to the splitting algorithm, it
is apparent that the computational effort varies drastically with the degree of splitting performed; one must therefore take into account the cost involved in generating each replication
of p̂n . We measure such cost in terms of the number of elementary function evaluation which
we will take to be simple addition, multiplication, comparison, and the generation of a single
uniform random variable. When we incorporate the computational cost per replication of
the estimator, (2) says that the total number of function evaluations needed has to keep pace
with the work-normalized squared coefficient of variation, i.e., cvn2 · Nn , where Nn is the cost
per replication of p̂n . We will show in section 5 that Nn is closely related to the expected
total number of the survival particles in a single run of the Splitting algorithm.
Coming back to the setting of Jackson networks. It is important to recognize that overflow probabilities in such a setting can be obtained by solving a system of linear equations.
So, a reasonable benchmark for any simulation based algorithm to be regarded “efficient” is
indeed how fast one can solve such a system of linear equation by a direct procedure. Jackson
networks are basically multidimensional simple random walks with constrained behavior on
the boundaries. In particular, they are Markov chains living on a countable state-space. The
overflow probabilities can be conveniently expressed as first passage time probabilities, which
in turn can be characterized as the solution to certain linear system of equations thanks to its
countable state-space Markov chain structure. We shall quickly review how to obtain such
linear system for a generic Markov chain Q = {Qk : k ≥ 0} living on a countable state-space
S with transition matrix {K (x, y) : x, y ∈ S}. Let A, B be two disjoint subsets of S, define
TA , inf{k ≥ 0 : X ∈ A}, TB , inf{k ≥ 0 : X ∈ B} and put p (x) = Px (TA ≤ TB ). A
simple conditioning argument on the first transition leads to
X
p (x) =
K (x, y) p (y)
(3)
y∈S

5
subject to the boundary conditions
p (x) = 1 for x ∈ A,

p (x) = 0 for x ∈ B.

In fact, p (·) is the minimum non-negative solution to the above system (see [5]).
Now, if Q describes the state of the embedded discrete time Markov chain corresponding
to a Jackson network with d stations then S = Z+d . The transition dynamics of a Jackson
network are specified as follows (see [21] p. 92). Inter-arrival times and the service times
are all independent and exponentially distributed random variables. The arrival rates are
given by the vector λ = (λ1 , ..., λd )T and service rates are given by µ = (µ1 , ..., µd )T . (By
convention all of the vectors in this paper are taken to be column vectors and T denotes
transposition.) A job that leaves station i joins station j with probability Pi,j and it leaves
the system with probability
d
X
Pi,j .
Pi,0 , 1 −
j=1

The matrix P = {Pi,j : 1 ≤ i, j ≤ d} is called the routing matrix. We shall consider open
Jackson networks, which satisfy the following conditions:
i) ∀i, either λi > 0 or λj1 Pj1 j2 ...Pjk i > 0 for some j1 , ..., jk .
ii) ∀i, either Pi0 > 0 or Pij1 Pj1 j2 ...Pjk 0 > 0 for some j1 , ..., jk .
iii) The network is stable (i.e. a stationary distribution exists).
These conditions simply require that each station will receive jobs either directly from the
outside or routed from other stations, and each job will leave the system eventually. Our main
interest lies in the evaluation of pn (x) assuming that B = {0} and An = {y : v T y = n} where
v is a binary vector which encodes a particular subset of the network (i.e., the i-th position of
the vector v is 1 if station i falls in the subset of interest, and 0 otherwise). We shall denote
by V (x) = xT v the mapping recording the total population in the stations corresponding to
the vector v. The case in which v = 1 = (1, 1, ..., 1)T corresponds to the total population of
the system. So, pn (x), or more precisely pVn (x), corresponds to the overflow probability in
the subset encoded by v within a busy period and starting from x and that. In this setting,
it follows (as we shall review in the next section) that
pVn (x) −→ 0 exponentially fast in n as

d
n ր ∞ and the system of equations (3) has O n unknowns. Gaussian elimination requires
O n3d function evaluations to find the solution of such system. But since each state of the
Markov chain in this case has possible interactions with only a small fraction of the entire
state-space, it is therefore possible to permutate the states (say in lexicographic order) so
that the system is banded (i.e. the associated matrix is sparse in the sense that
 its non-zero
d−1
entries fall to a diagonal band.) One can show that
the
bandwidth
is
O
n
, and therefore




2
solving such a banded linear system requires O nd · nd−1
= O n3d−2 operations (see,
e.g., [2]).

6
Estimators that possess weak efficiency (in a work-normalized sense) are guaranteed to
run at subexponential complexity. When comparing to the above polynomial algorithms
of solving systems of linear equations, the efficiency analysis of such estimators appears
to be insufficient. We will show in later analysis that the multilevel Splitting algorithm
suggested by Dean and Dupuis [8], applied to estimate the overflow probabilities in Jackson
networks, requires fewer function evaluations than directly solving the associated system of
linear equations.

3

Jackson Networks: Notation and Properties

As we mentioned in the previous section, a Jackson network is encoded by two vectors of
arrival and service rates, λ = (λ1 , ..., λd )T and µ = (µ1 , ..., µd )T , togetherPwith a routing
matrix P = {Pi,j : i ≤ i, j ≤ d}. Without loss of generality, we assume that di=1 (λi + µi ) =
1. The network is assumed to be open and stable so conditions i), ii), and iii) described in
the previous section are in place.
Given the stability assumption, the system of equations given by
φ i = λi +

d
X

φj Pji ,

∀i = 1, 2, ..., d

(4)

j=1

admits a unique solution φ = λT (I − P )−1 (see [3]). The traffic intensity at station i in the
system in equilibrium is given by ρi which is defined by

 T
λ (I − P )−1 i
φi
=
(5)
ρi =
µi
µi
and satisfies ρi ∈ (0, 1) for all i = 1, 2, ..., d. Define ρ∗ = max1≤i≤d ρi and let β be the
cardinality of the set {i : ρi = ρ∗ }.
We shall study the queueing network by means of the embedded discrete time Markov
chain Q = {Q(k) : k ≥ 0}, where Q(k) = (Q1 (k), . . . , Qd (k)). For each k, Qi (k) represents
the number of customers in station i immediately after the kth transition epoch of the system.
As mentioned before, the process Q lives in the space S = Z+d .
Let V (y) = y T v be the total population in the stations corresponding to the binary vector
v. We are interested in the overflow probability in any given subset of the Jackson network.
More precisely, we wish to estimate
pVn = P { total population in stations encoded by v reaches
n before returning to 0, starting from 0}.
In turn, pVn can be expressed in terms of the following stopping times,
T{0} , inf{k ≥ 1 : Q (k) = 0},
TnV , inf{k ≥ 1 : V (Q (k)) ≥ n}.

7
Indeed, if we use the notation Px (·) , P(·|Q(0) = x) then we can rewrite pVn as
pVn = P0 (TnV ≤ T{0} ).

(6)

pVn (x) = Px (TnV ≤ T{0} ).

(7)

Similarly,
pVn

The asymptotic analysis of
(x) can be studied by means of large deviations theory.
We shall indicate how this theory can be applied to specify an efficient splitting algorithm
in the next section. In the mean time, let us provide a representation for the dynamics of
the queue length process that will be convenient in order to motivate the elements of the
efficient splitting algorithm that we shall analyze.
As we mentioned earlier, Jackson networks are basically constrained random walks. The
constraints arise because the number of customers in each station must be non-negative.
Thinking about Jackson networks as constrained random walks facilitates the introduction
and motivation of the necessary large deviations elements behind the description of the
splitting algorithm. In order to specify the dynamics of the embedded discrete time Markov
chain in terms of a random walk type representation we need to introduce notation which
will be useful to specify the transitions at the boundaries induced by the non-negativity
constraints.
The state-space Z+d can be partitioned in 2d different regions which are indexed by all
the subsets E ⊆ {1, . . . , d}. The region encoded by a given subset E is defined as
∂E = {z ∈ Zd+ : zi = 0, i ∈ E, zi > 0, i ∈
/ E}.
The interior of the domain is given by ∂∅ and the origin is represented by ∂{1,2,...,d}. Subsets
other than the empty set represent the “boundaries” of the state-space and correspond to
system configurations in which at least one station is empty. The collection of all possible
values that the increments of the process Q can take depends on the current region at which
Q is positioned. However, in any case, such collection is a subset of
V , {ei , −ei + ej , −ej : i, j = 1, 2, ..., d},
where ei is the vector whose i-th component is one and the rest are zero. An element of
the form ei represents an arrival at station i, an element of the form −ei + ej represents a
departure from station i that flows to station j and an element of the form −ej represents a
departure from station i out of the system. The set of all possible departures from station i
is a subset of
V−
i , {w : w = −ei or w = −ei + ej for some j = 1, ..., d}.
Because of the nonnegativity constraints on the boundaries of the system we have to be
careful when specifying the transition dynamics. First we define a sequence of i.i.d. random
variables {Y (k) : k ≥ 1} so that for each w ∈ V

if w = ei ,
 λi
P (Y (k) = w) =
µi Pij
if w = −ei + ej ,
.

µi Pi0
if w = −ei .

8
The dynamics of the queue-length process admit the random walk type representation given
by
Q(k + 1) = Q(k) + ζ (Q(k), Y (k + 1)) ,
(8)
where ζ (·) is the constrained mapping and it is defined for x ∈ ∂E via

0
if w ∈ ∪i∈E V−
i
ζ (x, w) ,
.
w
otherwise
The large deviations theory associated to Jackson networks is somewhat similar (at least
in form) to that of random walks. One has to recognize, of course, that the non-smoothness of
the constrained mapping as a function of the state of the system creates substantial technical
complications, but we will leave aside this issue in our discussion because our objective is
simply to describe the form of the necessary large deviations results for our purposes. An
extremely important role behind the development of large deviations theory for light-tailed
random walks is played by the logmoment generating function of the increment distribution.
So, given the similarities suggested by the dynamics of (8) and those of a simple random
walk it is not surprising that the logmoment generating function of the increments, namely,


ψ (x, θ) , log E exp θT ζ (x, Y (k))
(9)

also plays a crucial role in the large deviations behavior of pVn (x) as n ր ∞.
In order to understand the large deviations behavior of pVn it is useful to scale space
by 1/n, thereby introducing a scaled queue length process {Qn (k) : k ≥ 0} which evolves
according to
1
Qn (k + 1) = Qn (k) + ζ (Qn (k), Y (k + 1)) .
n
Suppose that Qn (0) = y = x/n and note that
T{0} , inf{k ≥ 1 : Qn (k) = 0}, TnV , inf{k ≥ 1 : V (Qn (k)) ≥ 1}.
Note that using the scaled queue length process one can write


1
V
V
pn (y) = E pn (y + ζ (y, Y (1))) .
n

(10)

Large deviations theory dictates that
pVn (y) = exp (−nWV (y) + o (n))

(11)

as n ր ∞ for some non-negative function WV (·). In order to characterize WV (·) we can
combine the previous expression together with (10) and a formal Taylor expansion to obtain


1
1
V
E pn (y + ζ (y, Y (1)))
1= V
pn (y)
n
1
≈ E exp{−nWV [y + ζ (y, Y (1))] + nWV (y)}
n
= E exp{−∂WV (y)ζ (y, Y (1)) + o (1)}
= exp (ψ (y, −∂WV (y)) + o (1)) .

9
Sending n ր ∞ we formally arrive at the equation
ψ (y, −∂WV (y)) = 0

(12)

together with the boundary condition WV (y) = 0 if V (y) ≥ 1. The previous equation
is the so-called Isaacs equation which characterizes the large deviations behavior of pVn (·)
and it was introduced together with a game theoretic interpretation by Dupuis and Wang
in [10]. The solution to (12) is understood in a weak sense because the function WV (·) is
typically not differentiable everywhere. Nevertheles, it coincides with a certain calculus of
variations representation which can be obtained out of the local large deviations rate function
for Jackson networks (see [17]).
An asymptotic lower bound for WV (y) can be obtained by finding an appropriate subsolution to the Isaacs equation, in which the equality signs in (12) are appropriately replaced by
inequalities thereby obtaining a so-called subsolution to the Isaacs equation. In particular,
W V (·) is said to be a subsolution to the Isaacs equation if
ψ(y, −∂W V (y)) ≤ 0

(13)

subject to W V (y) ≤ 0 if V (y) ≥ 1. The subsolution property guarantees W V (y) ≤ WV (y),
which translates to an asymptotic logarithmic upper bound pVn (y). The subsolution is said
to be maximal at zero if W V (0) = WV (0). Not surprisingly, subsolutions are easier to construct than solutions and, as we shall discuss in the next section, beyond their use in the
development of asymptotic upper bounds they can be applied to the design of efficient simulation procedures. The use of subsolutions to the Isaacs equation for the design of efficient
simulation algorithms was introduced in [10]. A derivation of the subsolution equation (13)
following the same spirit leading to (12) using Lyapunov inequalities is given in [6].
As we mentioned in Section 2, the efficiency analysis of a rare-event simulation estimator
depends on the growth rate of its coefficient of variation. We are interested in an asymptotic
analysis that goes beyond the error term exp(o (n)) given by the large deviations approximation (11). So, we must enhance the large deviations approximations in order to provide a
more precise estimate for pVn . Developing such an estimate is the aim of the following proposition which follows as a direct consequence of Proposition 2 and the analysis in Section 5 in
[4]; see also Section 4 in this paper for a sketch of the proof.
Proposition 1. There exists K > 0 (independent of x and n) such that
pVn (x) ≤ KP {V (Q (∞)) = n}/P {Q (∞) = x},
where Q∞ is the steady state queue length. Moreover, if kxk ≤ c for some c ∈ (0, ∞) then2
pVn (x) = Ω[P {V (Q (∞)) = n}/P {Q (∞) = x}]
as n ր ∞.
2

Given two non-negative functions f (·) and g (·), we say f (n) = Ω (g (n)) if there exists c, n0 ∈ (0, ∞)
such that f (n) ≥ cg (n) all n ≥ n0 .

10
Remark It is important to keep in mind that we shall mostly work with the process Q (·)
directly, as opposed to the scaled version Qn (·) which is used in the analysis of [8].
The previous proposition provides the necessary means to estimate pVn up to a constant;
we just need to recall that the distribution of Q (∞) is computable in closed form (see [21]
p. 95). In particular, we have that
π (m1 , ..., md ) =

d
Y

P (Qj (∞) = mj )

j=1

=

d
Y

m

(1 − ρj ) ρj j ,

j = 1, ..., d, and mj ≥ 0.

j=1

We shall use π (·) to denote the stationary measure of Q. In simple words, the previous
equation says that the steady state queue length process has independent components which
are geometrically distributed. In particular, P (Qj (∞) = m) = ρm
j (1 − ρj ) for m ≥ 0. The
next proposition follows directly from standard properties of the geometric distribution (see
Proposition 3 in [4]).

−nγV βV −1
Proposition 2. P [V (Q (∞))
=
n]
=
Θ
e
n
, where γV = − log ρV∗ , in which ρV∗ =
P
max{ρi : vi = 1}; and βV = i I{ ρi = ρV∗ , vi = 1} is the number of bottleneck stations in the
target subset corresponding to v.

4

The Splitting Algorithm

The previous section discussed some large deviations properties required to guide the construction of an efficient splitting scheme using the theory developed in the work of Dean and
Dupuis [8]. In order to explain the construction suggested by Dean and Dupuis let us first
discuss the general idea behind the splitting algorithm that we shall analyze; a variation of
which was first applied to Jackson networks by by Villen-Altamirano and Villen-Altamirano
[18].
The strategy is to divide the state-space into a collection of regions {Cjn : 0 ≤ j ≤ ln (x)}
which are nested and that help define “milestone” events that interpolate between the initial
position of the process and the target set, which corresponds to the region C0n . That is, in
our setting we put C0n , {y ∈ S : V (y) ≥ n} and the remaining Cjn ’s are placed so that
n
. How to construct the level sets Cjn in order to induce efficiency will be
C0n ⊆ C1n ⊆ ... ⊆ CM
n
discussed below. An observation that is intuitive at this point, however, is that one should
have Mn = Θ (n) so that the next milestone event becomes accessible given the current level.
For the moment, let us assume that the Cjn ’s have been placed. The splitting algorithm
proceeds as follows.
Algorithm SA
3

Given two positive functions f (·) and g (·), recall that f (n) = Θ (g (n)) if f (n) = O (g (n)) and f (n) =
Ω (g (n)).

11
1.- Initiate the simulation procedure with a single particle starting from position x ∈ Ckn
for a given k ≥ 1. Let w1 = 1 be the initial weight associated to such particle.
n
2.- Evolve the initial particle until either it hits {0} or it hits level Ck−1
. If the particle
n
hits {0}, then the particle is said to die. If the particle reaches level Ck−1 then it is replaced
by r identical particles (for a given integer r > 1). The replacing particles are called the
immediate descendants or children of the initial particle, which in turn is said to be their
parent. The children are positioned precisely at the place where the parent particle reached
n
level Ck−1
. The weight wj associated to the j-th children (enumerate the children arbitrarily)
has a value equal to the weight of the parent particle multipled by 1/r.
3.- The procedure starting from step 1 is replicated for each of the offspring particles in
place; carrying over the value of each of the weights at each level for the surviving particles
(the weights of the particles that die can be disregarded).
4.- Steps 1 to 3 are repeated until all the particles either die or reach level C0n .
Dean and Dupuis in [8] show how to apply large deviations theory to select the Cjn ’s in
order to obtain a weakly efficient splitting algorithm. One needs to balance the number of
the Cjn ’s so that it is not unlikely for a given particle to reach the next level while keeping
the total number of particles controlled. We now provide a formal motivation for the use of
large deviations for constructing the Cjn ’s in a balanced way.
It is convenient, as we did in our formal large deviations discussion in the previous section,
to consider the scaled process Qn (·). Let us assume that the splitting mechanism indicated
in Algorithm SA is in place and that our initial position is set at level Q (0) = x, so that
Qn (0) = y = x/n. The Cjn ’s are typically constructed in terms of the level sets of a so-called
importance function which we shall denote by U (·). In particular, put Dn , {y ∈ n−1 S :
V (y) < 1} and set Cjn = nLzn (j) , where
Lz , {y ∈ Dn : U (y) ≤ z},

(14)

and the zn (j)’s are appropriately chosen momentarily. Then, define
ln (y) = min{j ≥ 0 : ny ∈ Cjn } = min{j ≥ 0 : x ∈ Cjn }.

(15)

The total weight corresponding to a particle that reaches level C0n given that it started at
level ln (y) is r −ln (y) . In order to have at least a weakly efficient algorithm we wish to achieve
two constraints. The first one imposes the total weight of a particle reaching level C0n to be
pVn (x) exp (o (n)); this would guarantee that the second moment of the resulting estimator
achieves asymptotic optimality. The second constraint dictates that the expected number of
particles that make it to C0n , which is roughly r ln (y) pVn (x) exhibits subexponential growth (i.e.
exp (o (n))); this would guarantee a cost per replication that is subexponential. Note that
both constraints lead to the requirement of r ln (y) pVn (x) = exp (o (n)). So, given a subsolution
W V (·) to the corresponding Isaacs equation, which implies that

pVn (x) ≤ exp −nW̄V (y) + o (n) ,

it suffices to ensure that

ln (y) log (r) − nW V (y) = o (n) .

(16)

12
The behavior of ln (y) as a n ր ∞ only relates to the properties of the function U (·) and it
is really independent of the large deviations behavior of the system. In particular, picking
zn (j) = ∆j/n, ∆ ∈ (0, 1] yields ln (y) = ⌈nU (y) /∆⌉ and therefore, equation (16) suggests
that one should select U (y) = ∆W V (y) / log (r) with W V (0) = WV (0) in order to obtain
a weakly efficient estimator for pVn . This is precisely the conclusion obtained in the work
of [8] who present a rigorous analysis that justifies the previous heuristic discussion. Our
development in the next section will sharpen the efficiency properties of the sampler proposed
in [8] when applied to Jackson networks. So, we content ourselves with the previous heuristic
motivation for the splitting method that we will analyze in the next section and which in
turn is based on the viscosity subsolution given by
W̄V (y) = ̺T y − log ρV∗ ,

(17)

where ̺i = − log ρi for i = 1, ..., d, see e.g., [11] and [8].
We close this section with a precise definition of the estimator that we will analyze. First,
given a constant ∆ ∈ (0, 1] (the level size) define W̄V (·) as indicated in (17) for each y = x/n
with x ∈ S. Then, select an integer r > 1 and define U (y) = ∆W V (y) / log (r). Given the
initial position x put y = x/n and define the sets {Cjn : 1 ≤ j ≤ ln (y)} as indicated above
(see equation (15)). Run Algorithm SA and let Nn be the number of particles that survive
up to C0n ; their corresponding final weight is 1/r ln (x) . Our estimator for pVn (x) is simply
Rn (x) = Nn (x) /r ln (x)

(18)

Now, for the sake of analytical convenience, when analyzing the second moment of Rn (x)
we will adopt the so-called fully branching representation of the previous estimator (see [8]).
Such fully branching representation is obtained by splitting death particles at level zero. It
is useful to think about fully branching conceptually in the following terms. Start with a
particle at position x in level Clnn (x) and run step 2 of Algorithm SA, but instead of killing the
particle if it hits {0} before reaching Clnn (x)−1 , just allow it to also split into r particles and
update the weight of the children as indicated in Algorithm SA when the particle reaches
Clnn (x)−1 . Step 3 continues, now the particles that are sitting in level Clnn (x)−1 will evolve and
the death particles will once again split and remain in state 0 (so state 0 is being populated
with death particles). After ln (x) iterations we have r ln (x) total particles labeled 1, 2, ..., r ln(x) ,
each with weight 1/r ln(x) . We define Ij as the indicator function of the event that the j-th
P ln (x)
particle is in C0n so that Nn (x) = rj=1 Ij . The fully branching representation of Rn (x) is
simply
ln (x)
rX
−ln (x)
Ij .
(19)
Rn (x) = r
j=1

5

Analysis of Splitting Estimators

We are now in a good position to perform a refined efficiency analysis for the estimator Rn (x).
We shall break our analysis into two parts. The first part corresponds to the expected number

13
of particles generated per run and the second part deals with the second moment of Rn (x).
First, we will review a technique studied in [4] based on the corresponding time reversed
chain associated to the queue length process Q. For both quantities we are able to obtain
an upper bound. We are then able to reach the conclusion that this multilevel Splitting
algorithm substantially outperforms the direct polynomial time algorithm for solving the
associated system of linear equations.
Our analysis takes advantage of the time reversed process associated to the underlying
Jackson network which we shall now define. Given the transition matrix {K (x, y) : x, y ∈ S}
of the process Q, we define the reversed Markov chain Q̃ = {Q̃ (k) : k ≥ 0} via the transition
matrix K̃ (·):
K̃ (y, x) = K (x, y) π (x) /π (y) ,
for x, y ∈ S. It turns out that Q̃ also describes the queue length process of an open stable
Jackson network with stationary distribution equal to π (·), (see [21] p. 95). We will use
P̃x (·) to denote the probability measure in path space associated to Q̃ given that Q̃ (0) = x.
The following result is similar to that of Proposition 2 in [4]. However, our representation
in (20) is slightly more useful for our purposes.
Proposition 3.
pVn (x) =
=



P̃π Q̃ (0) ∈ C0n , T̃{x} ≤ T̃{0} , T̃{x} < T̃C0n

(20)



P̃π Q̃ (0) ∈ C0n , T̃{x} ≤ T̃{0} < T̃C0n

(21)

π(x)Px (T{x} ≥ TC0n ∧ T{0} )

π(x)P0 (T{x} ≤ TC0n ∧ T{0} )

where T̃C0n , inf{k ≥ 1 : Q̃ (k) ∈ C0n }, T̃{x} , inf{k ≥ 0 : Q̃ (k) = x}, TC0n , inf{k ≥ 1 :
Q (k) ∈ C0n } and T{x} = inf{k ≥ 0 : Q (k) = x}. Moreover, there exists δ > 0 (independent
of x and n) such that
Px (T{x} ≥ TC0n ∧ T{0} ) ≥ Px (T||x|| ≥ T{0} ) ≥ δ,

(22)

where ||x|| is the L1 norm of x.
Proof. We assume that x 6= 0. The case x = 0 is included in the analysis of (21). First, we
observe that


pVn (x) = Px TC0n < T{0} , T{x} < TC0n ∧ T{0} + Px TC0n < T{0} , T{x} ≥ TC0n ∧ T{0}


= pVn (x) Px T{x} < TC0n ∧ T{0} + Px TC0n < T{0} , T{x} ≥ TC0n ∧ T{0} .

Therefore,


n < T{0} , T{x} ≥ TC n ∧ T{0}
P
T
x
C
0
0

.
pVn (x) =
Px T{x} ≥ TC0n ∧ T{0}

14
Following the same technique as in Proposition 2 in [4] we have that

π (x) Px TC0n < T{0} , T{x} ≥ TC0n ∧ T{0}
∞
X

π (x) Px TC0n < T{0} , T{x} ≥ TC0n ∧ T{0} , TC0n = k
=
k=0

=

∞
X
k=1

=

∞
X

π (x)

X

y0 =x,y1 ,..,yk

X

k=1 y0 =x,y1 ,..,yk

K (y0 , y1 ) × ... × K (yk−1, yk ) I k < T{0} , T{x} ≥ k ∧ T{0} , TC0n = k




K ′ (y1 , y0 ) × ... × K ′ (yk , yk−1) π (yk ) I k < T{0} , T{x} ≥ k ∧ T{0} , TC0n = k .

Letting yi′ = yk−i for i = 1, ...k we see that the summation in each of the terms above ranges
over paths y0′ , ..., yk′ satisfying that y0′ ∈ C0n , T̃{x} = k (so in particular yk′ = x) and also that
T̃{0} ≥ k, T̃C0n > k. So, we can interpret the previous sum as


P̃π Q̃ (0) ∈ C0n , T̃{x} ≤ T̃{0} , T̃{x} < T̃C0n .
This yields part (20). Part (21) corresponds to Proposition 2 of [4]. Finally, (22) also follows
as in Proposition 7 of [4].

Proposition 1 and 2 from Section 2 follow as a consequence of this result, the rest of the
details are given in Section 5 of [4].
Given the subsolution we proposed in Section 4, the importance function can be written
as


1 T
∆
∆
V
=
̺ x − log ρ∗
U (x/n) = W̄V (x/n)
log r
n
log r


1
= C ∆ − αT x∆ ,
(23)
n
where C = − log ρV∗ / log r, and α = ̺ / log ρV∗ . The level index function also simplifies to



1 T
nU (x/n)
n
⌉ = ⌈nC 1 − α x ⌉ = ⌈C n − αT x ⌉.
(24)
l (x) = ⌈
∆
n
We shall first look at the expected number of surviving particles of the splitting algorithm
which characterizes the stability of the algorithm. One shall keep in mind that when the
complexity of the splitting algorithm is concerned, what actually matters is the total function
evaluation involved in each run. An upper bound is obtained for this quantity, as measured
by the sum of all particles generated at interim levels weighted by the maximum remaining
function evaluations associated with each of them. We first have the following result.

Proposition 4. The expected terminal number of particles for the splitting algorithm specified
by (∆, U) above satisfies

E [Nn (x)] = Θ nβV −1
(25)

where βV , introduced in Proposition 2, denotes the number of bottleneck stations corresponding to the vector v.

15
Proof. Note that
E [Nn (x)] = r l

n (x)

pVn (x) .

V
−γV n βV −1
In Proposition 2 we know that pVn (x) = Θ
e
n
. Since e−γV = elog ρ∗ = e−C log r =

r −C , we can write pVn (x) = Θ r −nC nβV −1 . Hence, plug in ln (x) = ⌈C n − xT v ⌉, we have


E [Nn (x)] = Θ r −nC nβV −1 r ⌈nC⌉ = Θ nβV −1 .
As pointed out earlier, the number of terminal surviving particles, although serves as a
reasonable proxy to measure the stability of the algorithm, is not suitable for quantifying the
complexity. We also need to take into account the number of function evaluations required
to generate Rn (x). The next result addresses precisely this issue.
Proposition 5. The expected computational
effort per run required to generate a single

βV +1
replication of Rn (x) is O n
.

n
Proof. To see this, let Nm
, m = 0, ...., ln (x), be the number of particles that survive to level
Clnn (x)−m . Also let ηm,j be the remaining computational effort of the j-th particle at the start
of the m-th level until it either reaches the next level or it dies out. Put η̄m,j (xj ) to be
the expectation of ηm,j given that the position of the j-th particle at the start of level m is
xj . Note that the norm of the position of xj is less than c · m for a given constant c that
depends on the traffic intensities of the system but not on the position of the particle per-se.
Therefore, it is easy to see that

sup η̄m,j (xj ) ≤ c · m,

n
1≤j≤Nm

for some c ∈ (0, ∞). Intuitively, each particle at level m either advances to the next level, or
it dies out by hitting the zero level before moving to the next one, since it takes Θ (1) work
to cross one single layer, ηm,j is dominated by the work required to die out, and hence its
mean is bounded from above by c × m for some constant c. The expected total work per run
is then given by


" Nn
#
n
ln (x)−1 Nm
ln (x)−1
m
X X
X
X
E
ηm,j  =
E
η̄m,j (xj )
m=0

m=0

j=1

j=1

ln (x)−1

≤

X

n
E [Nm
]·c·m

m=0

≤c·

ln (x)−1 

X

m=0

≤c·

ln (x)−1 

X

m=0

m βV −1 V  m
ρ∗ C r m m
C


m βV
= O nβV +1 ,
C

for some positive constant c where the first inequality is due to independence.

16
To facilitate the analysis of the second moment of Rn (x) we add the following notations.
We follow the analysis in [8] to make our exposition here self-contained. For a given generation
m, denote by Qm,j the position of the jth particle; recall that the accumulated weight up
to the mth stage of such a particle is r m . Let χm,j be the disjoint grouping of particles in
the next generation (i.e., m + 1) according to their “parents” in generation m. For k ∈ χm,j ,
denote by dk the offsprings of this particle at the final stage ln (x). We then have the following
expansion of the second moment of Rn (x):


ln (x)
rX

Ex 

j=1

ln (x)−1

=

X

m=0

2 

Ij r −ln (x)  


rm
X
Ex 

(26)

X

X

j=1 k,l∈χm,j ,k6=l

Imk r −ln (x)

mk ∈dk

!

X

ml ∈dl

!



ln (x)
rX

Iml r −ln (x)  + Ex 

j=1



Ij r −2ln (x)  ,

where we define Imk to be the indicator function of the event that particle mk is in the set
C0n . The second term above is essentially the diagonal terms of the second moment (26),
and for the off-diagonal terms, for each generation, we categorize particles according to their
common ancestors, a technique used by [8]. For the first term, we have

!
!
ln (x)−1
rm
X
X
X
X
X
Iml r −ln (x) 
Imk r −ln (x)
Ex 
m=0

j=1 k,l∈χm,j ,k6=l

ln (x)

=

X

Ex

m=0

·

X

" rm
X

k,l∈χm,j ,k6=l

I (V (Qm,j ) > 0) r −m

j=1

ml ∈dl

mk ∈dk

2

1 X
Imk r −(ln (x)−m−1)
r m ∈d
k

k

!

!
X
1
Iml r −(ln (x)−m−1)  .
r m ∈d
l

l

17
Conditioning on the whole genealogy up to step m, we obtain
" rm
X
2
Ex
I (V (Qm,j ) > 0) r −m
j=1

·

1 X
Imk r −(ln (x)−m−1)
r m ∈d

X

k,l∈χm,j ,k6=l

= Ex

" rm
X

k

k

I (V (Qm,j ) > 0) r −m

j=1



· E
= Ex

X

k,l∈χm,j ,k6=l

" rm
X

!

2

1 X
Imk r −(ln (x)−m−1)
r m ∈d
k

k

!
X
1
Iml r −(ln (x)−m−1) 
r m ∈d
l

!

I (V (Qm,j ) > 0) r −2m

l

!
X
1
Im r −(ln (x)−m−1) 
r m ∈d l
l

l

j=1

·

1
E
r X̄m,j

X

k,l∈χm,j ,k6=l

X

Imk r −(ln (x)−m−1)

mk ∈dk

!

1
E
r X̄m,j

X

ml ∈dl

!!
.
Iml r −(ln (x)−m−1)

−(ln
Note that EQm,j
= pVn (Qm,j ), and W =
mk ∈dk Imk r
Adding over m we obtain
2 

ln (x)
rX
Ij r −ln (x)  
Ex 

P


(x)−m−1)

P

k6=l

r −2 = r/ (1 − r).

j=1

ln (x)−1

=W

Ex

X



r −m Ex I (V (Qm,1 ) > 0) pVn (Qm,j )2 .

m=0

m=0

I (V (Qm,1 ) > 0) r −2m pVn (Qm,j )2

j=1

ln (x)−1

=W

" rm
X

X

#

Combining this with the diagonal term in (26), which can be readily expressed as r −ln (x) pVn (x),
we arrive at the following expansion for the second moment of Rn (x):
ln (x)−1



Ex Rn (x)

2

=W

X

m=0



r −m Ex I (V (Qm,1 ) > 0) pVn (Qm,1 )2 + r −ln (x) pVn (x) .

(27)



The next result takes advantage of expression (27) to obtain an upper bound for Ex Rn (x)2 .
Proposition 6. The second moment of Rn (x) satisfies

where β =

Pd

i=1


E [Rn (x)]2 = pVn (x)2 O nβ .

I (ρi = ρ∗ ) is the number of bottleneck stations in the whole network.

(28)

18
In order to prove the previous result, we will show that the second moment of Rn (x)
is dominated by the first item on the righthand side of the equality in (27). In turn, the
asymptotic behaviour of such term hinges on the conditional distribution of the exact position
of the particle in generation m, Qm,1 in Clnn (x)−m .
Proof. We begin the proof with an important property implied by the splitting algorithm:
V (Qm,1 ) > 0 ⇔ Qm,1 ∈ Clnn (x)−m = nL(ln (x)−m)∆/n
⇔ Qm,1 ∈ {z ∈ nDn : U (z/n) ≤ (ln (x) − m) ∆/n}




∆
1 T
⌈C n − αT x ⌉ − m }
⇔ Qm,1 ∈ {z ∈ nDn : C ∆ − α z∆ ≤
n
n




1
1
⇒ Qm,1 ∈ {z ∈ nDn : C 1 − αT z ≤
C n − αT x − m + 1 }
n
n
m−1
⇔ Qm,1 ∈ {z ∈ nDn : αT z ≥ αT x +
}
C
⇔ Qm,1 ∈ {z ∈ nDn : ̺T z ≤ ̺T x − (m − 1) log r}
(29)
where we used the representations of U (·) and ln (x) in (23) and (24) and the definition of
Lz in (14). In other words, if a particle survives m generations then its current position
is beyond the mth level, which implies that the weighted sum of system population, with
weight given by the vector ̺, is bounded from above by that of the initial position adjusted
, inf{k ≥ 1 : αT Q (k) ≥
by a linear function in m. If we define the stopping time T̂ m
C
αT x + m−1
} = inf{k ≥ 1 : ̺T Q (k) ≤ ̺T x − (m − 1) log r}, the above property also implies
C
that V (Qm,1 ) > 0 ⇒ T̂ m
< T0 . Now, the expectation term in the sum of (27) can be
C
expressed as


Ex I (V (Qm,1 ) > 0) pVn (Qm,1 )2



≤Ex I ̺T Qm,1 ≤ ̺T x − (m − 1) log r pVn (Qm,1 )2


 
=Ex pVn (Qm,1 )2 |̺T Qm,1 ≤ ̺T x − (m − 1) log r Px T̂ m
<
T
0
C
i 


h 
T
T
m < T
<
T
,
̺
Q
≤
̺
x
−
(m
−
1)
log
r
P
(30)
T̂
=Ex P2x TnV < T0 |QT̂ m |T̂ m
0
m,1
x
0
C
C
C

where we used the property derived in (29). For the first item in (30), we have
h 
i

T
T
Ex P2x TnV < T0 |QT̂ m |T̂ m
<
T
,
̺
Q
≤
̺
x
−
(m
−
1)
log
r
0
m,1
C
C

 2 n
π (C0 ) T
T
|̺ Qm,1 ≤ ̺ x − (m − 1) log r
≤KE 2
π (Qm,1 )
h
i

n 2
T
≤c1 nβV −1 ρV∗
Eπ e−2̺ Qm,1 |̺T Qm,1 ≤ ̺T x − (m − 1) log r

(31)

where c1 , K are some constants independent of n. Here we used Propositions 1 and 2 for the
the last two inequalities respectively. As for the expectation term in (31), since the process

19
Q (·) has for each dimension an increment at most of unit size, we can write
h
i
−2̺T Qm,1 T
T
Eπ e
|̺ Qm,1 ≤ ̺ x − (m − 1) log r
h
i
T
= Eπ e−2̺ Qm,1 |̺T x − (m − 1) log r − δ ≤ ̺T Qm,1 ≤ ̺T x − (m − 1) log r

≤ c2 exp −2̺T x + 2 (m − 1) log r


−2 m−1
m−1
V
C
= c3 exp −2
,
log ρ∗ = c3 ρV∗
C

(32)

where
 c2 , c3 and
 δ are some positive constants. It remains to give a bound on the term
Px T̂ m
< T0 . As a result of Proposition 2,
C


Px T̂ m
<
T
0
C


1
P ̺T Q (∞) ≤ ̺T x − (m − 1) log r
≤
π (x)


1
(m − 1)
T
T
=
P α̂ Q (∞) ≥ α̂ x +
π (x)
Ĉ

where Ĉ = − log ρ∗ / log r, ρ∗ = maxi ρi ∈ (0, 1) and α̂ = ̺ /ρ∗ = (log ρ1 / log ρ∗ , ..., log ρd / log ρ∗ )T .
Note that α̂i ∈ (0, 1). To finish the proof we need the following Lemma.
Lemma 1.

 


m−1
(m − 1)
T
T
= Θ P Z (β, 1 − ρ∗ ) ≥
P α̂ Q (∞) ≥ α̂ x +
Ĉ
"
#Ĉ
β−1
m−1
m−1
=Θ
(ρ∗ ) Ĉ
Ĉ
where Z (n, p) denotes a NBin (n, p) (negative binomial) random variable.
Proof. Note that
α̂T Q (∞) = Q (∞)T
=

d
X

̺
log ρ∗

Qi (∞) I (ρi = ρ∗ ) +

i=1

= Z (β, 1 − ρ∗ ) + W.

d
X
i=1

Qi (∞) I (ρi 6= ρ∗ )

log ρi
log ρ∗

One direction is elementary, since α̂T Q (∞) ≥ Z (β, 1 − ρ∗ ), we clearly have




(m − 1)
(m − 1)
T
T
T
P α̂ Q (∞) ≥ α̂ x +
≥ P Z (β, 1 − ρ∗ ) ≥ α̂ x +
.
Ĉ
Ĉ

(33)

20
For the other direction, note that there exists a constant ρ̃ < ρ∗ such that
W =

d
X
i=1

d

Qi (∞) I (ρi 6= ρ∗ )

X
log ρi
≤
Qi (∞) I (ρi 6= ρ∗ ) = Z (d − β, 1 − ρ̃) .
log ρ∗
i=1

As a result,
α̂T Q (∞) ≤st Z (β, 1 − ρ∗ ) + Z (d − β, 1 − ρ̃) ,
where “ ≤st ” denotes that the left hand side is stochastically dominated by the right hand
side. But since 1 − ρ∗ < 1 − ρ̃, a similar argument as given by Proposition 3 in [4] allows us
to obtain




(m − 1)
(m − 1)
T
T
T
P α̂ Q (∞) ≥ α̂ x +
≤ c0 P Z (β, 1 − ρ∗ ) ≥ α̂ x +
,
(34)
Ĉ
Ĉ
for some finite constant c0 that is independent of m. Combining (33) and (34), we have


 

(m − 1)
(m − 1)
T
T
T
P α̂ Q (∞) ≥ α̂ x +
= Θ P Z (β, 1 − ρ∗ ) ≥ α̂ x +
.
(35)
Ĉ
Ĉ
Using again Proposition 3 of [4], we reach the conclusion that
#
"


β−1
m−1
(m
−
1)
m
−
1
=Θ
P α̂T Q (∞) ≥ α̂T x +
(ρ∗ ) Ĉ
Ĉ
Ĉ

Going back to (30), the above Lemma allows us to write

β−1


m−1
m−1
< T0 ≤ c4
Px T̂ m
(ρ∗ ) Ĉ .
C
Ĉ

(36)

If we combine (31), (32) and (36), we obtain the following upper bound for the expectation
term in the sum of expression (27):

−2 m−1

C
Ex I (V (Qm,1 ) > 0) pVn (Qm,1 )2 ≤ c pVn (x)2 ρV∗
= c pVn (x)2 ρV∗

− m−1
C

= c pVn (x)2 r m−1







m−1

Ĉ
β−1
m−1

m−1
Ĉ

β−1

Ĉ
β−1

(ρ∗ )

m−1
Ĉ

(37)

where for the first equality we use the fact that ρV∗ ≤ ρ∗ and Ĉ ≤ C, and for the second
equality we use ρV∗ = r −C . Putting the bound in (37) back to the sum in the first item of

21
(27), we have
ln (x)−1

X

m=0

≤ cr



r −m Ex I (V (m,1 ) > 0) pVn (m,1 )2
ln (x)−1

−1

X

m=0
2

pVn

(x)


= pVn (x) O nβ

2



m−1
Ĉ

β−1

(38)

Finally, note that the second item of (27) is dominated by (38), and it follows immediately
that

E [Rn (x)]2 = pVn (x)2 O nβ .
Equipped with these results, we are ready to summarize our discussions in the statement
of the following Theorem, which is the main result of this paper.
Theorem 1. To estimate the overflow probability pVn (x) using Rn (x),
 the number of function
βV +β+1
evaluations needed for a given level of relative error is O n
.

Proof. Recall from section 2 that the number of function evaluations sufficient to achieve a
pre-determined level of relative accuracy for the Splitting estimator is proportional to the
work-normalized squared coefficient of variation. This is therefore immediate by combining
the upper bound analysis of the computational effort per run in Proposition 5 along with
the upper bound of the second moment of Rn (x) available in Proposition 6.

A direct comparison to the O n3d−2 complexity of solving a system of linear equations
(see Section 2 ) yields the immediate conclusion that the Splitting algorithm is “efficient” in
the sense that it improves over the “benchmark” polynomial algorithm. Even in the worst
case when we look at the total population of the network and the network is totally symmetric, i.e., all stations are bottlenecks (βV = β = d > 3), the number of flops needed is off by a
substantial factor nd−3 . In the case where βV = β = 1, the algorithm only requires a number
of function evaluations that at most grows cubically in the level of overflow n. Furthermore,
if the number of bottlenecks is less than half of the total number of stations,
i.e. β < d/2, the

d
Splitting algorithm enjoys a running time of order smaller than O n , which is not worst
than storing the vector that encodes the solution to the associated linear system. If, on the
other hand, more than half of the stations are bottlenecks, faster importance sampling based
algorithms do exist at least
 for the case of tandem networks; see the analysis in [6], which
2(d−β)+1
implies that O n
function evaluations suffice to obtain an estimator with a given
relative precision. Overall, the analysis thus provides some sort of guidance on the choice
of simulation algorithms. It is meaningful to point out that the previous comparison is not
based on the sharpest analysis. In fact we only resort to a rather crude upper bound in the
analysis of the second moment of Rn (x) in (31). A sharper result is possible by bounding
the expectation term in (20) with more care. But as pointed out in the Introduction, even

22
though there is still room for a more refined analysis, we believe our work provides substantial insights leading to a better understanding of the relations between these two classes of
algorithms.

References
[1] V. Anantharam, P. Heidelberger, and P. Tsoucas. Analysis of rare events in continuous
time marked chains via time reversal and fluid approximation. IBM Research Report,
REC 16280, 1990.
[2] P. Arbenz and W. Gander. A survey of direct parallel algorithms for banded linear
systems. Technical Report 221, Department Informatik,ETH Zurich, 1994.
[3] S. Asmussen and P. Glynn. Stochastic Simulation: Algorithms and Analysis. SpringerVerlag, New York, NY, USA, 2008.
[4] J. Blanchet. Optimal sampling of overflow paths in jackson networks. forthcoming, 2009.
[5] J. Blanchet and P. Glynn. Efficient rare-event simulation for the maximum of a heavytailed random walk. Ann. of Appl. Probab., 18:1351–1378, 2008.
[6] J. Blanchet, K. Leder, and P. Glynn. Lyapunov functions and subsolutions for rare
event simulation. Preprint, 2009.
[7] J. Blanchet and M. Mandjes. Rare event simulation for queues. In G. Rubino and
B. Tuffin, editors, Rare Event Simulation Using Monte Carlo Methods, pages 87–124.
Wiley, West Sussex, United Kingdom, 2009. Chapter 5.
[8] T. Dean and P. Dupuis. Splitting for rare event simulation: A large deviation approach
to design and analysis. Stochastic Processes and Its Applications, (119):562–587.
[9] P. Dupuis, A. Sezer, and H. Wang. Dynamic importance sampling for queueing networks.
Ann. Appl. Probab., 17:1306–1346, 2007.
[10] P. Dupuis and H. Wang. Importance sampling, large deviations, and differential games.
Stoch. and Stoch. Reports, 76:481–508, 2004.
[11] P. Dupuis and H. Wang. Importance sampling for jackson networks. Preprint, 2008.
[12] P. Glasserman, P. Heidelberger, P. Shahabuddin, and T. Zajic. Multilevel splitting for
estimating rare event probabilities, 1999.
[13] P. Glasserman and S. Kou. Analysis of an importance sampling estimator for tandem
queues. ACM TOMACS, 5:22–42, 1995.

23
[14] S. Juneja and V. Nicola. Efficient simulation of buffer overflow probabilities in jackson
networks with feedback. ACM Trans. Model. Comput. Simul., 15(4):281–315, 2005.
[15] S. Juneja and P. Shahabuddin. Rare event simulation techniques: An introduction and
recent advances. In S. G. Henderson and B. L. Nelson, editors, Simulation, Handbooks in
Operations Research and Management Science. Elsevier, Amsterdam, The Netherlands,
2006.
[16] D. Kroese and V. Nicola. Efficient simulation of a tandem jackson network. ACM Trans.
Model. Comput. Simul., 12:119–141, 2002.
[17] K. Majewski and K. Ramanan. How large queues build up in a jackson network. To
Appear in Math. of O.R., 2008.
[18] M.Villen-Altamirano and J. Villen-Altamirano. Restart: A method for accelerating
rare even simulations. In J.W. Colhen and C.D. Pack, editors, Proceedings of the 13th
International Teletraffic Congress. In Queueing, performance and control in ATM, pages
71–76. Elsevier Science Publishers, 1993.
[19] V. Nicola and T. Zaburnenko. Efficient importance sampling heuristics for the simulation
of population overflow in jackson networks. ACM Trans. Model. Comput. Simul., 17(2),
2007.
[20] S. Parekh and J. Walrand. Quick simulation of rare events in networks. IEEE Trans.
Automat. Contr., 34:54–66, 1989.
[21] P. Robert. Stochastic Networks and Queues. Springer-Verlag, Berlin, 2003.
[22] M. Villén-Altamirano and J. Villén-Altamirano. Restart: a straightforward method for
fast simulation of rare events. In Winter Simulation Conference, pages 282–289, 1994.

