Rare event estimation using stochastic spectral embedding

Estimating the probability of rare failure events is an essential step in the reliability assessment of engineering systems. Computing this failure probability for complex non-linear systems is challenging, and has recently spurred the development of active-learning reliability methods. These methods approximate the limit-state function (LSF) using surrogate models trained with a sequentially enriched set of model evaluations. A recently proposed method called stochastic spectral embedding (SSE) aims to improve the local approximation accuracy of global, spectral surrogate modelling techniques by sequentially embedding local residual expansions in subdomains of the input space. In this work we apply SSE to the LSF, giving rise to a stochastic spectral embedding-based reliability (SSER) method. The resulting partition of the input space decomposes the failure probability into a set of easy-to-compute \rev{conditional} failure probabilities. We propose a set of modifications that tailor the algorithm to efficiently solve rare event estimation problems. These modifications include specialized refinement domain selection, partitioning and enrichment strategies. We showcase the algorithm performance on four benchmark problems of various dimensionality and complexity in the LSF.


Introduction
Ensuring the reliability of structures and systems is a core task in many engineering disciplines.This includes the reliability of machines (Bertsche, 2008), medical devices (Fries, where the indicator function 1 D f (x) takes the value 1 in the failure domain D f and 0 everywhere else: 1, if g(x) ≤ 0, 0, if g(x) > 0. (2) Analytical computation of the failure probability is rarely possible in practice and direct numerical integration (e.g., via quadrature) is often hindered by the inherently small scale of the failure probability and the potentially high input dimensionality M .For these reasons, initial efforts in the reliability literature focused on developing methods that approximate the limit-state surface (Basler, 1960;Hasofer and Lind, 1974;Rackwitz and Fiessler, 1978;Zhang and Der Kiureghian, 1995;Hohenbichler et al., 1987;Breitung, 1994;Tvedt, 1990;Cai and Elishakoff, 1994).These approximation methods remain competitive for a certain class of reliability problems even today, but there exist well known examples where the shape of the limit-state surface leads to gross errors in the estimated failure probability (e.g., high-dimensional, non-linear problems Valdebenito et al. (2010)).In those cases the direct estimation of the failure probability by means of probabilistic simulation methods is the method of choice.The most widely used methods for this are Monte Carlo simulation (MCS, Fishman (2011); Rubinstein and Kroese (2016)), approximation methods in conjunction with importance sampling (IS, Hohenbichler and Rackwitz (1988); Melchers (1999)), its adaptive variants (sequential IS, Papaioannou et al. (2016) and cross entropy-based IS, Kurtz and Song (2013); Papaioannou et al. (2019)) and subset simulation (SuS, Au and Beck (2001); Zuev et al. (2012)).Other noteworthy methods that can be used to directly estimate the failure probability are radial-based importance sampling (RBIS, Harbitz (1986)), directional simulation (DS, Bjerager (1988)) and line sampling (LS, Koutsourelakis et al. (2004); Papaioannou and Straub (2021)).
Computationally expensive models remain challenging to analyse, even with the most efficient simulation methods.In this setting, cheap to evaluate surrogate models can lead to significant computational savings, especially with so-called active learning reliability methods.These methods approximate the limit-state function g using a sequence of increasingly accurate surrogate models.The first method in this family was efficient global reliability analysis (EGRA, Bichon et al. (2008)).Arguably the best known method today is adaptive Kriging MCS (AK-MCS, Echard et al. (2011)) and its variants (Echard et al., 2013;Huang et al., 2016).A recent review of the field is given in Moustapha et al. (2021).
In this paper we propose a novel active learning reliability method that utilizes the recently proposed stochastic spectral embedding method (SSE, Marelli et al. (2021)) and more specifically the active learning sequential partitioning approach developed for Bayesian inverse problems in Wagner et al. (2021).The proposed approach, termed stochastic spectral embedding-based reliability (SSER), benefits from a handful of modifications to the original SSE, namely new refinement domain selection, partitioning and sample enrichment schemes, which we present in Section 2. SSER constructs a set of spectral representations of g with higher accuracy in proximity of the limit-state surface.The proposed adaptive partitioning strategy results in a set of subdomains in the vicinity of the limit-state surface.Provided that g is sufficiently well-behaved in these subdomains, our approach is capable of efficiently addressing problems of varying complexity (i.e., non-linear computational models, high dimensional problems, small failure probabilities).We showcase this property by comparing SSER to other state-of-the-art methods on four applications in Section 3.
2 Stochastic spectral embedding-based reliability

Stochastic spectral embedding
Stochastic spectral embedding (SSE, Marelli et al. (2021)) is a function approximation technique that works by sequentially expanding residuals R k S (X) in subdomains of the parameter space D k X ⊆ D X .At each level in the sequence, SSE constructs local spectral expansions of the residual of the current approximation at a set of subdomains, as illustrated in Figure 1.
As for any other metamodeling techniques, this approach can be used in reliability analyses to directly approximate the limit-state function: In this expression, K ⊆ N 2 is a set of index pairs with elements k = ( , p) used to identify the residual expansions in the SSE representation.The first element of the index pair ∈ {0, • • • , L} denotes the level at which the residual expansion resides, where L is the total number of levels.The second element p ∈ {1, • • • , P } denotes the index of the residual expansion at a given level, with P denoting the total number of expansions at the -th level.
1 D k X is the indicator function that is equal to 1 in the k-th subdomain and 0 everywhere else.For further derivations, it is beneficial to distinguish between standard and terminal domains among the card(K) domains used in the SSE representation.A terminal domain is any domain that has not been split.It can be located at any level of the SSE representation and the set of all terminal domains with indices gathered in T constitutes a complete partition of the input domain, i.e., D X = k∈T D k X (top domains highlighted in orange in Figure (a)).With this, Eq. ( 3) is expanded to In principle, any spectral expansion technique could be used to expand the residuals (Marelli et al., 2021).In the present contribution, we consider only polynomial chaos expansions (PCEs) such that where A k is a multi-index set describing the dimension-wise polynomial degrees, {Ψ k α } α∈A k are mutually orthogonal polynomial basis functions and a k α denote the corresponding coefficients (Xiu and Karniadakis, 2002;Blatman and Sudret, 2011;Lüthen et al., 2021).
For computational and bookkeeping ease, we manage the SSE construction in what we call the quantile space, i.e., the M -dimensional unit hypercube denoted by D U .As discussed in Marelli et al. (2021), for every original or real input space D X there exists an isoprobabilistic mapping to D U (Figure (a)), such that the partitioning of the input space can be conducted in the quantile space without any loss of generality.
The SSE representation shown in Figure 1 ultimately partitions the input space into card(T ) disjoint subdomains such that k∈T D k X = D X .Using this representation, the failure probability from Eq. ( 1) can be rewritten as a sum of conditional failure probabilities where V k is the domain-wise probability mass that is readily available from the construction in the quantile space, and P k f is the k-th conditional failure probability given by This conditional failure probability is estimated here with suitable probabilistic simulation methods utilizing the SSE representation of the limit-state function (more details are provided in Section 2.3).Assuming the simulation method is accurate enough, the accuracy of the failure probability estimate depends only on that of the SSE representation itself.To devise an adaptive algorithm, it is necessary to devise a suitable error measure on each conditional failure probability estimate.As shown in Marelli et al. (2021), the SSE representation provides domain-wise prediction error estimators, e.g., by means of crossvalidation (e.g., leave-one-out error).Such estimators can in principle be used as a proxy to assess the accuracy of P k f .However, these estimators are not robust in the presence of very localized behaviour that is typical for limit-state functions.An improved estimator that takes into account the point-wise prediction accuracy can be derived by means of the local error estimation of bootstrap PCE (bPCE), originally introduced in a reliability context in Marelli and Sudret (2018).

Bootstrap SSE
The idea of bPCE is to use bootstrap resampling (Efron, 1979) on the original experimental design to create a set of B experimental designs and construct B sparse PCEs with them.
The individual bootstrap PCE realizations can then be used to estimate point-wise statistics of the PCE prediction.This idea can be applied to the residual expansions within the SSE approach.Due to the sequential embedding of the residual expansions, the prediction error of the full SSE approximation can be assessed solely based on the prediction errors in the terminal domains (Marelli et al., 2021).In every terminal domain, we therefore construct B residual expansions to obtain B bootstrap predictions By analogy with Marelli and Sudret (2018), we can define an estimate of the point-wise misclassification probability as: where 1 D f,SSE and 1 D f,SSE are evaluated with the mean and bootstrap predictors of SSE respectively (see also Eq. ( 2)).
The bootstrap replications can also be used to directly obtain statistics of the conditional failure probability estimates, such as variance or confidence bounds.As an example, the conditional failure probability variance is given by with the failure probability of the b-th replication given by P k,(b) By means of Eq. ( 6) this conditional failure probability variance can be used to write an expression for the total failure probability variance as This expression assumes independence between the conditional failure probability estimates.
Due to the ancestry relation between subdomains, we cannot in general rule out dependence between those quantities, however, we assume that it is negligible because all expansions in terminal domains are constructed with independently enriched samples.
Similarly, the bootstrap replications can be used to define confidence bounds on the total failure probability.Let P (b) be the total failure probability estimated with the b-th replication from all terminal domains.The equal tail confidence bounds P k f , P k f are then defined by with α ∈ [0, 0.5] and γ def = 1 − 2α is called the symmetrical confidence level.It is common practice to take γ = 95% and α = 2.5%.

Dependent input parameters
Generally, the input random vector X may have mutually dependent marginals that we treat with copula theory in this work (Nelsen, 2006;Joe, 2015).A typical approach to address dependence in reliability problems is to use isoprobabilistic mappings to an independent space, e.g., the standard normal space.However, sparse spectral techniques acquire their strength from the sparsity of effects principle (Montgomery, 2019), which states that in many engineering models the majority of the output is attributable to low-interaction order terms.A mapping to the standard normal space might therefore reduce the efficiency of those techniques.
Nonetheless, dependence is challenging for the construction of the spectral basis functions Ψ k α in Eq. ( 5).However, as extensively studied in Torre et al. (2019a), polynomial chaos expansions seem to be most accurate for predictive purposes when dependence is simply ignored and the basis is derived by tensor product of univariate bases orthogonal to the marginals instead.Applying this to SSE suggests it is best to construct the domain-wise polynomial basis by assuming the domain-wise marginals to be independent.
As mentioned before, the partitioning into orthogonal subdomains is defined in the independent quantile space (Marelli et al., 2021).Therefore, the transformation of those domains to the dependent physical space results in non-orthogonal domains (see Figure 2).
Constructing the univariate polynomial bases requires knowledge of the marginal bounds in the physical space, which cannot be computed in a straightforward fashion.
To approximate these bounds, we therefore (1) randomly sample the (M −1)-dimensional boundary of the hypercube in the quantile space, (2) transform those points to the real space using a suitable isoprobabilistic transform (Rosenblatt, 1952;Torre et al., 2019b) and proceed to (3) compute the basis with an orthogonal enveloping hypercube in the physical space around those transformed points (see Figure (b)).It is interesting to note that for prediction purposes, the resulting multivariate basis will only be evaluated inside the subdomain boundary and not the subdomain boundary envelope, as the region outside the subdomain boundary belongs to neighbouring subdomains.

An updated sequential partitioning algorithm
Failure domains typically occupy only a small fraction of the input space.The original equal probability partitioning strategy presented in Marelli et al. (2021); Wagner et al. (2021) is therefore expected to converge only slowly to those domains.We therefore propose modifications that make the algorithm more efficient for reliability problems.The structure of the proposed algorithm follows the sequential partitioning algorithm used in the adaptive SSLE approach presented in Wagner et al. (2021).It is repeated here in an abbreviated form.

Refinement domain selection
At every step, the adaptive sequential partitioning algorithm chooses a refinement domain from the set of terminal domains T to be split and enriched with sample points.To choose a domain, the terminal domains are ranked based on their importance to the overall approximation.In the present setting, the goal of refinement is the reduction of uncertainty associated with the failure probability.Viewing Eq. ( 12) as a variance decomposition of the estimator variance, it is natural to refine the terminal domain that has the largest contribution to this variance estimator.More formally, out of the terminal domains we refine the domain with the largest This equation depends on the domain-wise probability mass V k and the variance of the conditional failure probability estimator Var P k f .In practice, Var P k f is estimated with Eq. ( 10), which in turn depends on the B estimates of . These estimates are computed with accurate simulation methods using the bootstrap SSE approximations (see Section 2.3).
If the SSE approximation in the current subdomain is not sufficiently accurate, it can happen that all bootstrap replications miss existing failure regions, ultimately resulting in a gross underestimation of E k .To avoid overlooking such domains as possible refinement domains, the algorithm additionally prioritizes them based on the probability mass V k .The refinement domain is ultimately chosen as 16) is met, which depends on the following intermediate re-prioritization criterion where P (i) is the total failure probability estimator at the i-th iteration of the algorithm and the threshold is heuristically chosen to ε Pf = 0.1%.This criterion is triggered when the failure probability estimator does not change significantly in three successive iterations.

Partitioning strategy
Once a refinement domain D k refine X is selected, the partitioning strategy determines how it is split.Refinement domains are split into two subdomains with unequal probability mass.
For notational simplicity, we omit the superscript k refine denoting the refinement domain in this section and emphasize instead that the following equations are valid for all k ∈ T .
For illustrative purposes, the partitioning strategy is exemplified for a simple problem in Figure 3.
To maximize enrichment efficiency, we base our splitting strategy on separating regions that predict failure/safety correctly from those that do so incorrectly.A measure of prediction accuracy in this respect is given by the misclassification probability P m (Echard et al., 2011) defined for bootstrap SSE in Eq. ( 9).To this end, we define two auxiliary conditional random vectors in the quantile space D U as These vectors identify regions of zero (Z 0 ) and non-zero (Z + ) misclassification probability in the current refinement domain.
Based on these random vectors, we compute the location of the split.We choose a set of splitting locations in each dimension υ i ∈ D Ui , i = 1, • • • , M that ideally completely separate the support of Z + from Z 0 (see Figure 3).Due to the likely disjoint support of each auxiliary random vector, the split should instead confine a maximum of Z + 's probability mass to one side of the split, while confining a maximum of Z 0 's probability mass to the other side of the split.Denoting by Z + i and Z 0 i the marginals of the auxiliary random vectors in the i-th dimension, we pick the splitting location where the objective function L i characterizes the split properties by returning the maximum of the respective auxiliary probability masses in the resulting subdomains.The splitting location υ i that maximizes L i thus splits the initial domain D X into two subdomains D (1)

Xi and D
(2) Xi that fulfil the initially stated goal of optimally separating regions that correctly predict failure/safety from regions that do so incorrectly.The objective function L i is bounded between [0, 1], where a value of 1 indicates a perfect split that creates one subdomain with only non-zero and one subdomain with all zero misclassification probability regions.
To ultimately choose a splitting direction d ∈ {1, • • • , M }, we compare the values of the objective functions L i ( υ i ) and split along the dimension that achieves the best split, i.e., In practice, a sufficiently large sample of the auxiliary random vectors Z 0 and Z + is used to conduct the computations of this section.This sample is readily available, as the SSE surrogate model can be evaluated at a negligible computational cost.

No misclassified sample points
It might occur that the the algorithm does not detect misclassified sample points at the time of partitioning.When this happens, the misclassification probability P m (Eq.( 9)) is zero everywhere and consequently the auxiliary random vector Z + is not defined.Because of the preceding refinement selection criterion in Eq. ( 15) this can only happen at the first step of the algorithm or after the intermediate re-prioritization criterion (Eq.( 16)) has been triggered.
To let the algorithm proceed in this case, the definition of the auxiliary random vectors is updated to discriminate based on the probability to predict values within a predefined empirical quantile around the limit-state surface.More formally, P m in Eq. ( 17) is then replaced with where 1 Dt is defined with an empirical quantile t such that P The value ε t = 0.01 has proven to be a good choice in practice.

Sample enrichment
After partitioning a selected refinement domain, in the original SSE algorithm points are added such that every subdomain contains exactly a total of N ref points.Because the partitioning strategy now allows the creation of domains with considerably different probability mass, the original sample enrichment approach is problematic.As an example, splitting in the tails of the input distribution will result in one large domain with N ref points from the previous step and one smaller domain with 0 points.In this case, no new points would be added to the large domain.The subsequently constructed residual expansion would be based on (almost) the same information as the parent domain expansion and be prone to over-fitting.We therefore employ a different strategy here, which always adds N ref points to every domain, independent of the existing points inside that domain.We refer to those points as the sample budget here.
The second change pertains to the placement of the sample budget.The original algorithm used random sampling to place points uniformly in the refinement subdomain.To exploit the fact that the SSE for reliability applications needs to be more accurate near the limit-state surface, we propose a slightly modified strategy: we split the sample budget and randomly place points uniformly in the refinement domain quantile space, where N curr is the number of existing points in the refinement domain.The remainder of the sample budget is used to randomly place sample points in the subset of the refinement domain with non-zero misclassification probability P m .This corresponds to sampling the auxiliary random vector Z + defined in Eq. ( 17).We do this to ensure stability of the residual expansions that are known to deteriorate in accuracy when the experimental design clusters in a subspace.
In the extreme case of no existing points (i.e., N curr = 0), this places half of the sample budget uniformly and the other half according to Z + .In the other extreme case of N curr ≥ N ref , the whole sample budget is devoted to sampling Z + .
At the first iteration, to construct the initial expansion, the enrichment occurs randomly in the entire input quantile space.We observed from repeated tests that the algorithm is less stable if the initial expansion is less accurate.

Stochastic spectral embedding-based reliability algorithm
Taking all presented ingredients, the full stochastic spectral embedding-based reliability (SSER) algorithm is written as follows: Input: • Input parameters X, limit-state function g S (X) of g(X) in the full domain D 0,1 X , retrieve its approximation error E 0,1 and initialize T = {(0, 1)} iii.Retrieve E +1,s from Eq. ( 14) (c) Retrieve Pf from Eq. ( 6) and its variance from Eq. ( 12) (d) If stopping criterion in Eq. ( 23) is met or less than two new expansions were created, terminate the algorithm, otherwise go back to 2 3. Termination

Stopping criterion
The algorithm can be terminated based on any stopping criterion from the active learning literature (Moustapha et al., 2021).In the present implementation, we opt for a criterion based on the stability of the reliability index bounds, also known as beta bounds.The generalized reliability index can be computed from the failure probability as β = −Φ −1 (P f ), where Φ −1 is the inverse normal CDF Ditlevsen (1979).Using the 95% confidence intervalls on β from the bootstrap replications as defined for the failure probability in Eq. ( 13), the stopping criterion is given by β where the failure threshold is heuristically set to ε β = 3%.For stability, we abort the algorithm only if this criterion is fulfilled in three consecutive iterations.
Additionally, the algorithm is terminated if the computational budget N ED , corresponding to the total number of admissible limit-state evaluations, has been exhausted.

Conditional failure probability estimation
In the SSE approach, the total failure probability is decomposed into a sum of weighted, conditional failure probabilities (see Eq. ( 6)).As a side effect of the proposed partitioning strategy (Section 2.2.2), these conditional failure probabilities typically end up being quite large or negligibly small.
Because of the SSE surrogate, simulation methods can be used at a low computational cost.The case of large conditional failure probabilities can be treated easily with simple Monte Carlo simulation.The case of small conditional failure probabilities, however, is more challenging even with surrogate models.Small failure probabilities should not be underestimated by the algorithm in large domains, because they can contribute significantly to the total failure probability.In cases where the Monte Carlo simulation does not detect failure domains with the provided computational budget, we therefore switch to subset simulation (Au and Beck, 2001) that is known to effectively estimate very small failure probabilities (Moustapha et al., 2021).

Applications
In this section, the proposed SSER algorithm is applied to four problems of varying complexity.These applications were selected to benchmark the algorithm performance in different types of reliability problems: the first two applications have analytical limit-state functions that are challenging for many standard reliability methods.They are the well-known four-branch function with two input parameters and four distinct failure regions, and a twodimensional piecewise linear reliability problem that is known to be challenging for subset simulation.The third and fourth applications are engineering models that involve finite element simulations.The third application is an engineering model of a five-story frame with M = 21 mutually dependent input parameters.Lastly, the fourth application is a continuum mechanical problem of a plate with a hole, where the Young's modulus is parametrized by a high dimensional (M = 869) random field.
All computations were done in MATLAB with the UQLab uncertainty quantification framework Marelli and Sudret (2014).Specifically, we used the input module Lataniotis et al. ( 2021) for modeling the random inputs, the PCE module Marelli et al. (2021) for the domain-wise spectral expansions with bootstrap replications and the reliability module Marelli et al. (2021) for computing the conditional failure probabilities.
A summary of the results and comparisons to other state-of-the-art reliability algorithms is given in Table 1, and a detailed rundown of each is provided in the following.Table 1: Applications summary: if not given with an external reference, the values were computed with UQLab Marelli and Sudret (2014); Marelli et al. (2021).Values in square brackets correspond to the 95% confidence intervals.N denotes the number of limit-state evaluations needed to compute the reported values.For SSER confidence intervals are also given for the number of limit-state evaluations until convergence from a set of independent runs.The used methods are Monte Carlo simulation (MCS), subset simulation (SuS), adaptive kriging Monte Carlo simulation with U learning function (AK-MCS+U), sparse polynomial chaos expansions with Monte Carlo simulation (sPCE+MCS), improved CE method with failure-informed dimension reduction (iCEred).
a Estimated from 100 independent SuS runs with reported CV Pf = 0.215 as Pf • (1 ± 2CV Pf / √ 100).b Reported value Pf = 3.9 • 10 −7 does not lie within reference confidence intervals.c Additional 675 gradient evaluations.d From independent algorithm runs to reach convergence.confidence intervals from the mean predictions of the independent SSER runs.mal distribution, i.e., X ∼ N (0, 1 2 ), and the limit-state function is given analytically by We obtain a reference failure probability of P f = 3.2 • 10 −5 (β = 4) for this problem from Monte Carlo simulation with N = 10 8 samples.
We perform 50 independent runs of SSER, with different random seeds.The number of refinement samples was set to N ref = 40 and the maximum polynomial degree for the residual expansions to p max = 6.The domain-wise accuracy was estimated with B = 500 bootstrap replications.
The convergence of SSER to the reference reliability index is shown in Figure 7.As intended by this limit-state function, most runs of the algorithm initially overestimate the reliability index.Due to its refinement domain selection strategy, however, SSER later returns to the high failure probability domains and accurately estimates the reliability index.
The median number of steps required to reach the stopping criterion is 5, corresponding to     In the following steps, it focuses again on the big region that incorporates the larger failure domain and accurately approximates it at Step 4. Two steps later, at Step 6, the stopping criterion is fulfilled and SSER is terminated.

Five-story frame
In this example we consider the five-story structural frame described in Figure 9.This example has been investigated before in Liu and Der Kiureghian (1991); Wei and Rahman (2007); Blatman and Sudret (2010).It consists of eight different structural elements with uncertain dimensions and material properties.The structure is subject to uncertain, horizontal loads P 1 , P 2 and P 3 .
The input parameters are collected in a random vector X = {P 1 , • • • , A 21 } and sorted according to Figure (c).These parameters are not mutually independent.Instead, we impose a Gaussian copula (Nelsen, 2006) such that the joint input distribution function can be written as   confidence intervals from the mean predictions of the independent SSER runs.
Monte Carlo simulation and N = 10 8 evaluations of the limit-state function.
We perform 50 independent runs of SSER on this problem, with different random seeds.
The number of refinement samples is set to

Plate with a hole
This last example was previously analysed in Uribe et al. (2020) and its description here follows closely the description there.Consider a two-dimensional, square steel plate under plane stress conditions with a hole in the middle (see Figure 11).The plate is defined on forces results in the following system of partial differential equations (PDEs) to the Karhunen-Loève expansion Ghanem and Spanos (1991) and write where λ k and ϕ k (ξ) are the ordered (i.e., λ k ≥ λ k+1 ) eigenvalues and eigenfunctions of the covariance operator C(ξ, ξ The mean and variance of the underlying Gaussian field are µ E and σ 2 E respectively and θ ∈ R k denotes the standard Gaussian vector of random coefficients.The expansion is truncated at K = 868 terms, which account for 92.5% of the spatial average of the variance of the Gaussian random field.The eigenpairs are estimated with the Nyström method based on 100 Gauss-Legendre points in each direction, disregarding points inside the hole. Additionally, we model the applied external load as a normally distributed random variable with mean µ q = 60 MPa and standard deviation σ q = 12 MPa such that q ∼ N (µ q , σ q ).Assuming further independence between the standard normal vector θ and q, the probabilistic input vector of this problem is gathered in X = (θ, q).The quantity of interest is the principal stress at a control point marked by a red star in Figure 11 defined as It depends on the stress field σ(ξ) at the control point that is obtained from the elastic constitutive equations after solving Eq. ( 28) with the finite element method for the displacement field u(ξ).Failure in this problem is defined as this principle stress exceeding σ th = 320 MPa, which can be formalized in the following limit-state function: Due to the complexity of this problem, it is not possible to obtain a reference solution with Monte Carlo simulation in a reasonable time.We use instead directly the reference solution from Uribe et al. ( 2020) that was computed with 100 independent subset simulation (Au and Beck, 2001) runs with N = 3,000 samples per level and 6 levels, leading to a reference failure probability of P f = 3.75 • 10 −6 (β = 4.48).
On this problem, we perform 10 independent runs of SSER with different random seeds.
The number of refinement samples was set to N ref = 100 and the maximum polynomial degree for the residual expansions to p max = 4.To allow PCEs to be constructed efficiently in such high dimensions, it was necessary to restrict the maximum interactions of polynomial basis functions to 1 (Marelli et al., 2021).The domain-wise accuracy was estimated with B = 100 bootstrap replications.
The convergence of SSER to the reference reliability index is shown in Figure 12.The median number of steps required to reach the stopping criterion is 6, corresponding to

Conclusions
Leveraging the flexibility of the recently proposed stochastic spectral embedding formalism (Marelli et al., 2021), we show that the adaptive sequential partitioning approach introduced in (Wagner et al., 2021) can be efficiently modified to an active learning reliability method.
The introduced modifications pertain to the refinement domain selection, partitioning and sample enrichment strategies.
The SSER algorithm is shown to accurately estimate failure probabilities spanning multiple orders of magnitude at a fraction of the cost of plain Monte Carlo simulation.The  confidence intervals from the mean predictions of the independent SSER runs.
method is suitable for identifying multiple failure regions (four-branch function, Section 3.1) and is able to deal with problems that involve non-smooth limit-state functions (piecewise linear function, Section 3.2).It also performs reasonably well in moderate-dimensions with dependent input vectors (five-story frame, Section 3.3 and is competitive to a state-of-the-art reliability method iCEred (Uribe et al., 2020) on a high-dimensional example (plate with a hole, Section 3.4).
In general, the efficiency of SSER depends on two factors: (1) whether the algorithm identifies all failure modes and (2) whether the chosen spectral expansion technique can efficiently approximate g at the terminal domains that lie in proximity of the limit-state surface.Regarding the first factor, the refinement domain selection criterion in Eq. ( 15) should guarantee that SSER eventually identifies all modes of failure (e.g., four-branch function and piecewise linear function).However, in problems with a large number of disjoint limit-states, this sequential process might require a considerable number iterations and consequently model evaluations.Also, if the stopping criterion is triggered before all limit-states have been identified, the total failure probability might be underestimated.The second factor depends on g's complexity at the terminal domains in proximity of the limit-state surface.
Generally, spectral expansion techniques are known to converge slowly in highly non-linear and high-dimensional problems.Hence, we expect a decrease in efficiency of the algorithm in such problems.As an example, SSER is expected to perform poorly in first excursion problems that are both high-dimensional and highly non-linear in the limit state-function in proximity of the limit-state surface.
In some examples (most notably in Section 3.4) there persists a small bias towards a higher reliability index compared to the reference solution.This means that the method systematically misses small parts of failure regions.A possible explanation for this is the partitioning strategy that 'unknowingly' cuts away undiscovered parts of failure regions and incorrectly classifies them as safe.This process is facilitated by the boundary discontinuities and will be explored in future works.
An interesting side result of the proposed method is the domain partition that hints at the location of failure domains.This information, in conjunction with the available conditional failure probabilities, could be directly used for design purposes.
The modularity of the sequential partitioning algorithm for constructing SSEs makes it well suited to for future modifications.Possible improvements could focus on different spectral expansion or even surrogate modelling techniques, as well as simultaneous partitioning in multiple dimensions.A particularly interesting idea is the introduction of rotated domain partitions by means of partial least squares Papaioannou et al. (2019).These partitions could be more flexibly tuned to fit the limit-state surface and thereby increase the efficiency of SSER.

Figure 1 :
Figure 1: Stochastic spectral embedding: Illustration of SSE representation in Eq. (3) in the quantile and real space.Red points denote the experimental design and orange areas denote terminal or unsplit domains, that partition the input space into a set of disjoint subdomains such that D X = k∈T D k X .
For each subdomain: i. Enrich sample ii.Create residual expansion iii.Update domain-wise error

Figure 2 :
Figure 2: Stochastic spectral embedding: Illustration of how dependence is handled in the algorithm.For two mutually dependent input parameters X 1 and X 2 , the plots show a set of sample points in the real and quantile space with a subset of points inside an exemplary subdomain.The parameter dependence leads to a loss of orthogonality of the subdomains in the real space, which necessitates the introduction of an enveloping hypercube in the real space to construct the spectral basis functions Ψ k α .

Figure 3 :
Figure 3: Partitioning strategy exemplified on a one dimensional function showing the bootstrap replications and the resulting support of Z + and Z 0 as well as the objective function L. The splitting location υ is highlighted as a vertical, dashed red line.The plots are shown in the quantile space D U .

2.
For ( , p) = k refine from Eq. (15) (a) Split the current subdomain D ,p X in 2 sub-parts D +1,{s1,s2} X according to Eq. (19) , where s 1 and s 2 are unique indices for domains on the -th level , and update T (b) For each split s ∈ {s 1 , s 2 } i. Enrich sample X with min{N ref , N ED − |X |} new points inside D +1,s X ii.If N ref new points inside X A. Create the truncated expansion R +1,s S (X) of the residual R +1 (X) in the current subdomain using a subset of X inside D +1,s X B. Update the residual on the next level in the current subdomain R +2 (a) Return Pf and its variance A graphical sketch of the algorithm for a two-dimensional toy problem is shown in Figure 4.The figure shows the first two iterations of the SSE construction in the quantile and real space, highlighting the terminal domains T at each step in orange.The initialization step in Figure (a) shows the initial experimental design and global domain D 0,1 X and their counterpart in the quantile space D 0,1 U .At the first iteration, the algorithm chooses the only available refinement domain k refine = (0, 1), partitions it into two subdomains, enriches the experimental design and constructs two additional residual expansions.At the second iteration in Figure (c), the algorithm refines D (1,2) X out of two terminal domains T = {(1, 1), (1, 2)}, because it has the largest contribution to the total failure probability variance.Refinement then continues until the stopping criterion is met.

Figure 4 :
Figure 4: Graphical representation of the first steps of the active learning algorithm described in Section 2.2.4 for a two-dimensional toy problem with independent inputs.Upper row: partitioning in the quantile space; lower row: partitioning in the unbounded real space with f X contour lines in dashed blue.Red dots show the adaptive experimental design with N ref = 4.The terminal domains T are highlighted in orange.The red area denotes the failure domain.In this illustrative example the failure probability is Pf ≈ 4.2 • 10 −2 .

Figure 5 :
Figure 5: Four-branch function: Convergence to the reference reliability index β as a function of limit-state evaluations and algorithm steps.The 95% confidence intervals in Figure (a) are based on the bootstrap replications.This figure shows the convergence until the stopping criterion in Eq. (23) is met at Step 9 of SSER.The bounds in Figure (b) are the 90%, 75% and 50% 2N ref initial +(4 • 2)N ref = 400 limit-state evaluations.

Figure 8
Figure 8 reveals why SSER initially overestimates the reliability index and underestimates

Figure 6 :
Figure 6: Four-branch function: Misclassification of failure/safe points for selected steps of SSER.The white areas were classified correctly, while the gray areas were misclassified.The thick black line indicates the limit-state surface.The domain bounds identified by the SSER are shown as thin black lines.The black points are the used experimental design and the red points are the points added at the shown step of the algorithm.The convergence of the reliability index estimator is shown in Figure (a).

Figure 7 :
Figure 7: Piecewise linear function: Convergence to the reference reliability index β as a function of limit-state evaluations and algorithm steps.The 95% confidence intervals in Figure (a) are based on the bootstrap replications.This figure shows the convergence until the stopping criterion in Eq. (23) is met at Step 6 of SSER.The bounds in Figure (b) are the 90%, 75% and 50% confidence intervals from the mean predictions of the independent SSER runs.

)
In this expression FXi are the cumulative distribution functions (CDF) of the parameters as listed in Figure (c), Φ 21 (U ; R) is the CDF of a 21-variate Gaussian distribution with

Figure 8 :
Figure 8: Piecewise linear function: Misclassification of failure/safe points for selected steps of SSER.The white areas were classified correctly, while the gray areas were misclassified.The thick black line indicates the limit-state surface.The domain bounds identified by SSER are shown as thin black lines.The black points are the used experimental design and the red points are the points added at the shown step of the algorithm.The convergence of the reliability index estimator is shown in Figure (a).

Figure 10 :
Figure 10: Five-story frame: Convergence to the reference reliability index β as a function of limit-state evaluations and algorithm steps.The 95% confidence intervals in Figure (a) are based on the bootstrap replications.This figure shows the convergence until the stopping criterion in Eq. (23) is met at Step 6 of SSER.The bounds in Figure (b) are the 90%, 75% and 50% 40 and the maximum polynomial degree for the residual expansions to p max = 4.The domain-wise accuracy is estimated with B = 500 bootstrap replications.The convergence to the reference reliability index for a single run of SSER and all 50 runs is shown in Figure 10.The median number of steps required to reach the stopping criterion is 3, corresponding to 2N ref initial +(2 • 2)N ref = 240 limit-state evaluations.The stopping criterion for the shown run in Figure (a) is met at Step 6.

aFigure 11 :
Figure 11: Plate with a hole: Finite element discretization with 282 eight-node serendipity quadrilateral elements (Uribe et al., 2020).The red star * marks the location of the control node for the stress exceedance criterion in Eq. (31).
) where u(ξ) is the 2D displacement field and ν = 0.29 is the steel Poisson's ratio.As sketched in Figure 11 a Dirichlet boundary condition is applied on the left domain border such that u(ξ) = 0 for ξ ∈ Γ 1 .Additionally, a Neumann boundary condition is applied on the right domain border such that σ 11 (ξ) = q for ξ ∈ Γ 2 , where σ 11 is the stress component in the first direction.The Young's modulus E(ξ) is modelled as lognormal random field with mean µ E = 2 • 10 5 MPa and standard deviation σ E = 3 • 10 4 MPa.We use an isotropic exponential kernel k(ξ, ξ ) = exp (− −1 ||ξ − ξ || 2 ) as the autocorrelation function of the underlying Gaussian random field with correlation length = 0.04 m.To parametrize this random field, we resort

Figure 12 :
Figure 12: Plate with a hole: Convergence to the reference reliability index β as a function of limit-state evaluations and algorithm steps.The 95% confidence intervals in Figure (a) are based on the bootstrap replications.This figure shows the convergence until the stopping criterion in Eq. (23) is met at Step 7 of SSER.The bounds in Figure (b) are the 90%, 75% and 50% To account for this, we start out with an initial experimental design of size 2N ref and switch to the experimental design size of N ref in all subsequent iterations.
Finally, we note that N ref is a tuning parameter of the algorithm.Larger values of N ref lead to more accurate local expansions, resulting in increased stability of the algorithm at the cost of slower convergence rates.Vice-versa, smaller values of N ref can significantly increase the convergence rates, but can potentially produce inaccurate limit-state function approximations that lead to biased failure probability estimates.