A direct Hamiltonian MCMC approach for reliability estimation
AA D
IRECT H AMILTONIAN
MCMC A
PPROACH FOR R ELIABILITY E STIMATION
A P
REPRINT
Hamed Nikbakht
Department of Civil and Environmental EngineeringThe Pennsylvania State UniversityUniversity Park, PA 16802 [email protected]
Konstantinos G. Papakonstantinou
Department of Civil and Environmental EngineeringThe Pennsylvania State UniversityUniversity Park, PA 16802September 11, 2019 A BSTRACT
Accurate and efficient estimation of rare events probabilities is of significant importance, since oftenthe occurrences of such events have widespread impacts. The focus in this work is on preciselyquantifying these probabilities, often encountered in reliability analysis of complex engineeringsystems, by introducing a gradient-based Hamiltonian Markov Chain Monte Carlo (HMCMC)framework, termed Approximate Sampling Target with Post-processing Adjustment (ASTPA). Thebasic idea is to construct a relevant target distribution by weighting the high-dimensional randomvariable space through a one-dimensional likelihood model, using the limit-state function. To samplefrom this target distribution we utilize HMCMC algorithms that produce Markov chain samplesbased on Hamiltonian dynamics rather than random walks. We compare the performance of typicalHMCMC scheme with our newly developed Quasi-Newton based mass preconditioned HMCMCalgorithm that can sample very adeptly, particularly in difficult cases with high-dimensionality andvery small failure probabilities. To eventually compute the probability of interest, an original post-sampling step is devised at this stage, using an inverse importance sampling procedure based on thesamples. The involved user-defined parameters of ASTPA are then discussed and general defaultvalues are suggested. Finally, the performance of the proposed methodology is examined in detailand compared against Subset Simulation in a series of static and dynamic low- and high-dimensionalbenchmark problems. K eywords Hamiltonian MCMC · Quasi-Newton · Rare Event Probability · High-dimensional Parameter Space · Reliability Estimation.
In this work, we investigate Hamiltonian Markov Chain Monte Carlo (HMCMC) schemes for estimation of rare eventsprobabilities, a commonly encountered important problem in several engineering and scientific applications [1, 2],most often observed in the form of failure probability, or alternatively, reliability estimation. Calculating such smallprobabilities with accuracy presents many numerical and mathematical challenges, particularly in cases with highdimensional random spaces and/or expensive computational models, that practically limit the afforded number of modelcalls. The well known gradient based First Order Reliability Method (FORM), and variants, have a very long history inreliability estimation problems, with numerous successes [3, 4, 5, 6]. Such asymptotic approximation methods naturallyhave of course limitations, however, in general settings. Hence, numerous sampling based methods have been alsosuggested in the literature to tackle the problem in its utmost generality, e.g. [7]. The current state-of-the-art samplingmethod for problems of this type is termed Subset Simulation (SuS) [8] and belongs to the family of MCMC techniques.Within the context of Subset Simulation, various random-walk and non-random-walk-based MCMC proposal steps[9, 10] have been explored and suggested, to improve the sampling efficiency of SuS, including Hamiltonian steps [11]. a r X i v : . [ s t a t . C O ] S e p n this work we completely deviate from SuS and we introduce a gradient-based Hamiltonian Markov Chain MonteCarlo (HMCMC) sampling framework, termed Approximate Sampling Target with Post-processing Adjustment (ASTPA)[12], that is directly used for rare events probabilities estimation. The basic idea of ASTPA is to construct a relevanttarget distribution to sample from, by weighting the high-dimensional random variable space through a one-dimensionallikelihood model, using the limit-state function, and to then utilize an original post-sampling step, using an inverseimportance sampling procedure based on the acquired samples. Hamiltonian MCMC schemes are employed toperform the sampling. The Hamiltonian Monte Carlo (HMC) method, originally developed by [13], and more recentlypopularized mainly through the works of [14, 15, 16, 17], is characterized by scalability [15, 18, 17], fast mixing rates,weak sample auto-correlation, even in complex high-dimensional parameter spaces [19, 20, 21], and has achievedbroad-spectrum successes in most general settings e.g. [22, 23, 24, 25]. Herein, we compare the performance of thetypical HMCMC scheme with our newly developed Quasi-Newton based mass preconditioned HMCMC algorithmthat also exploits the information about the localized geometry of the failure region, through an inexpensive BFGSapproximation. The involved user-defined parameters of ASTPA are also discussed in the paper and general defaultvalues are suggested. The performance of the proposed methodology is finally examined and compared successfullyagainst Subset Simulation, in a series of static and dynamic, low- and high-dimensional benchmark problems.
In HMCMC methods, Hamiltonian dynamics are used to produce distant state steps for the Metropolis proposals,thereby avoiding the slow exploration of the state space that results from the diffusive behavior of simple random-walkproposals. Given a parameter of interest θ with (unnormalized) density π Θ (.), the Hamiltonian Markov Chain MonteCarlo method introduces an auxiliary momentum variable z and samples from the joint distribution characterized by: π ( θ , z ) ∝ π Θ ( θ ) π Z | Θ ( z | θ ) (1)where π Z | Θ ( . | θ ) is proposed to be a symmetric distribution. With π Θ ( θ ) and π Z | Θ ( z | θ ) being uniquely described upto normalizing constants, the functions U ( θ ) = − log π Θ ( θ ) and K ( θ , z ) = − log π Z | Θ ( z | θ ) are introduced as thepotential energy and kinetic energy, owing to the physical laws which motivate the Hamiltonian Markov Chain MonteCarlo algorithm. The total energy H ( θ , z ) can be thus expressed as: H ( θ , z ) = U ( θ ) + K ( θ , z ) (2)and is often termed the Hamiltonian H . The kinetic energy function is unconstrained and can be formed in various waysbased on the implementation. In most typical cases, the momentum is given by a zero-mean normal distribution [15, 19],and accordingly the kinetic energy can be written as: K ( θ , z ) = − log π Z | Θ ( z | θ ) = − log π Z ( z ) = z T M − z , wherethe M is a symmetric, positive-definite covariance (mass) matrix.HMCMC generates a Metropolis proposal on the joint state-space ( θ , z ) by sampling the momentum and simulatingtrajectories of Hamiltonian dynamics in which the time evolution of the state ( θ , z ) is governed by Hamilton’s equations,expressed typically by: d θ dt = ∂H∂ z = ∂K∂ z = M − z , d z dt = − ∂H∂ θ = − ∂U∂ θ = ∇ θ L ( θ ) (3)where L ( θ ) denotes the log-density of the target distribution. Hamiltonian dynamics prove to be an effective proposalgeneration mechanism because the distribution π ( θ , z ) is invariant under the dynamics of Eq. (3). These dynamicsenable a proposal state, obtained by an approximate solution of Eq. (3), to be distant from the current state, yet havinghigh probability of acceptance. The solution to Eq. (3) is in general analytically intractable and thus the Hamiltonianequations need to be numerically solved by discretizing time, using some small step size, ε . A symplectic integratorthat can be used for the numerical solution is the leapfrog one, as follows: z t+ ε /2 = z t − ( ε ∂U∂ θ ( θ t ) , θ t+ ε = θ t + ε ∂K∂ z ( z t+ ε /2 ) , z t+ ε = z t+ ε /2 − ( ε ∂U∂ θ ( θ t+ ε ) (4)The main advantages of using the leapfrog integrator are its simplicity, its volume-preserving feature, and its reversibility,due to its symmetry, by simply negating z , facilitating a valid Metropolis proposal. See [15], [19] and[26] for details onenergy-conservation, reversibility and volume-preserving integrators and their connections to HMCMC. It is noted thatin the above leapfrog integration algorithm, the computationally expensive part is to acquire the ∂U∂ θ term at the updatedlocation θ . Taking L = τ /ε steps of the leapfrog integrator approximates the evolution ( θ (0) , z (0)) −→ ( θ ( τ ) , z ( τ )) ,where τ is the trajectory length or path length, and provides the exact solution in the limit ε −→ .As discussed, the typical HMCMC version is based on a Gaussian momentum π Z | Θ ( z | θ ) = π Z ( z ) ∼ N ( , M ) (or z ∼ N ( , M ) ). The mass matrix M is often set to the identity matrix, I , but can also be adapted to precondition the2 lgorithm 1 Hamiltonian Markov Chain Monte Carlo procedure HMCMC( θ , ε , L , L ( θ ) , NIter ) for m = 1 to N Iter do z ∼ N (0 , I ) (cid:46) momentum sampling from standard normal distribution θ m ← θ m − , ˜ θ ← θ m − , ˜ z ← z for i = 1 to L do ˜ θ , ˜ z ← Leapfrog( ˜ θ , ˜ z , ε ) (cid:46) leapfrog integration end for with probability : α = min (cid:26) exp( L ( ˜ θ ) − z . ˜ z )exp( L ( θ m − ) − z . z ) (cid:27) (cid:46) Metropolis step θ m ← ˜ θ , z m ← - ˜ z end for end procedure sampler when relevant information about the target distribution is available (see Section 4). A standard procedure fordrawing NIter samples via HMCMC is described in Algorithm 1, where L ( θ ) is the log-density of the target distributionof interest. θ are the initial values for the θ , and L is the number of leapfrog steps, as explained before. For eachHMCMC step, we first resample the momentum and then implement the L leapfrog updates (Leapfrog( ˜ θ , ˜ z , ε )) beforewe accept or reject the Metropolis proposal at the pertinent step.The efficiency of HMCMC relies significantly on selecting suitable values for ε and L . In this work we select thestepsizes ε in such a way that the corresponding average acceptance rates are approximately 65 % , as values between60 % and 80 % are typically assumed optimal [15, 16, 18]. The dual averaging algorithm of [16] was adopted here tofind these stepsizes. To determine the value of L , we estimate the trajectory length τ so as to have a sufficient so callednormalized Expected Square Jumping Distance ( ESJD ) τ − / E (cid:107) θ ( t +1) ( τ ) − θ ( t ) ( τ ) (cid:107) , as introduced in [27], and thenwe randomly perturb each trajectory length τ ( t ) in the range [0 . τ, . τ ] to avoid periodicity ( t denotes the t -th iterationof HMCMC). In all our experiments we determine L and control the trajectory length in this manner, as we have foundit to work well in practice. The role of these parameters ( ε and τ (or L )) and techniques for determining them have beenquite extensively studied and for more details we refer the readers to [15, 16, 18]. The failure probability P F for a system, that is the probability of a defined unacceptable system performance, can beexpressed as a d -fold integral, as: P F = E [ I F ( θ )] = (cid:90) g ( θ ) ≤ I F ( θ ) π θ ( θ ) d θ (5)where θ is the random vector [ θ , ..., θ d ] T ; F ⊂ R d is the failure event in the parameter space; g( θ ) is the limit-statefunction that can include one or several distinct failure modes and defines the failure of the system by g( θ ) ≤ I(.) is the indicator function with: I F ( θ ) = 1 if θ ∈ g( θ ) ≤ I F ( θ ) = 0 otherwise; E is the expectation operator, and π θ is the joint probability density function (PDF) for Θ . It is common practice in reliability analysis to have the jointPDF of Θ be 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 Θ , e.g. [28]. When thisis not the case however, but the probabilistic characterization of X can be defined in terms of marginal distributionsand correlations, the Nataf distribution (equivalent to Gaussian copula) can be used to model the joint PDF, and themapping to the standard normal space can be then accomplished [29].The main idea of our approach to calculate the failure probability is to construct an appropriate approximate targetdistribution to sample from, based on Hamiltonian MCMC methods that can quickly reach regions of interest and cankeep the number of model calls to a minimum, and to then utilize a post-sampling step to acquire the exact probabilityestimation, without any additional model calls. We construct this approximate target distribution by combining themultidimensional parameter space Θ with a one-dimensional likelihood function, using the limit-state expression.This one-dimensional likelihood function is expressed as a Gaussian PDF with mean = µ g ( θ ) = 0 , where g ( θ ) is the3 igure 1: The above figures represent the analytical target distribution, the simulated target distribution samples based on our HMCMC-based method, and the fitted Gaussian Mixture Model describing the simulated samples, from left to right, respectively. limit-state function, and a dispersion factor σ : N (cid:18) g ( θ ) g c µ g ( θ ) = 0 , σ (cid:19) , g c = (cid:26) g ( ) , if g ( ) > or g ( ) < , otherwise (6)where g c is a normalizing constant. The reason for this normalization, g ( θ ) /g c , is to control the suggested upper andlower bounds of σ . The target PDF is then defined as:Target probability distribution ∝ N (cid:18) g ( θ ) g c µ g ( θ ) = 0 , σ (cid:19) × (cid:18) θ ∼ N ( , I ) (cid:19) (7)Having the total number of model calls in mind, as well as the coefficient of variation of the estimator (C.O.V), thesuggested value for σ is in the range [0 . . . Fine tuning σ in that range is not generally necessary. It is recommended,in general, to use higher σ values (0 . − . in nonlinear high-dimensional problems and multi-modal cases when ≤ g ( ) ≤ , since a larger σ usually allows longer state jumps and fewer required model calls. On the other hand, alower σ generally increases the accuracy of the estimator, at the expense of a slightly increased number of model calls.Fig. 1 concisely portrays the overall approach by using a bimodal target distribution. The gray curves represent theparabolic limit-state function g ( θ ) of this problem, with the failure domain being outside g ( θ ) . The left figure displaysthe constructed target distribution, by adopting the previously described approach, which in this simple 2D case canbe visualized. The middle figure shows drawn samples from the target distribution by our suggested HamiltonianMCMC variant, described in Section 4. For their initial stage, our HMCMC samplers have an adaptive annealed phase,mainly in order to automatically tune parameters and reduce the computational cost, overall, and then follow the typicalHamiltonian approach, except in our Quasi-Newton case (Section 4) the mass matrix is appropriately preconditioned.As such, during the burn-in period, we initialize the spread of likelihood, σ , equal to that then follows an exponentialdecay throughout the burn-in period, while at the end of this initial period, σ takes its constant value, as described above,for the stationary phase of the algorithms.To finally compute the failure probability we have to adjust Eq. (5) accordingly, since the samples have been sampledbased on our constructed approximate target distribution. An original post-sampling step is devised at this stage usingour inverse importance sampling procedure, i.e. having the samples, choose a pertinent Importance Sampling Density(ISD) automatically, based on the samples. Given that, the probability of failure after some algebra (see [12] for details)can be computed as follows: P F = (cid:90) I F ( θ ) π θ ( θ ) d θ = (cid:90) I F ( θ ) C . ˜ h ( θ ) (cid:96) ( θ ) d θ (8)where C = N (cid:80) Ni =1 ˜ h ( θ i ) Q ( θ i ) ; ˜ h ( . ) denotes the non-normalized target PDF, (cid:96) ( θ ) is our likelihood function, and Q ( . ) isa computed Gaussian Mixture Model (GMM), based on the already available samples and the generic ExpectationMaximization (EM) algorithm, as indicatively seen in the right plot of Fig. 1.Our described newly proposed method is termed ASTPA (Approximate Sampling Target with Post-processing Adjust-ment) and, as a summary, comprises of constructing a target distribution model, performing HMCMC sampling, andfinally applying a post-sampling step. For more details on supplementary justifications about this method, we referreaders to [12]. 4 QUASI-NEWTON EXTENSIONS AND CONNECTIONS TO HMCMC
In high-dimensional problems, the computational cost of the typical HMCMC sampler may increase considerably and aprohibitive number of model calls per leapfrog step may be required. In this work, we address this issue in a developedNewton-type context, where the Hessian information is approximated without any required additional model calls perleapfrog step. To this end, the well-known BFGS approximation [30] is used in our Quasi-Newton type HamiltonianMCMC approach. Let θ ∈ R d , consistent with the previous section. Given the k - th estimate W k , where W k is anapproximation to the inverse Hessian at θ k , the BFGS update W k +1 can be expressed as: W k +1 = ( I − s k y Tk y Tk s k ) W k ( I − y k s Tk s Tk y k ) + s k s Tk s Tk y k (9)where I is the identity matrix, s k = θ k +1 − θ k , and y k = ∇ f ( θ k +1 ) − ∇ f ( θ k ) where f : R d −→ R denotes anyrelevant target distribution function in this case. Our developed Quasi-Newton preconditioned Hamiltonian MarkovChain Monte Carlo (QNp-HMCMC) method is presented in detail in Algorithm 2. In the burn-in phase we are stillsampling the momentum from an identity mass matrix but the ODEs of Eq. (3) now become: ˙ θ = WM − z , ˙ z = W ∇ θ L ( θ ) . (10)where W ∈ R d × d is the symmetric positive definite matrix of Eq. (9) and being the inverse Hessian matrix provides aninformed approximation of the local geometry of the parameter space, accelerating exploration of the domain. Thefinal estimation of the approximated inverse of the Hessian matrix, W , from the burn-in phase is then used to definethe preconditioned covariance matrix to sample the momentum variable for the stationary, non-adaptive stage of thechain. It can be shown that all utilized dynamics in both phases of the algorithm enable us to maintain the desired targetdistribution as the invariant one. In Section 5 we empirically evaluate and compare the QNp-HMCMC performance invarious settings. For further details on the QNp-HMCMC method, its performance in different settings, and its validity,see [12]. In this section, four numerical examples are implemented to illustrate the efficiency of the proposed methods. Inall examples, the tuning parameters ( ε , τ , σ ) are systematically used as mentioned in Secs. 2 and 3. In the contextof reliability problems, we use the default value τ = 0 . as a starting point and then employ the ESJD metric [27]as described in Section 2. The burn-in period is chosen to be on average 15% of the total number of model calls,while the upper bound of the burn-in size is limited to 20%. The described methods are compared to the Component-wise Metropolis-Hastings based Subset Simulation (CWMH-SuS). For the sake of comparison, we use two proposaldistributions in CWMH-SuS, a uniform distribution of width 2 and a standard normal one. The parameters of SubsetSimulation are chosen as n s = n s isthe number of samples for each subset level, and p = 0 . , where p is the percentile of the samples that determinesthe intermediate subsets [8]. Comparisons are illustrated in terms of accuracy and computational cost. In particular,the tables show the P F estimation, including the mean number of limit-state function calls in order to calculate thevalue and gradient of the target distribution, the analytical gradients are provided in all examples, in the HMCMC-based algorithms, and the value of the limit-state function in SuS. In all examples, the number of limit-state functionevaluations for all methods has been set to be roughly the same to each other for comparison purposes. Results are basedfor all examples on 500 independently performed simulations, so that the sample mean and C.O.V of the results can beacquired. It should be noted that the ASTPA parameters are carefully chosen for all examples but are not optimized forany one. Hence, comparative and perhaps improved alternate performance might be achieved with a different set ofparameters. The first example is expressed by the following limit state function for two standard normal random variables [31]: g ( θ ) = r − θ − κ ( θ − e ) (11)where r , κ and e are deterministic parameters chosen as r = 6 , κ = 0 . and e = 0 . . The probability of failure is 3.95E-5 and the limit-state function consists of two design points (failure modes), as seen in Fig. 1. For the HMCMC-basedalgorithms, the likelihood dispersion factor, σ , is 0.7 and the burn-in sample size is taken as 200. Consistent to thediscussion in Section 3, the trajectory length is set to τ = 1 . Table 1 compares the number of model calls, the coefficientof variation and the E [ ˆ P F ] obtained by all tested methods. The Subset Simulation results are based on n s = E [ ˆ P F ] . Fig. 1 also demonstrates that the QNp-HMCMC samples accurately describe the two important failure regions.5 lgorithm 2 Quasi-Newton preconditioned Hamiltonian Markov Chain Monte Carlo procedure QN P -HMCMC( θ , ε , L , L ( θ ) , BurnIn , NIter ) W = I for m = 1 to N Iter do if m ≤ BurnIn then z ∼ N ( , M ) (cid:46) where M = I θ m ← θ m − , ˜ θ ← θ m − , ˜ z ← z , B ← W for i = 1 to L do ˜ θ , ˜ z ← Leapfrog-BurnIn( ˜ θ , ˜ z , ε , B ) Update W using Eq. (9) end for with probability : α = min (cid:26) exp( L ( ˜ θ ) − z . ˜ z )exp( L ( θ m − ) − z . z ) (cid:27) θ m ← ˜ θ , z m ← - ˜ z (cid:46) If proposal rejected: W ← B else (cid:46) If m > BurnIn z ∼ N ( , M ) (cid:46) where M = W − θ m ← θ m − , ˜ θ ← θ m − , ˜ z ← z for i = 1 to L do ˜ θ , ˜ z ← Leapfrog( ˜ θ , ˜ z , ε , M ) end for with probability : α = min (cid:26) exp( L ( ˜ θ ) − z . M − . ˜ z )exp( L ( θ m − ) − z . M − . z ) (cid:27) θ m ← ˜ θ , z m ← - ˜ z end if end for end procedure function L EAPFROG -B URN I N ( ˜ θ , ˜ z , ε, B ) ˜ z ← z + ( ε/ B ∇ θ L ( θ ) ˜ θ ← θ + ε B ˜ z ˜ z ← z + ( ε/ B ∇ θ L ( ˜ θ ) return ˜ θ , ˜ z . end function function L EAPFROG ( ˜ θ , ˜ z , ε, M ) ˜ z ← z + ( ε/ ∇ θ L ( θ ) ˜ θ ← θ + ε M − ˜ z ˜ z ← z + ( ε/ ∇ θ L ( ˜ θ ) return ˜ θ , ˜ z . end function able 1: Performance of various methods for the parabolic/concave limit-state function σ = 0 . τ = 1
500 Independent Simulations CWMH-SuS HMCMC QNp-HMCMC U ( − , N (0 , Number of model calls 4,559 4,565 4,391 4,926C.O.V 0.62 0.65 0.35 0.39 E [ ˆ P F ] (Exact P F ∼ Table 2: Performance of various methods for the four-branch series system σ = 0 . τ = 1 n s =
500 Independent Simulations CWMH-SuS HMCMC QNp-HMCMC U ( − , N (0 , Number of model calls 2,841 2,852 2,867 2,887C.O.V 0.26 0.30 0.29 0.26 E [ ˆ P F ] (Exact P F ∼ σ = 0 . τ = 1 n s = E [ ˆ P F ] (Exact P F ∼ This example is a well-known benchmark system reliability problem, defined by the following limit-state function inthe standard normal space: g ( θ ) = min . θ − θ ) − ( θ − θ ) / √
23 + 0 . θ − θ ) + ( θ − θ ) / √ / √
2) + ( θ − θ )(7 / √
2) + ( θ − θ ) (12)The trajectory length is chosen as τ = 1 and the likelihood dispersion factor, σ , is fixed to 0.7. The burn-in is set to 200samples. Table 2 shows that the SuS with uniform proposal gives more accurate P F estimation with smaller C.O.Vthan the HMCMC-based methods for the case of n s = n s = In this example, a nonlinear undamped single-degree-of-freedom (SDOF) oscillator subjected to a rectangular impulseload is analysed, as described in [32, 33]. The limit-state function is given as: g ( k , k , M, r, T , F ) = 3 r − (cid:12)(cid:12)(cid:12)(cid:12) F M ω sin( ω T (cid:12)(cid:12)(cid:12)(cid:12) (13)where ω = (cid:112) ( k + k ) /M is the natural frequency of the oscillator, T is the duration of the impulse load, M isthe mass, k and k are the stiffnesses of the primary and secondary springs, r is the displacement at which one ofthe springs yields, and F is the amplitude of the force. The description of all random variables is listed in Table 3.SuS results for both proposals are based on n s = µ F , of F . For the HMCMC-basedmethods, the trajectory length and the likelihood dispersion factor are chosen as τ = 0 . and σ = 0 . respectively. Theburn-in sample size is set to 500. It is shown in this example that the QNp-HMCMC approach provides significantlymore accurate and stable results in terms of the C.O.V. and E [ ˆ P F ] . Particularly for the lowest failure probability level,QNp-HMCMC approach noticeably outperforms all other methods. As results indicate, the QNp-HMCMC methodis roughly insensitive to the failure probability level and there is no negative influence on the method when changing µ F . For the two SuS variants, it is noteworthy to say here that the SuS with the standard normal proposal distributionindicates reasonably better performance in this example than the one with the uniform proposal.7 a) (b) Figure 2: (a) Simulated samples from the target distribution, (b) Analytical target distribution.Table 3: Random variables of the undamped oscillator
Variable Distribution Mean C.O.V M Gaussian 1 0.05 k Gaussian 1 0.1 k Gaussian 0.1 0.1 r Gaussian 0.5 0.1 T Gaussian 1 0.2 F Gaussian 0.6-0.45 Table 4: Performance of various methods for the undamped oscillator example σ = 0 . τ = 0 . µ F = 0 .
500 Independent Simulations CWMH-SuS HMCMC QNp-HMCMC U ( − , N (0 , Number of model calls 5,170 5,160 5,132 5,119C.O.V 0.67 0.51 0.14 0.11 E [ ˆ P F ] (Exact P F ∼ σ = 0 . τ = 0 . µ F = 0 . Number of model calls 7,583 7,617 7,523 7,515C.O.V 0.77 0.70 0.21 0.15 E [ ˆ P F ] (Exact P F ∼ In this last example, we consider a SDOF oscillator, initially at rest, with natural frequency ω = 7 . rad/s anddamping ratio ξ = 0 . , subjected to a Gaussian white noise ( W ( t ) ) excitation with spectral density of magnitude S = 1 . The response of the system is computed at discrete time instants { t j = ( j − t : j = 1 , ...n } with ∆ t = 0 . , and the duration of study is T = 5 sec . Thus, the number of time instants is equal to n = T / ∆ t + 1 = 101 .The state vector θ is the sequence of i.i.d. standard normal random variables that generate the W ( t j ) = (cid:113) πS ∆ t θ j atthe discrete time instants, resulting in involved random variables in this example. Failure is characterized by thepositive displacement response exceeding a threshold level R : g ( θ ) = R − max { Y ( t ) } .The burn-in sample size is taken as 1,000 for the HMCMC-based methods. SuS results are based on n s = E [ ˆ P F ] . Compared to the HMCMC approach, this example agrees with the additional results in [12] and confirms thatthe application of QNp-HMCMC in high-dimensional reliability problems is in general more attractive. By decreasing8 able 5: Performance of various methods for SDOF oscillator under white noise σ = 0 . τ = 0 . R = 1 .
500 Independent Simulations CWMH-SuS HMCMC QNp-HMCMC U ( − , N (0 , Number of model calls 11,000 11,011 11,063 11,059C.O.V 0.32 0.35 0.30 0.24 E [ ˆ P F ] (Exact P F ∼ σ = 0 . τ = 0 . R = 2 Number of model calls 13,578 13,646 13,644 13,618C.O.V 0.41 0.48 0.34 0.29 E [ ˆ P F ] (Exact P F ∼ the target failure probability, results also reveal that QNp-HMCMC gives us a substantially improved estimation incomparison to all other methods. A novel approach for estimation of rare event probabilities termed Approximate Sampling Target with Post-processingAdjustment (ASTPA), is presented in this paper, suitable for low- and high-dimensional problems, very small proba-bilities and multiple failure modes. ASTPA can provide an accurate unbiased estimation of the failure probabilitieswith an efficient number of limit-state function evaluations. The basic idea of ASTPA is to construct a relevant targetdistribution by weighting the high-dimensional random variable space through a one-dimensional likelihood model,using the limit-state function. To sample from this target distribution we utilize gradient-based HMCMC schemes,including our newly developed Quasi-Newton based mass preconditioned HMCMC algorithm (QNp-HMCMC) that cansample very adeptly, particularly in difficult cases with high-dimensionality and very small failure probabilities. Finally,an original post-sampling step is also devised, using an inverse importance sampling procedure based on the samples.The performance of the proposed methodology is examined and compared very successfully herein against SubsetSimulation in a series of static and dynamic low- and high-dimensional benchmark problems. As a general guideline,QNp-HMCMC is recommended to be used for problems with more than 20 dimensions, where traditional HMCMCschemes may not perform that well. However, even in lower dimensions QNp-HMCMC performs reasonably well andis still a competitive algorithm. Since we are utilizing gradient-based sampling methods, all of our analyses and resultsare based on the fact that analytical gradients can be computed. In cases where numerical schemes are needed for thegradient evaluations, then HMCMC methods will not be competitive in relation to SuS. It should also be pointed outthat different combinations of the HMCMC and QNp-HMCMC algorithms can be possible, based on problem-specificcharacteristics. Some of the ongoing and future work is directed towards exploring various ASTPA variants, and onestimating first-passage problems under numerous settings and high-dimensional parameter spaces.
REFERENCES [1] Ali Bakhshi and Hamed Nikbakht. Loading pattern and spatial distribution of dynamic wind load and comparisonof wind and earthquake effects along the height of tall buildings. In
Proceedings of the 8th International Conferenceof Structural Dynamics, EURODYN , pages 1607–1614, 2011.[2] Hamed Esmaeili, Ali Kheyroddin, Mohammad Ali Kafi, and Hamed Nikbakht. Comparison of nonlinear behaviorof steel moment frames accompanied with rc shear walls or steel bracings.
The Structural Design of Tall andSpecial Buildings , 22(14):1062–1074, 2013.[3] Rüdiger Rackwitz. Reliability analysis—A review and some perspectives.
Structural Safety , (4):365–395, 2001.[4] Armen Der Kiureghian. First-and second-order reliability methods. Engineering Design Reliability Handbook , , 2005.[5] Pei-Ling Liu and Armen Der Kiureghian. Optimization algorithms for structural reliability. Structural Safety , (3):161–177, 1991.[6] Karl Breitung. 40 years FORM: Some new aspects? Probabilistic Engineering Mechanics , :71–77, 2015.[7] Gerhart I Schuëller and Helmuth J Pradlwarter. Benchmark study on reliability estimation in higher dimensions ofstructural systems–an overview. Structural Safety , (3):167–182, 2007.[8] Siu-Kui Au and James L Beck. Estimation of small failure probabilities in high dimensions by Subset Simulation. Probabilistic Engineering Mechanics , (4):263–277, 2001.99] Iason Papaioannou, Wolfgang Betz, Kilian Zwirglmaier, and Daniel Straub. MCMC algorithms for SubsetSimulation. Probabilistic Engineering Mechanics , :89–103, 2015.[10] Konstantin M Zuev. Subset Simulation method for rare event estimation: An introduction. Encyclopedia ofEarthquake Engineering , pages 1–25, 2015.[11] Ziqi Wang, Marco Broccardo, and Junho Song. Hamiltonian Monte Carlo methods for Subset Simulation inreliability analysis.
Structural Safety , :51–67, 2019.[12] Hamed Nikbakht and Konstantinos G. Papakonstantinou. Hamiltonian MCMC methods for estimating rare eventsprobabilities in high-dimensional problems. Journal of Computational Physics , under review.[13] Simon Duane, Anthony D Kennedy, Brian J Pendleton, and Duncan Roweth. Hybrid Monte Carlo.
Physics LettersB , (2):216–222, 1987.[14] Radford M Neal. Bayesian learning for neural networks. PhD Thesis, University of Toronto , 1995.[15] Radford M Neal. MCMC using Hamiltonian dynamics.
Handbook of Markov Chain Monte Carlo , (11):2, 2011.[16] Matthew D Hoffman and Andrew Gelman. The No-U-Turn Sampler: Adaptively setting path lengths in Hamilto-nian Monte Carlo. Journal of Machine Learning Research , (1):1593–1623, 2014.[17] Mark Girolami and Ben Calderhead. Riemann manifold langevin and Hamiltonian Monte Carlo methods. Journalof the Royal Statistical Society: Series B (Statistical Methodology) , 73(2):123–214, 2011.[18] Alexandros Beskos, Natesh Pillai, Gareth Roberts, Jesus-Maria Sanz-Serna, and Andrew Stuart. Optimal tuningof the Hybrid Monte Carlo algorithm.
Bernoulli , (5A):1501–1534, 2013.[19] Michael Betancourt. A conceptual introduction to Hamiltonian Monte Carlo. arXiv preprint arXiv:1701.02434 ,2017.[20] Mohammad Mahdi Kamani, Farshid Farhat, Stephen Wistar, and James Z Wang. Shape matching using skeletoncontext for automated bow echo detection. In , pages901–908. IEEE, 2016.[21] Mohammad Mahdi Kamani, Farshid Farhat, Stephen Wistar, and James Z Wang. Skeleton matching withapplications in severe weather detection. Applied Soft Computing , 70:1154–1166, 2018.[22] Andrew Gelman, Hal S Stern, John B Carlin, David B Dunson, Aki Vehtari, and Donald B Rubin.
Bayesian DataAnalysis . Chapman and Hall/CRC, 2013.[23] John Kruschke.
Doing Bayesian Data Analysis: A tutorial with R, JAGS, and Stan . Academic Press, 2014.[24] Cole C Monnahan, James T Thorson, and Trevor A Branch. Faster estimation of Bayesian models in ecologyusing Hamiltonian Monte Carlo.
Methods in Ecology and Evolution , (3):339–348, 2017.[25] Elena Akhmatskaya and Sebastian Reich. GSHMC: An efficient method for molecular simulation. Journal ofComputational Physics , (10):4934–4954, 2008.[26] Nilesh Tripuraneni, Mark Rowland, Zoubin Ghahramani, and Richard Turner. Magnetic Hamiltonian Monte Carlo. arXiv preprint arXiv:1607.02738 , 2016.[27] Ziyu Wang, Shakir Mohamed, and Nando Freitas. Adaptive Hamiltonian and Riemann manifold Monte Carlo. In International Conference on Machine Learning , pages 1462–1470, 2013.[28] Michael Hohenbichler and Rudiger Rackwitz. Non-normal dependent vectors in structural safety.
Journal of theEngineering Mechanics Division , (6):1227–1238, 1981.[29] Armen Der Kiureghian and Pei-Ling Liu. Structural reliability under incomplete probability information. Journalof Engineering Mechanics , (1):85–104, 1986.[30] Jorge Nocedal and Stephen Wright. Numerical Optimization . Springer, 2006.[31] Armen Der Kiureghian and Taleen Dakessian. Multiple design points in first and second-order reliability.
StructuralSafety , (1):37–49, 1998.[32] Christian G Bucher and Ulrich Bourgund. A fast and efficient response surface approach for structural reliabilityproblems. Structural Safety , (1):57–66, 1990.[33] Roland Schöbi and Bruno Sudret. Structural reliability analysis for p-boxes using multi-level meta-models. Probabilistic Engineering Mechanics ,48