Splitting Technique for Rare Event Simulation

• Rare event simulation by splitting technique is a Monte Carlo method that decomposes small probability events into a sequence of manageable intermediate events.
• It leverages large deviation theory and variational analysis to design optimal splitting levels, thereby controlling estimator variance and ensuring computational efficiency.
• Applications in queueing networks, reliability engineering, and both Markovian and non-Markovian systems demonstrate its broad adaptability and practical impact.

Rare event simulation by splitting technique refers to a large family of Monte Carlo algorithms designed to efficiently estimate extremely small probabilities of unlikely events—typically of the form that a Markov process reaches a rare set before entering another, more probable set. The essential principle is to decompose the rare event estimation task into a sequence of more tractable intermediate events by replicating (“splitting”) trajectories at predefined thresholds or levels, thereby amplifying sample representation in the important regions of the state space. Modern splitting techniques are grounded in large deviation theory, variational calculus, and precise control of estimator variance, providing both rigorous efficiency guarantees and concrete algorithmic design strategies.

1. Fundamental Principles of Splitting Techniques

Splitting (or multilevel splitting) transforms a rare event probability estimation problem—where straightforward Monte Carlo requires an infeasible number of samples—into a series of less rare and hence more computable events. The standard setting considers a Markov process $\{X_t\}$, and events of the form:

$$p(x) = P \Bigl\{ X_j \in B,\;\; X_i \notin A,\;\; 0 \leq i \leq j<\infty \;\bigg| X_0 = x \Bigr\}$$

where $A$ is a “bad” set and $B$ the rare “target.”

To estimate $p(x)$, the state space is partitioned into nested level sets

$C_0 = B \subset C_1 \subset \cdots \subset C_K$

typically defined via an importance function or level set function $V(\cdot)$: $C_i = \{ y: V(y) \leq z_i \}$ for thresholds $z_i$. When a simulated trajectory reaches a new level, it is split into several offspring, each inheriting a discounted statistical weight. The final rare event probability is estimated as a weighted sum over all particles that reach $B$.

This approach is fundamentally different from importance sampling (which changes probability measures); splitting modifies the trajectory ensemble to maintain unbiasedness while dramatically reducing variance.

2. Large Deviation Theory and Variational Formulation

The underlying assumption enabling rigorous analysis is that the rare event probability admits a large deviation scaling:

$$\lim_{n \to \infty} -\frac{1}{n} \log p^n(x_n) = W(x)$$

where $W(x)$ arises as a solution to a calculus of variations or optimal control problem, with

$$W(x) = \inf_{\phi:\phi(0)=x,\;\phi(t)\in B} \int_0^t L(\phi(s), \dot{\phi}(s)) ds$$

or equivalently as the value function solving a Hamilton–Jacobi–Bellman (HJB) PDE.

This large deviation rate function provides not just the exponential scaling, but also a methodology for designing both optimal importance sampling and splitting algorithms; the optimal set of splitting levels and the importance function structure emanate from $W(x)$.

3. The Subsolution Concept and Stability Criteria

A central analytic tool is the subsolution notion. Given the variational (or HJB) problem for $W(x)$, a function $\bar{W}: D \to \mathbb{R}$ is a subsolution if it satisfies:

• Boundary: $\bar{W}(x) \leq 0$ for all $x \in B$
• Difference inequality: $\bar{W}(x) - \bar{W}(y) \leq \mathcal{J}(x, y)$ for $x, y \notin A \cup B$, where $\mathcal{J}(x, y)$ is the transition cost.

A (maximal) subsolution obeying equality with $W(x)$ yields an optimally efficient splitting algorithm. The stability of the algorithm—meaning that the expected number of generated particles grows subexponentially (as opposed to exponentially) in the system size parameter—is equivalent to the existence of such a subsolution for the weighted importance function:

$\bar{W}(x) = \frac{\log E[r(M)]}{\Delta} V(x)$

This provides necessary and sufficient conditions for algorithm design.

4. Variance Analysis and Asymptotic Optimality

Under the subsolution condition, the estimator variance admits a precise asymptotic characterization:

$$\lim_{n \rightarrow \infty} -\frac{1}{n} \log E \left[ (s_{SA}^n)^2 \right] = W(x) - V(x) \frac{\log E\left[ \sum_{i=1}^{r(M)} w_i(M)^2 \right]}{\Delta}$$

As $W(x)$ is the exponential rate governing the rare event probability itself, the second moment optimality is achieved when $\bar{W}(x) = W(x)$, giving a variance decay rate approaching the ideal value $2W(x)$ for unbiased estimators.

Choosing strict subsolutions (i.e., scaling $\bar{W}(x)$ down by a constant $< 1$) trades off a minor reduction in performance for a substantial decrease in the number of particles, which is sometimes necessary to avoid computational bottlenecks.

5. Design Procedure for Asymptotically Optimal Schemes

The optimal design of splitting algorithms follows a systematic procedure:

• Select level spacing $\Delta$ and mean offspring number $u = E[r(M)]$.
• Choose $V(x)$ such that $\bar{W}(x) = \frac{\log u}{\Delta} V(x)$ is a maximal subsolution, i.e., ideally $\bar{W}(x) = W(x)$.
• Implement splitting with constant weights, for instance:

$$\begin{aligned} q_1 &= \lceil u \rceil - u,\quad q_2 = 1 - q_1 \ r(1) &= \lceil u \rceil, \quad r(2) = \lfloor u \rfloor \ w_i(j) &= 1/u \end{aligned}$$

If the number of generated particles is excessive, a strict subsolution can be employed.

This principle formalizes the tradeoff between computational cost and statistical performance and directly links simulation strategy to the large deviation structure of the underlying process.

6. Numerical Examples and Applications

The analysis is supported by numerical experiments across queueing networks (tandem Jackson networks with both shared and separate buffers), Markov- and non-Markov-modulated systems, and even finite-time sample mean problems. In each case, the appropriate subsolution is derived—e.g., for a tandem Jackson network with shared buffer:

$\bar{W}(x) = \langle p, x \rangle + \mathrm{const}$

with $p$ chosen in alignment with large deviation optimal transitions; for the two-buffer case:

$\bar{W}(x) = -\,_1x_1 -\,_2x_2$

Simulation estimates converge to the large deviation predictions, and the growth rate of particle numbers remains subexponential, verifying theoretical claims.

Tables in the referenced work provide detailed comparisons of theoretical and simulated values, including standard errors, confidence intervals, and computational statistics.

7. Broader Impact and Theoretical Implications

Grounding rare event simulation algorithms in large deviation principles and subsolution analysis provides a rigorous basis for algorithm design, ensuring efficiency (in both variance and computational resource usage) and offering a concrete methodology for a diverse array of applications, including non-Markovian systems, queueing theory, statistical physics, and reliability engineering.

The framework unites tools from variational analysis, PDE theory (via HJB equations), and modern simulation methodology, facilitating the cross-fertilization of ideas between large deviations theory and rare event computation. A logical avenue for future research includes extending the subsolution-based methodology to broader classes of stochastic processes, incorporating adaptive schemes for particle control, and exploring further connections to partial differential equations and stochastic control.

The theory underlying rare event simulation by splitting technique establishes strong necessary and sufficient stability conditions, optimality characterizations for estimator variance, and links simulation efficiency directly to the analytic large deviation rate structure, setting a foundation for principled, high-reliability estimation in practical complex systems.