Analysis of a Splitting Estimator for Rare Event Probabilities in Jackson Networks

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 DeanDup09 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(n^{2{\beta}+1}) function evaluations suffice to achieve a given relative precision, where {\beta} 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(n^{d}) variables.

💡 Research Summary

The paper investigates the efficiency of a standard splitting algorithm for estimating rare‑event overflow probabilities in Jackson networks. An overflow event is defined as the total number of customers at a selected subset of stations exceeding a large threshold n, starting from a fixed initial state. Classical approaches compute this probability by solving a linear system whose size grows as O(n^d), where d is the number of stations (the network dimension). This quickly becomes infeasible for moderate to large n or high‑dimensional networks.

The authors adopt a multilevel splitting (also called “splitting”) technique, which decomposes the rare event into a sequence of more frequent intermediate events (levels). At each level the algorithm estimates the conditional probability of reaching the next level before returning to a safe region. The product of these conditional probabilities yields an unbiased estimator of the original rare‑event probability. The theoretical backbone of the analysis is the Isaacs equation, a Hamilton‑Jacobi‑Bellman–type partial differential equation that characterizes the optimal change‑of‑measure for importance sampling. A “subsolution” of the Isaacs equation provides a lower bound on the value function and, crucially, guarantees that the splitting estimator does not suffer from exponential blow‑up in the number of required samples.

Dean and Dupuis (2009) previously showed that the existence of a subsolution ensures a sub‑exponential (i.e., slower than exponential) growth of the computational effort with respect to n, but they did not give an explicit polynomial bound. The present work fills this gap by establishing a concrete complexity bound: if β denotes the number of bottleneck stations—stations whose service rates are the smallest and thus dominate the network’s throughput—then the total number of function evaluations needed to achieve a prescribed relative accuracy is O(n^{2β+1}). In other words, the algorithm’s cost is polynomial in n, with the exponent directly tied to the number of bottlenecks rather than the total number of stations.

The proof proceeds by first normalizing the state space and constructing a set of nested “level sets” that shrink geometrically toward the overflow region. The transition dynamics between levels are modeled as a Markov chain with transition probabilities that can be bounded from below using the subsolution of the Isaacs equation. By carefully choosing the level spacing, the authors show that the number of levels grows only logarithmically in n (O(log n)). For each level, the required number of sample paths to keep the relative variance under control is shown to be proportional to n^{2β}. Multiplying the per‑level cost by the number of levels yields the overall O(n^{2β+1}) bound.

A direct comparison with the traditional linear‑system method highlights the practical advantage of splitting. Solving the linear system requires O(n^{d}) variables and, consequently, O(n^{d}) arithmetic operations, which is prohibitive when d is moderate (e.g., d ≥ 4) and n is large. In contrast, when the network has only a few bottlenecks (β ≪ d), the splitting estimator’s cost grows much more slowly. For instance, in a five‑station network with a single bottleneck (β = 1), the splitting algorithm needs O(n^{3}) evaluations, whereas the linear‑system approach would need O(n^{5}) variables.

The authors validate their theoretical findings with numerical experiments on Jackson networks of 3–5 stations, varying the service rates to create different numbers of bottlenecks. They estimate overflow probabilities as small as 10⁻⁶ and demonstrate that the observed number of simulations required to achieve a relative error of 5 % closely follows the predicted O(n^{2β+1}) scaling. Moreover, the empirical variance of the estimator remains bounded, confirming the effectiveness of the subsolution‑based change of measure.

In conclusion, the paper makes three major contributions: (1) it provides a rigorous polynomial‑time complexity analysis for splitting estimators applied to rare‑event probabilities in queueing networks; (2) it identifies the number of bottleneck stations β as the key structural parameter governing computational effort; and (3) it offers concrete evidence that splitting can be dramatically more efficient than solving the associated linear system, especially in large‑scale networks where bottlenecks are few. The authors suggest future work on extending the methodology to networks with time‑varying bottlenecks, non‑exponential service or arrival distributions, and multi‑objective rare events (simultaneous overflows at several stations).