A direct importance sampling-based framework for rare event uncertainty quantification in non-Gaussian spaces

This work introduces a new framework for directly estimating rare event probabilities in complex, high-dimensional non-Gaussian spaces. Rare event uncertainty quantification is a pivotal component in contemporary decision-making across diverse domains. It is particularly crucial for assessing and mitigating risks tied to infrequent yet highly impactful rare event occurrences, such as structural failures induced by extreme loads like earthquakes and hurricanes. This field, commonly referred to as reliability estimation in engineering, has a long history, and comprehensive insights into its complexities can be read in [1, 2, 3, 4, 5, 6]. In many practical applications, rare event probabilities are fortunately exceptionally low, often within the range of $10^{-4}$ to even $10^{-9}$ and lower, posing significant numerical and mathematical challenges. Further complicating these challenges is the reliance on computationally expensive models of the involved systems with high-dimensional random variable spaces. Achieving accurate estimates while minimizing model evaluations is thus a critical parameter.

Numerous existing solution techniques are specifically designed for Gaussian spaces and typically seek to transform non-Gaussian random variables to Gaussian ones. However, when such Gaussian transformations are not possible, these techniques often encounter challenges in non-Gaussian spaces due to potential asymmetries, complex dependencies, intricate geometries, and constrained spaces. Approximation methods, such as First Order Reliability Method (FORM) [7, 8], Second Order Reliability Method (SORM) [9], and Large deviation theory (LDT) [10, 11, 12], provide asymptotic expressions for estimating rare event probabilities [13]. In general, these methods lose accuracy and/or practicality for arbitrary non-Gaussian distributions, and may also exhibit considerable errors in problems with moderately to highly nonlinear performance functions [7, 14]. Sampling methods can tackle the problem in its utmost generality, however, with a higher computational cost. In addressing the known limitations of the crude Monte Carlo approach [15], numerous sophisticated sampling methods have been proposed. Among others, the Subset Simulation (SuS) approach, originally presented in [16], has been widely used and continuously enhanced for Gaussian spaces [17, 18, 19, 20], given its potential to handle high-dimensional problems, despite some limitations [21]. Several more recent works have also shown enhanced performance in non-Gaussian spaces by mainly employing advanced samplers within SuS, such as Hamiltonian MCMC [22, 23], Riemannian Manifold Hamiltonian Monte Carlo (RMHMC) [24, 25], Hamiltonian Neural Networks (HNN) [26], and affine invariant ensemble MCMC sampler [27, 28]. It should be noted that Subset Simulation also has similarities with multilevel splitting methods [29, 30, 31]. Another sampling method capable of efficiently estimating rare event probabilities is importance sampling (IS), a well-known variance-reduction approach. Since this work builds upon the theoretically rigorous foundations of importance sampling, we provide a comprehensive review and discussion of IS techniques for rare event estimation in Section 2. Combining surrogate models, that approximate computationally expensive original models, with sampling-based approaches has also demonstrated significant savings in computational cost across various applications [32, 33, 34]; however, their applicability to high-dimensional, complex problems still faces practical challenges.

In this paper, we introduce a novel framework for efficiently and precisely estimating rare event probabilities, directly in non-Gaussian spaces, building on our foundational Approximate Sampling Target with Post-processing Adjustment (ASTPA) approach, recently developed for Gaussian spaces in [35]. As shown in this work, the developed framework demonstrates excellent performance in challenging high-dimensional non-Gaussian scenarios, thereby overcoming a significant limitation of numerous existing methods. In essence, rare event probability estimation is a normalizing constant estimation problem for the joint distribution of the involved random variables truncated over the rare event domain, also known as the optimal sampling target; a demanding problem traditionally approached through sequential techniques. Conversely, the ASTPA approach directly and innovatively addresses this challenge by decomposing the problem into two less demanding estimation problems. Hence, the proposed methodology first constructs an unnormalized sampling target, appropriately designed for non-Gaussian spaces, relaxing the optimal importance sampling distribution. To sample this constructed target adequately, our Quasi-Newton mass preconditioned Hamiltonian MCMC (QNp-HMCMC) algorithm is employed, particularly suitable for complex spaces, as discussed later. The obtained samples are subsequently utilized not only to compute a shifted estimate of the sought probability but also to guide the second ASTPA stage. Post-sampling, the normalizing constant of the approximate sampling target is computed through the inverse importance sampling (IIS) technique [35], which utilizes an importance sampling density (ISD) fitted based on the samples drawn in the first stage. A practical approach augmenting the IIS stability in challenging non-Gaussian settings, by avoiding anomalous samples, is also devised in this work. The rare event probability is eventually computed by utilizing this computed normalizing constant to correct the shifted estimate of the first stage. The statistical properties of the proposed estimator are thoroughly derived, including its proven unbiasedness and analytical coefficient of variation (C.o.V). Notably, the derived analytical C.o.V expression, accompanied by an applied implementation technique based on the effective sample size, has demonstrated accurate agreement with numerical results.

Hamiltonian Markov Chain Monte Carlo (HMCMC) is one of the most successful and powerful MCMC algorithms [22, 36], based, as the name suggests, on Hamiltonian dynamics [37]. Its performance largely depends on tuning three parameters: the mass matrix, the simulation step size, and the length of the simulated trajectories. Its efficiency in high-dimensional complex spaces can be significantly improved by selecting an appropriate mass matrix [24, 38, 39]. Therefore, in this work, we rigorously analyze a Quasi-Newton mass preconditioned HMCMC (QNp-HMCMC) sampling algorithm, an approach we initially conceptualized in our earlier work [35]. QNp-HMCMC effectively and uniquely leverages the geometric information of the target distribution. During the burn-in stage, a skew-symmetric preconditioning scheme for the Hamiltonian dynamics is adopted [40], incorporating Hessian information merely approximated based on the already available gradients. The estimated Hessian is subsequently utilized to provide an informed, preconditioned mass matrix, better describing the scale and correlation of the involved random variables. To reinforce its theoretical basis, we prove here that simulating the skew-symmetric preconditioned dynamics in QNp-HMCMC leads to the correct stationary target distribution. To effectively implement QNp-HMCMC, we provide detailed recommendations for approximating the Hessian information in complex non-Gaussian spaces. The simulation step size in QNp-HMCMC is automatically tuned using the dual averaging algorithm [41, 42]. Equipped with the preconditioned dynamics and mass matrix, QNp-HMCMC is characterized by large informed sampling steps, enabling it to work efficiently even with single-step simulated Hamiltonian dynamics, avoiding the challenging task of tuning its trajectory length in general, intricate domains, and enabling increased computational efficiency by minimizing the number of model calls (gradient evaluations). As proven in this work, this QNp-HMCMC version with a single-step implementation is equivalent to an original preconditioned Metropolis-adjusted Langevin algorithm (MALA). Finally, an optimization approach is also devised in this paper to initiate the QNp-HMCMC sampling effectively.

This paper is structured as follows. Section 2 offers a comprehensive discussion on existing importance sampling-based methods for rare event probability estimation, highlighting their common characteristics and, notably, their connections and differences with our ASTPA approach. The developed framework is subsequently presented in Section 3, and Section 4 discusses in detail the single-step QNp-HMCMC algorithm. Section 5 introduces an effective optimization approach for rare event domain discovery, specifically using the Adam optimizer [43], aiming to provide good initial states for the MCMC sampler. Given that it is more practical for MCMC and HMCMC algorithms to operate in unconstrained spaces, Section 6 discusses transforming bounded random variables to unconstrained ones. A summary of the QNp-HMCMC-based ASTPA framework for non-Gaussian domains is then presented in Section 7. To fully examine the capabilities of the proposed methodology, in Section 8 its performance is demonstrated and successfully compared against the state-of-the-art Subset Simulation, and related variants, in a series of challenging low- and high-dimensional non-Gaussian problems. The paper then concludes in Section 9.

Let $\bm{X}$ be a continuous random vector taking values in $\mathcal{X}\subseteq\mathbb{R}^{d}$ and having a joint probability density function (PDF) $\pi_{\bm{X}}$. We are interested in rare events $\mathcal{F}\coloneqq\{\bm{x}\,:\,g(\bm{x})\leq 0\}$, where $g:\mathcal{X}\rightarrow\mathbb{R}$ is a performance function, also known as the limit-state function, defining the occurrence of a rare event. In the context of reliability estimation, which considers failures as rare events, the limit state function, for example, can be formulated as $g(\bm{X})=\lambda-\mathcal{M}(\bm{X})$, where $\mathcal{M}(\bm{X})$ denotes a model output (a quantity of interest). Here, the predefined threshold $\lambda$ determines the rarity of the failure event by controlling when $g(\bm{X})$ starts to take non-positive values, i.e., rare event occurrence. In this work, we aim to estimate the rare event probability $p_{\mathcal{F}}$:

where $I_{\mathcal{F}}:\mathcal{X}\rightarrow\{0,1\}$ is the indicator function, i.e., $I_{\mathcal{F}}(\bm{x})=1$ if $\bm{x}\in\mathcal{F}$, and $I_{\mathcal{F}}(\bm{x})=0$ otherwise, and $\mathop{\mathbb{E}}$ is the expectation operator.

Our objective in this work is to estimate the described integration in Eq. 1 under these challenging yet realistic settings: (i) the analytical calculation of Eq. 1 is generally intractable, as the limit-state function can be complicated, e.g., involve the solution of a differential equation; (ii) the computational effort for evaluating $I_{\mathcal{F}}$ for each realization $\bm{x}$ is assumed to be quite significant, often relying on computationally intensive models, necessitating the minimization of such function evaluations (model calls); (iii) the rare event probability $p_{\mathcal{F}}$ is extremely low, e.g., in order of $p_{\mathcal{F}}\sim 10^{-4}-10^{-9}$; (iv) the random variable space, $\mathcal{X}\subseteq\mathbb{R}^{d}$, is assumed to be high-dimensional, e.g., $d=10^{2}$ and more; (v) the joint probability density, $\pi_{\bm{X}}(\bm{x})$, is non-Gaussian, with transformations to the Gaussian domain assumed inapplicable, necessitating Eq. 1 evaluation directly in non-Gaussian spaces.

Under these settings, several sampling methods become highly inefficient and fail to address the problem effectively. For instance, the standard Monte Carlo estimator of Eq. 1 using $\{\bm{x}_{i}\}_{1}^{N}$ draws from $\pi_{\bm{X}}$, $\hat{p}_{\mathcal{F}}=\frac{1}{N}\sum_{i=1}^{N}I_{\mathcal{F}}(\bm{x}_{i})$, has a coefficient of variation $\sqrt{(1-p_{\mathcal{F}})/(Np_{\mathcal{F}})}$, rendering this estimate inefficient by requiring a prohibitively large number of samples to accurately quantify small probabilities; see setting (iii) above. Importance sampling (IS), as a variance-reduction method, can efficiently compute the integral in Eq. 1, through sampling from an appropriate importance sampling density (ISD), $h(\bm{\bm{X}})$. The ISD places greater importance on rare event domains in the random variable space compared to $\pi_{\bm{X}}(\bm{x})$, satisfying the sufficient condition $\text{supp}(I_{\mathcal{F}}\,\pi_{\bm{X}})\subseteq\text{supp}(h)$, with $\text{supp}(\cdot)$ denoting the support of the function. This results in an increased number of rare event samples, thereby reducing the number of required model calls. Eq. (1) can then be written as:

leading to the unbiased IS estimator:

The efficiency of IS relies heavily on the careful selection of both the ISD and the sampling algorithm. A theoretically optimal ISD $h^{*}$, providing zero-variance IS estimator, can be given as [44]:

Yet, it is evident that this ISD cannot be utilized due to the fact that the normalizing constant in Eq. (4) is the sought target probability $p_{\mathcal{F}}$. Notably, finding a near-optimal ISD, approximating $h^{*}$, is an ongoing research pursuit. The existing approaches can be largely classified into two categories. The first directly follows the IS estimator in Eq. 3 and seeks to approximate $h^{*}$ by a parametric or non-parametric probability density function (PDF), $h$. An early approach in this category utilizes densities centered around one or more design points [45], provided by an initial FORM analysis, which has later been deemed ineffective for high-dimensional, complex spaces [16]. A more recent popular approach is the cross entropy-based importance sampling (CE-IS). It alternatively constructs a parametric ISD $h$ by minimizing the Kullback-Leibler (KL) divergence between $h^{*}$ and a chosen parametric family of probability distributions [46, 47, 48, 49, 50]. Recent CE-IS variants incorporate information from influential rare event points, e.g., design points, to effectively initialize the cross entropy process [51, 52]. Within this first category of methods, several different noteworthy approaches have also been proposed to approximate $h^{*}$, e.g., [53, 54, 55, 56, 57].

Before delving into the second category of methods, we define ISD $h={\nicefrac{{\tilde{h}(\bm{X})}}{{C_{h}}}}$, where $\tilde{h}$ is an unnormalized importance distribution, and $C_{h}=\int_{\mathcal{X}}\tilde{h}(\bm{x})d\bm{x}$ is its normalizing constant. Following this definition, $\tilde{h}^{*}=I_{\mathcal{F}}(\bm{X})\pi_{\bm{X}}(\bm{X})$ denotes an unnormalized optimal importance sampling distribution, and $C_{h^{*}}=p_{\mathcal{F}}$ is its normalizing constant. Thus, in the second category of methods, sequential variants of importance sampling are employed, mainly utilizing a number of intermediate unnormalized distributions. These unnormalized distributions are generally more flexible than the PDFs usually utilized within the first category, as discussed above, thus being capable of providing a better approximation of $h^{*}$. Employing an unnormalized importance sampling distribution $\tilde{h}$, Eq. 2 can be written as:

Existing methods currently compute the normalizing constant, $C_{h}$, using consecutive intermediate unnormalized distributions $\tilde{h}_{i}$’s, $i=1,2,...,n$, following Eq. 6. These distributions should be designed to link between $\pi_{\bm{X}}$ and $h^{*}$, with $\tilde{h}_{1}$ and $\tilde{h}_{n}=\tilde{h}$ being the closest to $\pi_{\bm{X}}$ and $h^{*}$, respectively.

where $C_{h_{i}}$ is the normalizing constant of $h_{i}$, and $C_{\pi}=1.0$, typically, is the normalizing constant of the original distribution $\pi_{\bm{X}}$. $C_{\pi}$ is left in Eq. 6 for generality purposes, particularly to cover scenarios when an unnormalized original distribution $\tilde{\pi}_{\bm{X}}$ is instead provided (substitute $\pi_{\bm{X}}={\nicefrac{{\tilde{\pi}_{\bm{X}}}}{{C_{\pi}}}}$ in Eq. 5 to verify that). Eq. 6 can be computed using the sequential importance sampling (SIS) approach as [58, 59]:

Alternatively, bridge sampling can evaluate Eq. 6 as [60, 61, 62]:

where $h_{i}^{B}(\bm{X})$ is a bridge distribution [61] between $\tilde{h}_{i}$ and $\tilde{h}_{i-1}$, and $C_{h_{i}^{B}}$ is its associated normalizing constant. This latter technique aims to provide a more accurate estimate even if the two consecutive distributions, $\tilde{h}_{i}$ and $\tilde{h}_{i-1}$, are not close enough. Eq. 6 can also be computed using path sampling [63, 64]. A key to the success of these sequential methods is to design nearby intermediate distributions and to sufficiently sample them to accurately compute Eq. 6, a non-trivial and computationally demanding task. A discussion on methods within this second category can be found in [65]. Interestingly, the subset simulation method [16] can also be seen as a special case of this second category, wherein the ratio ${\nicefrac{{C_{h_{i}}}}{{C_{h_{i-1}}}}}$ is the conditional probability between every two consecutive intermediate events, $\mathcal{F}_{i}$ and $\mathcal{F}_{i-1}$.

In this work, we are originally estimating Eq. 5 in a non-sequential manner, using a single unnormalized distribution $\tilde{h}$, approximating $h^{*}$, i.e., $n=1$. Therefore, the proposed approach significantly reduces the computational cost associated with employing multiple intermediate distributions. The discovery of the rare event domain is alternatively achieved through the utilization of advanced samplers and optimization techniques, as detailed in Section 5. It should also be noted that the discussed importance sampling-based methods in the current literature are mostly tailored for Gaussian spaces. Consequently, the ability of our proposed framework to operate directly in non-Gaussian spaces is an additional substantial benefit.

ASTPA estimates the rare event probability in Eq. 1 through an innovative non-sequential importance sampling variant evaluating Eq. 5. To this end, ASTPA constructs a single unnormalized approximate sampling target $\tilde{h}$, relaxing the optimal ISD in Eq. 4. Post-sampling, its normalizing constant is computed utilizing our devised inverse importance sampling. For clarity, we rewrite Eq. 5 as:

Computing $p_{\mathcal{F}}$ is thus decomposed into two (hopefully simpler) problems. The first involves constructing $\tilde{h}$ and sampling $\{\bm{x}_{i}\}_{i=1}^{N}\sim h$, to compute the unbiased expectation of the weighted indicator function (shifted probability estimate) as:

An empirical estimation for the normalizing constant, $\hat{C}_{h}$, is then computed using the inverse importance sampling, as discussed in Section 3.2. The ASTPA estimator for the sought rare event probability can then be computed as:

Fig. 1 concisely portrays the proposed framework using a correlated Gumbel distribution and a quadratic limit state function, as detailed in Section 8.1, with $p_{\mathcal{F}}\sim 2.51\times 10^{-7}$. The construction of the approximate sampling target in ASTPA is discussed in Section 3.1, including its recommended parameters for problems described in non-Gaussian spaces. Then, we present the inverse importance sampling approach in Section 3.2. Subsequently, Section 3.3 discusses thoroughly the statistical properties of the proposed estimator presented in Eq. 11.

In this section, we discuss sampling the approximate target $\tilde{h}$ in ASTPA using gradient-based Markov Chain Monte Carlo (MCMC) methods. Specifically, we focus on the state-of-the-art Hamiltonian MCMC (HMCMC) methods, including our developed Quasi-Newton mass preconditioned HMCMC (QNp-HMCMC), aimed at improving sampling efficiency for intricate probability distributions. We first present the standard HMCMC algorithm, followed by our developed QNp-HMCMC, emphasizing its theoretical foundations and practical implementation techniques for complex non-Gaussian spaces. Connections with the preconditioned Metropolis-adjusted Langevin algorithm (MALA) are illustrated, highlighting the uniqueness of QNp-HMCMC, even with a single-step integration.