Hamiltonian MCMC Methods for Estimating Rare Events Probabilities in High-Dimensional Problems

In this work, we develop a framework for estimation of rare events probabilities, a commonly encountered important problem in several engineering and scientific applications, often observed in the form of probability of failure ($P_{F}$) estimation or, alternatively, reliability estimation. In many practical applications, failure probabilities are fortunately very low, from $10^{-4}$ to even $10^{-9}$ and lower, and calculating such small probabilities presents many numerical and mathematical challenges, particularly in cases with high dimensional random spaces and/or expensive computational models, that practically limit the afforded number of model calls. The number of model calls is thus of great importance in these problems and one of the critical parameters that limits or prohibits use of several available techniques in the literature.

The reliability estimation problem has a long history in the engineering community [1, 2, 3, 4, 5]. One of the significant early successes was the discovery of the so called First Order Reliability Method (FORM) [6, 7], long investigated by Der Kiureghian, Ditlevsen and co-workers [8, 9], and many others, e.g., Shinozuka [10], providing also several enhancements, including second order effects (SORM) by Breitung [11]. In FORM/SORM methods, the search for the most probable failure point (MPP) is usually performed by gradient-based optimization methods [12, 13]. Although these asymptotic approximation methods are usually computationally inexpensive, they have several limitations and may involve considerable errors, particularly in high-dimensional problems or in problems with highly nonlinear limit-state functions [6, 14]. As such, various sampling-based methods have also been suggested by the reliability community, e.g., Schuëller and Pradlwarter [15], to tackle the problem in its utmost generality, with crude Monte Carlo approaches being prohibitive for this type of problems due to their excessive computational demands. Only some of many notable contributions can be seen in [16, 17, 18, 19, 20], describing and studying the state-of-the art-Subset Simulation and its enhancements, originally presented in [21], and in [22, 23, 24] utilizing importance sampling schemes, often also combined with the cross-entropy method [25, 26, 27, 28, 29]. Alternative approaches include directional and line Sampling [30, 31, 32, 33], the PHI2 method for time-variant reliability [34], and asymptotic sampling strategies [35, 36], among others. The problem of estimating rare event probabilities has also attracted a great deal of attention in other relevant communities and in mathematical literature, with several suggested methods sharing many similarities with FORM/SORM approaches, e.g.,[37], and Subset Simulation, such as in approaches involving Sequential Monte Carlo samplers for rare events [38, 39, 40], interacting particle system methodologies [41, 42], multilevel splitting methods [43, 44] and forward flux sampling [45], to name but a few.

In this paper, we are presenting a new solution approach to the problem by combining gradient-based approaches, already familiar to the engineering community and often available in computational tools, with Markov Chain Monte Carlo (MCMC) sampling methods in the form of Hamiltonian MCMC. MCMC methods [46] are plausibly the most broadly accepted ones to generate samples from target distributions, in cases where direct sampling is not possible. Despite notable successes, many MCMC methods scale poorly with the number of dimensions and can become inefficient. For complicated multivariate models, classic methods such as random-walk Metropolis-Hastings [47] and Gibbs sampling [48] may require an unacceptably long time and number of samples to adequately explore the target distribution. Originally developed by Duane et al. [49] and to a large extent understood and popularized through the works of Neal [50, 51], Hamiltonian Markov Chain Monte Carlo (HMCMC), usually called Hamiltonian Monte Carlo (HMC) in the literature, produces Markov chain samples based on Hamiltonian dynamics principles, is characterized by much better scalability [51, 52] and faster mixing rates, is capable of generating samples with much weaker auto-correlation, even in complex high-dimensional random spaces [53], and has enjoyed broad-spectrum successes in most general settings [54]. Balanced against these features and achievements is of course the need for multiple gradient evaluations in each HMCMC iteration, making the method more computationally intensive per iteration than other algorithms, such as random-walk Metropolis-Hastings and Metropolis-adjusted Langevin [55, 56], for example. Girolami and Calderhead introduced a Riemannian Manifold Hamiltonian Monte Carlo (RMHMC) approach in [57, 58] that has demonstrated significant successes in many challenging problems but requires computing higher-order derivatives of the target distribution. Overall, the two impediments to using Hamiltonian MCMC methods are the required gradients, since analytical formulas are not always available and numerical techniques are computationally costly, particularly in high dimensions, and the heedful tuning of the involved parameters [51]. The first issue can in certain cases be solved by automatic differentiation (e.g. [54, 59]) and stochastic gradient approaches [60], while for the second a fully automated state-of-the-art HMCMC algorithm has been developed by Hoffman and Gelman, known as the No-U-Turn Sampler (NUTS [61]). NUTS introduces, among others, an expensive tree building procedure, in order to trace when the Hamiltonian trajectory turns back on itself. Many of these approaches are not however relevant and/or applicable to the analyzed problem in this work, since, in general, many rare event and reliability problems involve complex, computationally expensive models, complicating and/or precluding use of automatic differentiation and data-based stochastic gradient techniques, as well as methodologies that require a considerably high number of model calls per problem, such as NUTS.

A new computationally efficient sampling framework for estimation of rare events probabilities is thus presented in this work, having exceptional performance in quantifying low failure probabilities for any type of reliability problems described in both low and high dimensional stochastic spaces. The introduced methodology is termed Approximate Sampling Target with Post-processing Adjustment (ASTPA) and comprises a sampling and a post-processing phase. The sampling target in ASTPA is constructed by appropriately combining the multi-dimensional random variable space with a cumulative distribution function that utilizes the limit-state function. Having acquired the samples, an adjustment step is then applied, in order to account for the fact that the samples are drawn from an approximate target distribution, and to thus correctly quantify the rare event probability. An original inverse importance sampling procedure is devised for this adjustment step, taking its name from the fact that the samples are already available. Although the ASTPA framework is general and can be combined with any appropriate Monte Carlo sampling method, it becomes substantially efficient when directly supported by gradient-based Hamiltonian Markov Chain Monte Carlo (HMCMC) samplers. To address the scalability issues a typical HMCMC sampler may manifest in high-dimensional spaces, a new Quasi-Newton mass preconditioned HMCMC approach is also developed. This new sampling scheme follows an approximate Newton direction and estimates the pertinent Hessian in its burn-in stage, only based on gradient information and the BFGS approximation, and eventually utilizes the computed Hessian as a preconditioned mass matrix in the main non-adaptive sampling stage. An approximate analytical expression for the uncertainty of the computed estimation is also derived, showcasing significant accuracy with numerical results, and all involved user-defined parameters of ASTPA are thoroughly analyzed and general default values are suggested. Finally, to fully examine the capabilities of the proposed methodology, its performance is demonstrated and compared against Subset Simulation in a series of challenging low- and high-dimensional problems.

The failure probability P_F for a system, that is the probability of a defined unacceptable system performance, can be expressed as a $d$-fold integral, as:

where $\theta$ is the random vector $[\theta_{1},...,\theta_{d}]$^T, $F\subset\mathbb{R}^{d}$ is the failure event, g($\theta$) is the limit-state function that can include one or several distinct failure modes and defines the system failure by g($\theta$)$\leq$ 0, I(.) denotes the indicator function with $I_{F}$ ($\theta$) = 1 if $\theta$ $\in$ g($\theta$)$\leq$ 0 and $I_{F}$($\theta$) = 0 otherwise, $\mathop{\mathbb{E}}$ is the expectation operator, and $\pi_{\Theta}$ is the joint probability density function (PDF) for $\bm{\Theta}$. As is common practice for problems of this type, in this work the joint PDF of $\bm{\Theta}$ is the standard normal one, due to its rotational symmetry and exponential probability decay. In most cases this is not restrictive, since it is uncomplicated to transform the original random variables X to $\bm{\Theta}$, e.g. [62]. When this is not the case, however, and the probabilistic characterization of X can be defined in terms of marginal distributions and correlations, the Nataf distribution (equivalent to Gaussian copula) is commonly used to model the joint PDF, and the mapping to the standard normal space can be then accomplished [8, 63].

The focus in this work is to analyze the described integration in Eq.˜1 under very general settings, including the following challenging sampling context: (i) Computation of Eq.˜1 can only be done in approximate ways; (ii) the relationship between $\theta$ and $I_{F}$ is not explicitly known and for any $\theta$ we can merely check whether it belongs to the failure set or not, i.e. calculate the value $I_{F}$($\theta$); (iii) the computational effort for evaluating $I_{F}$($\theta$) for each value of $\theta$ is assumed to be quite significant, so that it is essential to minimize the number of such function evaluations (model calls); (iv) the probability of failure P_F is assumed to be very small, e.g. in order of $\textit{P\textsubscript{F}}\sim 10^{-4}-10^{-9}$; (v) the parameter space $\mathbb{R}^{d}$ is assumed to be high-dimensional, in the order of $10^{2}$ and more, for example. Under these general settings, several sampling methods, including direct Monte Carlo approaches and NUTS [61], become highly inefficient and fail to address the problem effectively. Subset Simulation (SuS) [21] has however proven successful and robust in dealing with problems of this type and is shown to outperform other relevant methods in numerous papers, e.g. [24, 31]. SuS relies on a modified component-wise Metropolis MCMC method that can successfully work in high dimensions and does not require a burn-in sampling stage. A notable adaptive conditional sampling (aCS) methodology within the SuS framework is also introduced by Papaioannou et al. in [16], providing important advantages and enhanced performance in several cases. Relevant SuS variants are thus utilized in this work, for validation and comparison purposes with our presented methodology that completely deviates from SuS and is efficiently supported by a direct Hamiltonian MCMC sampling approach [64].

In complex high-dimensional problems, the performance of the typical HMCMC sampler, presented as Algorithm˜1, may deteriorate and a prohibitive number of model calls could be required. A variety of methods have been proposed in the literature to address this issue. Among others, a Riemannian Manifold Hamiltonian Monte Carlo (RMHMC) has been suggested in [57] that takes advantage of the manifold structure of the variable space, at the cost of calculating second- and third-order derivatives of distributions and using a generalized leapfrog scheme, requiring additional model calls per leapfrog step. Although possible in some cases regarding second-order derivatives, e.g. [67], in the majority of cases higher-order derivatives are not provided by computational models, such as finite element models. In addition, the computational cost still increases importantly and extra model calls per leapfrog step are usually restrictive for computationally expensive models.

In this work, we instead address the high-dimensionality performance issue in a different, Newton-type context, without needing additional model calls per leapfrog step and by only using the Hessian information of the target distribution, either the exact one, when relevant information is provided freely by the computational model, or, even more general, an approximate one that does not increase the computational cost. An approximate Hessian can be given in a systematic manner based on already available gradient information, similar to Quasi-Newton methods used in nonlinear programming [68]. The well-known BFGS approximation [68] is thus utilized for our Quasi-Newton type Hamiltonian MCMC approach and all numerical examples in this work are analyzed accordingly, based solely on this most general approximate Hessian case. Let $\bm{\theta}\in\mathbb{R}^{d}$, consistent with the previous sections. Given the $k$-$th$ estimate $\textbf{W}_{k}$, where $\textbf{W}_{k}$ is an approximation to the inverse Hessian at $\bm{\theta}_{k}$, the BFGS update $\textbf{W}_{k+1}$ can be expressed as:

where I is the identity matrix, $\textbf{{s}}_{k}=\bm{\theta}_{k+1}-\bm{\theta}_{k}$, and $\textbf{{y}}_{k}=-\nabla\mathcal{L}(\bm{\theta}_{k+1})+\nabla\mathcal{L}(\bm{\theta}_{k})$ where $\mathcal{L}:\mathbb{R}^{d}\longrightarrow\mathbb{R}$ denotes the log density of the target distribution, as before. There is a long history of efficient BFGS updates for very large systems and several numerical techniques can be used, including sparse and limited-memory approaches.

Two relevant studies in the literature on Quasi-Newton extensions and connections to MCMC algorithms can be found in [69, 70]. Our developed method, however, has fundamental differences that are summarized in that we are focusing on Hamiltonian methods, we are using two integrated coupled phases, an adaptive and a non-adaptive, and finally we consistently incorporate the Quasi-Newton outcomes in both stages of momentum sampling and simulation of Hamiltonian dynamics. In more detail, in the adaptive burn-in phase of the algorithm we are still sampling the momentum from an identity mass matrix, $\textbf{M}=\textbf{I}$, but the ODEs of Eq.˜3 now become:

which is equivalent to the implicit linear transformation $\bm{\theta}^{\prime}=\textbf{W}\bm{\theta}$, and W is given by Eq.˜5. Accordingly, the leapfrog integrator is then reformulated as:

Hence, this new dynamic scheme efficiently and compatibly adjusts both the z and $\bm{\theta}$ evolutions based on the local structure of the target distribution, and also features a Quasi-Newton direction for the momentum variables, allowing large jumps across the state space. The final estimation of the approximated inverse Hessian matrix, W, in the adaptive burn-in phase is then used in the subsequent non-adaptive phase of the algorithm as a preconditioned mass (inverse covariance) matrix, $\textbf{M}=\textbf{W}^{-1}$, used to sample the Gaussian momentum $\textbf{z}\sim\textbf{N}(\textbf{0},\textbf{M})$. As such, typical Hamiltonian dynamics are now used, albeit with this properly constructed mass matrix that takes into account the scale and correlations of the position variables, leading to significant efficiency gains, particularly in high-dimensional problems. The BFGS procedure in Eq.˜5 normally provides a symmetric, positive-definite W matrix in an optimization context. However, in our case we are using BFGS under different settings that may not satisfy the curvature condition $\textbf{{s}}_{k}^{T}\textbf{{y}}_{k}>0$, resulting in occasional deviations from positive-definiteness. Several standard techniques can be then implemented to ensure positive-definiteness, such as a damped BFGS updating [68] or the simple addition $\textbf{W}_{new}=\textbf{W}_{old}+\delta\textbf{I}$, where $\delta\geq 0$ is some appropriate number. A straightforward method to determine $\delta$ is to choose it larger than the absolute value of the minimum eigenvalue of $\textbf{W}_{old}$. Another technique involves utilizing a semidefinite programming approach to identify an optimized diagonal matrix to add to $\textbf{W}_{old}$. Alternatively, W can be updated only when the curvature condition is satisfied, which directly guarantees positive definiteness. To further ensure the stability of the sampler, a positive threshold can be introduced to the curvature condition instead of zero, e.g., $\textbf{{s}}_{k}^{T}\textbf{{y}}_{k}>10$. This latter approach has been used and has worked well in this work. Since the final estimation of W in the adaptive burn-in phase is then extensively utilized in the subsequent non-adaptive phase, we suggest use of a directly provided positive-definite matrix W at this step. This can simply be accomplished by adding one more burn-in iteration step at the end of the burn-in phase, until an appropriate sample, directly supported by a positive definite W matrix is identified.

Our derived Quasi-Newton mass preconditioned Hamiltonian Markov Chain Monte Carlo (QNp-HMCMC) method is concisely summarized and presented in Algorithm˜2. Overall, QNp-HMCMC is a practical, efficient approach that only requires already available gradient information and provides important insight about the geometry of the target distribution, eventually improving computational performance and enabling faster mixing.

In order for an appropriate number of samples to discover and enter the relevant regions, contributing to the rare event probability estimation, a suitable approximate target distribution is constructed in this work, as analyzed in Section˜5.1, then sampled by Hamiltonian MCMC methods that can effectively reach regions of interest. For their initial stage, our HMCMC samplers have an adaptive annealed phase. This adaptive phase will be thoroughly explained later in Section˜5.3. To estimate the pertinent probability, Eq.˜1 needs to be then adjusted accordingly, since the samples are sampled from our constructed target distribution and not $\pi_{\Theta}(\bm{\theta})$. An original post-sampling step is devised at this stage, using an inverse importance sampling procedure, described in Section˜5.2, i.e. first sample, then choose the importance sampling density automatically, based on the samples. Fig.˜1 concisely portrays the overall approach by using a bimodal target distribution with $P_{F}\sim 3.95\times 10^{-5}$. The gray curves represent the parabolic limit-state function $g(\bm{\theta})$ of this problem, with the failure domain being outside $g(\bm{\theta})$. The left figure displays the constructed target distribution, which in this simple 2D case can be visualized. The middle figure shows drawn samples from the target distribution by our suggested QNp-HMCMC algorithm, described in Section˜4, and the right figure demonstrates the inverse importance sampling step. All these different steps will be discussed in detail in the following sections.

Several numerical examples are studied in this section to examine the performance and efficiency of the proposed methods. In all examples, input parameters follow the provided guidelines in Section˜5.5. To compute the normalizing constant $C_{h}$, IIS samples around $20\%$ of the total number of samples have been drawn from a computed GMM with diagonal covariance matrices and, generally, $10$ and $1$ Gaussian components for low and high dimensional problems, respectively. Results are compared with the Component-Wise Metropolis-Hastings based Subset Simulation (CWMH-SuS) [21], with two proposal distributions, a uniform one of width 2 and a standard normal one. The adaptive Conditional Sampling (aCS) SuS variant introduced in [16] is also used in all examples, implemented based on the online provided code by the authors [72]. In all examples, the SuS parameters are chosen as $n_{s}$ = 1,000 and 2,000 for low- and high-dimensional problems, when needed, respectively, with $n_{s}$ the number of samples in each subset level, and $p_{0}=0.1$, where $p_{0}$ is the samples percentile for determining the appropriate subsets. Comparisons are provided in terms of accuracy and computational cost, reporting the P_F estimation and the mean number of limit-state function calls. Analytical gradients are provided in all examples, hence one limit-state/model call can provide both the relevant function values and gradients. The number of limit-state function evaluations for the HMCMC-based methods has been determined based on reported C.o.V values $\in[0.1,\,0.40]$, as estimated by 500 independent simulations. C.o.V values estimated by the analytical expression in Eq.˜26 are also reported in parenthesis, for both HMCMC and QNp-HMCMC methods. The reference failure probabilities are obtained based on the mean estimation of 100 independent simulations using $10^{8}$ crude Monte Carlo samples, where applicable for larger reference probabilities, and SuS method with $n_{s}$ = 100,000 for smaller reference probabilities. The problem dimensions are denoted by $d$ and all ASTPA parameters are carefully chosen for all examples but are not optimized for any one. Comparative and perhaps improved alternate performance might thus be achieved with a different set of parameters.