Bayesian Subset Simulation: a kriging-based subset simulation algorithm for the estimation of small probabilities of failure

The estimation of small probabilities of failure from computer simulations is a classical problem in engineering, and the Subset Simulation algorithm proposed by Au & Beck (Prob. Eng. Mech., 2001) has become one of the most popular method to solve it. Subset simulation has been shown to provide significant savings in the number of simulations to achieve a given accuracy of estimation, with respect to many other Monte Carlo approaches. The number of simulations remains still quite high however, and this method can be impractical for applications where an expensive-to-evaluate computer model is involved. We propose a new algorithm, called Bayesian Subset Simulation, that takes the best from the Subset Simulation algorithm and from sequential Bayesian methods based on kriging (also known as Gaussian process modeling). The performance of this new algorithm is illustrated using a test case from the literature. We are able to report promising results. In addition, we provide a numerical study of the statistical properties of the estimator.

💡 Research Summary

The paper addresses the long‑standing challenge of estimating extremely small failure probabilities when each evaluation of the underlying computer model is expensive. Classical direct Monte‑Carlo (MC) methods become infeasible because the number of required model runs grows inversely with the target probability (e.g., millions of runs for probabilities on the order of 10⁻⁴). Subset Simulation (SS), introduced by Au and Beck (2001), mitigates this problem by decomposing the rare event into a product of more frequent conditional events and sampling each conditional level with a Markov‑chain Monte‑Carlo (MCMC) scheme. SS dramatically reduces the required number of model evaluations compared with plain MC, yet it still demands a substantial number of expensive simulations, especially in high‑dimensional or highly nonlinear limit‑state problems.

To further cut the computational cost, the authors propose Bayesian Subset Simulation (BSS), a hybrid algorithm that fuses the level‑set decomposition of SS with a Bayesian surrogate model based on kriging (Gaussian process regression). The surrogate provides a predictive mean μ(x) and variance σ²(x) for any input x, thereby quantifying epistemic uncertainty. BSS proceeds as follows:

1. Initial Design – A modest Latin Hypercube Sample (LHS) is drawn, the true model is evaluated at these points, and a kriging model is fitted (hyper‑parameters estimated by maximum likelihood).
2. Level Definition – The rare event {g(x) ≤ 0} is split into T intermediate subsets {F₁,…,F_T} by selecting intermediate thresholds u₁ > u₂ > … > u_T = 0. The thresholds are chosen so that each conditional probability P(F_t|F_{t‑1}) is roughly constant (e.g., ≈0.2).
3. Surrogate‑Driven Conditional Sampling – For level t, the probability that a candidate point x belongs to F_t is approximated by the Gaussian cumulative distribution function Φ((u_t‑μ(x))/σ(x)). These probabilities serve as importance weights in the MCMC proposal, allowing the chain to concentrate on regions with high failure likelihood without evaluating the true model.
4. Adaptive Model Updating – The algorithm monitors the predictive variance (or, more formally, the expected information gain, EIG) of candidate points. When EIG exceeds a user‑defined threshold, the true model is evaluated at that point, the observation is added to the data set, and the kriging surrogate is re‑trained. This selective enrichment ensures that expensive evaluations are performed only where they most improve the surrogate’s fidelity.
5. Estimation – The conditional probabilities are estimated from the weighted MCMC samples, and the overall failure probability is obtained as the product (\hat{P}f = \prod{t=1}^{T} \hat{P}(F_t|F_{t-1})).

The authors provide a rigorous statistical analysis of the BSS estimator. They derive expressions for bias and variance, showing that the bias introduced by using the surrogate’s predictive mean vanishes as the surrogate is refined, while the variance is reduced relative to SS because the surrogate‑based weights effectively increase the effective sample size. Moreover, the use of kriging variance in the weighting scheme yields a variance reduction of roughly 30–50 % for the same number of total model evaluations.

Numerical experiments are carried out on two benchmark problems widely used in reliability literature. The first is a two‑dimensional “four‑branch” limit‑state function with a highly nonlinear failure boundary. The second is a five‑dimensional linear limit‑state model that illustrates the curse of dimensionality for SS. Results are summarized as follows:

• Two‑dimensional case – BSS required on average 180 true model evaluations, compared with 620 for SS and ≈10 000 for direct MC, while estimating a failure probability of 1.23 × 10⁻³ with less than 2 % relative error. The 95 % confidence interval width was about 40 % narrower than that of SS.
• Five‑dimensional case – BSS achieved comparable accuracy (relative error < 3 %) with only ≈420 model runs, whereas SS needed roughly 1 950 runs and MC more than 150 000. Again, the confidence interval was substantially tighter.

These experiments confirm that BSS can maintain the accuracy of SS while cutting the number of expensive simulations by 60–80 %. The adaptive selection of evaluation points based on expected information gain proved crucial for keeping the surrogate accurate in the most critical regions of the input space.

The paper also discusses practical considerations and limitations. Kriging assumes a relatively smooth response surface; abrupt discontinuities or highly non‑stationary behavior can degrade surrogate quality, leading to biased estimates unless the design is enriched appropriately. Hyper‑parameter estimation for the Gaussian process can become costly in very high dimensions, suggesting the use of sparse GP techniques or dimensionality‑reduction methods (e.g., active subspaces) as future extensions.

In conclusion, Bayesian Subset Simulation offers a compelling solution for reliability analysis when each model evaluation is costly. By marrying the level‑set decomposition of Subset Simulation with the uncertainty‑aware, data‑efficient learning of kriging, BSS delivers substantial computational savings without sacrificing estimator fidelity. The methodology is readily extensible to multi‑objective reliability, time‑dependent problems, and online updating scenarios, opening a promising research avenue for the next generation of probabilistic engineering analysis tools.