#  Reliability analysis of discrete-state performance functions via adaptive
sequential sampling with detection of failure surfaces

Miroslav Vořechovský  Corresponding author.
Institute of Structural Mechanics, Brno University of Technology,
Veveří 331/95, 602 00 Brno, Czech Republic,
e-mail<EMAIL_ADDRESS>
###### Abstract

The paper presents a new efficient and robust method for rare event
probability estimation for computational models of an engineering product or a
process returning categorical information only, for example, either success or
_failure_. For such models, most of the methods designed for the estimation of
failure probability, which use the numerical value of the outcome to compute
gradients or to estimate the proximity to the failure surface, cannot be
applied. Even if the performance function provides more than just binary
output, the state of the system may be a non-smooth or even a discontinuous
function defined in the domain of continuous input variables. This often
happens because the mathematical model features non-smooth components or
discontinuities (e.g., in the constitutive laws), bifurcations, or even
domains in which _no reasonable model response_ is obtained. In these cases,
the classical gradient-based methods usually fail. We propose a simple yet
efficient algorithm, which performs a sequential adaptive selection of points
from the input domain of random variables to extend and refine a simple
distance-based surrogate model. Two different tasks can be accomplished at any
stage of sequential sampling: (i) estimation of the failure probability, and
(ii) selection of the best possible candidate for the subsequent model
evaluation if further improvement is necessary. The proposed criterion for
selecting the next point for model evaluation maximizes the expected
probability classified by using the candidate. Therefore, the perfect balance
between _global exploration_ and _local exploitation_ is maintained
automatically. If there are more rare events such as failure modes, the method
can be generalized to estimate the probabilities of all these event types.
Moreover, when the numerical value of model evaluation can be used to build a
smooth surrogate, the algorithm can accommodate this information to increase
the accuracy of the estimated probabilities. Lastly, we define a new simple
yet general geometrical measure of the global sensitivity of the rare-event
probability to individual variables, which is obtained as a by-product of the
proposed refinement algorithm.

###### keywords:

$\psi$ criterion , failure probability , failure surface , Gaussian space ,
extension of experimental design, gradient-free optimization , categorical
limit state function , binary surrogate model , importance sampling

††journal: Computer Methods in Applied Mechanics and Engineering

accurate estimation of failure probability with error prediction, resistance
to multiple design points, and strong nonlinearities

the method needs at least a categorical limit state function (including no
response output) for continuously distributed inputs

adaptive sequential extension of experimental design using the $\psi$
criterion, minimizing the number of limit state evaluations

the novel criterion automatically balances between exploration and
exploitation in the input space

efficient importance sampling in important rings, excluding the safe ball and
unimportant exterior

global sensitivity indices quantifying the importance of individual variables
to rare event probability

adaptive expansion and refinement near the failure surface avoiding clusters
of points (no redundant information)

##  1 Introduction

As summarized in a 2001 review by Rackwitz [1], the issue of the determination
of failure probability attracted considerable attention in the second half of
the last century. Currently, the interest of industry and researchers shows
that the importance of the topic is even increasing. Freudenthal, whose
landmark paper [2] dates back to 1947, is often considered the father of
modern structural reliability. The standard problem setting is a limit state
function $g(\boldsymbol{x})$ (also called a performance function or
computational model) with continuously distributed random inputs
${\boldsymbol{X}}$ forming the random vector with dimension
${N_{\mathrm{var}}}$. The probability density function
$f_{{\boldsymbol{X}}}({\boldsymbol{x}})$ is assumed as known. The spectrum of
existing methods for reliability estimation is very rich and they can be
classified based on the way how and how much they use and also make targeted
modifications to the three groups of input input information (i) the
information about the limit state function $g(\boldsymbol{x})$, (ii) the
density of the input vector ${\boldsymbol{X}}$, and (iii) the geometry and
topology of the input space. Crude Monte Carlo sampling estimation of
probabilistic integrals does not modify or make special use of any of these
inputs which makes it very robust method, but also quite inefficient. Other
methods utilize more or less one or more pieces of information and also make
various kinds of assumptions about these inputs to improve their efficiency.

The development of the mathematical theory of the first order reliability
method (FORM) [3, 4] and the second-order reliability theory (SORM) [5, 6, 7,
8] with the related FORM/SORM asymptotic approximations for multinormal
integrals using Taylor series expansions of the first/second order [9, 10] was
completed in the eighties [11, 12, 13]. These foundational methods are still
perhaps the most important results on which modern design codes for engineers
are built. In 1983, two papers [14, 15] brought importance sampling (IS)
concepts to the attention of the reliability community. The great advantage of
having error estimation alongside probability estimation made IS an excellent
tool for reliability estimation updating, and it also improved FORM/SORM
estimations [16, 17, 18]. Importance sampling around the mean value and the
method of asymptotic sampling [19] both exploit the fact that in most
engineering reliability problems the failure regions are located further from
the central region of the joint density, and the sampling density increases
the spread compared to the joint density. Another version of importance
sampling used a sampling density centered on the design point, which expresses
the premise that the most central failure point is surrounded by a small
neighborhood, which contains most of the failure probability $p_{\pazocal{F}}$
(a result derived by Breitung and Hohenbichler [7]; see also Breitung [10]).
Indeed, a design point in the space of independent standard normal variables
is a failure point that maximizes the joint normal density, or equivalently,
minimizes the distance to the origin. Moving away from this point into the
failure domain as well as moving along the boundary between failure and safe
regions decreases the probability density. These points are searched based on
the premise that the value of the performance function describing the state of
the analyzed product or process decreases towards them in the safe region and
becomes negative in the failure region. FORM usually operates in the standard
normal space where the “design points” have to be found.

As pointed out by Rackwitz [1], in practical engineering the _reliability
analysis_ of an engineering product (be it a structure, a mechanical system,
or a process) is itself perhaps less interesting than the _optimization_ [20]
of those products or processes with reliability constraints. However,
reliability optimization methods encompass reliability analysis and call it
repeatedly [21, 22].

With ever-growing computational power, numerical tools are now, and more than
ever, being widely applied to the representation and solution of complex
problems in engineering and the sciences. Unfortunately, in many cases, the
problem of evaluating reliability via the estimation of failure probability
remains intractable. This is often because of the large computational expenses
associated with the model evaluations that are needed to run a selected
sampling method for failure probability estimation, such as importance
sampling [23, 24, 25, 26], line sampling [27, 28], directional simulation [29,
30], asymptotic sampling [19] or subset simulation [31, 32, 33]. High-fidelity
computational models describing the performance of engineering systems are
computationally intensive, and the application of advanced (and often
adaptive) sampling methods in combination with a high-fidelity model to obtain
sufficiently accurate estimations of failure probability is not feasible.

In order to build an approximate representation of a true performance
function, which is inexpensive to analyze via sampling strategies, a variety
of surrogate models have been developed. Smart strategies exploit the
information from already analyzed points from the design domain, and they can
adaptively refine the surrogate model; see, e.g. [34, 35]. Many _adaptive_
reliability analysis methods have been developed, some using radial basis
functions [36, 37], support vector regression surrogates [38, 39, 40, 41, 42],
artificial neural networks [43, 44, 45], sparse polynomial chaos [46, 47] or
Kriging [48, 49]. Another alternative for probabilistic calibration as new
information emerges is to use Bayesian updating with Kriging [50]. Authors
have invented various combinations of methods, such as a combination of
Importance Sampling and Kriging [51]. Unfortunately, most of these methods are
developed with the assumption that the surrogate approximates a well-behaved
function defined over the whole design space. Moreover, methods based on
Kriging (a.k.a. Gaussian process regression) are often criticized for the
strong dependence of the results on the selection of the kernel, which is a
user choice, often subjective, and hardly supported by sufficient data.
Recently, the sample space partitioning strategy has been proposed along with
an adaptive Kriging model [52], which relaxes this assumption. However, most
adaptive strategies often use various kinds of learning functions in which the
numerical value of the performance function plays an important role.

In most of the existing methods, it is assumed that the performance of a
system is a smooth function $g(\boldsymbol{x})$ defined over the whole input
domain $\pazocal{D}$. Traditionally, the performance function
$g(\boldsymbol{x})$ is defined such that its negative values signal a
_failure event_ and positive values success. The boundary between the safe and
failure sets is called the _failure surface_ and will be denoted as
$\partial{\pazocal{F}}$. We assume it is an
$({N_{\mathrm{var}}}-1)$-dimensional object, and it is often illustrated as
the zero-valued contour of $g(\boldsymbol{x})$. When the failure surface
$\partial{\pazocal{F}}$ is a differentiable function and the function is
almost linear in the standard normal space, the failure probability has a
simple relationship to the shortest distance $\beta$ from the origin to the
failure surface: ${p_{\pazocal{F}}}\approx\Phi(-\beta)$. The rotational
symmetry of the standard normal density makes this FORM approximation simple
because the failure domain is approximated by a half-space fitted to the true
failure domain at the distinct _nearest failure point_. Many optimization
algorithms have been developed [53] to find the nearest failure point for
smooth failure surfaces. The most important is the method by Hasofer and Lind
[3] for second-moment reliability analysis, which Rackwitz and Fiessler [4]
later extended to include distribution information (the famous “HL-RF”
iterative scheme). We remark that while this point marks the failure event
with the maximum Gaussian density, it may not be the most likely failure point
in terms of the original non-Gaussian density due to the nonlinearity of the
probability transformation. The main task then becomes the location of the
most central failure point (often called the design point or the $\beta$
point).

On the other hand, methods guided by the numerical values of the performance
function might be misled into heading in an incorrect search direction. Not
only the “design point search” needed for FORM and SORM relies on the
assumption that the performance function decreases towards the most central
failure point. Subset Simulation (SuS) [31] or some sequential adaptive
variants of Importance Sampling methods [17, 26] are also based on this
assumption. In this way, SuS, for instance, can be viewed as a stochastic
version of the gradient (downhill) optimization method and, as exemplified by
Breitung [54] using numerous examples, the method can fail for functions with
complicated landscapes. Similar to the design point search for FORM, such
downhill optimization methods may not be successful in global optimization
because the search may be caught in local minima. Repeated SuS runs may not
help if the search is always initiated from the origin and the information
about the geometry of the performance function is not stored. The problem of
multiple design points has not been satisfactorily solved even in FORM.
Moreover, various different formulations of $g(\boldsymbol{x})$ that lead to
the same failure set may alter the result of $p_{\pazocal{F}}$ estimation
because the formulations change the evolution of the estimation process, which
is obviously wrong.  The idea of extrapolation with a sequence of modified
problems, which was focused on modifications of the sampling density in
Asymptotic Sampling [19], can be seen in an analogy to making modification to
the limit state function [55, 56], thus making assumptions about the role of
supposedly smooth landscape of $g(\boldsymbol{x})$. A different group of
reliability methods which can be termed “moment methods” also use the
numerical values of $g(\boldsymbol{x})$: the basic idea is to fit a proper
probability distribution to the output variable of the limit state function
based on the knowledge of its estimated moments of certain type (see, e.g.,
integer moment based methods [57, 58], fractional moment based methods [59,
60, 61], moment-generating function (or Laplace transform) based methods [62,
63, 64]).

Moreover, it sometimes happens that complex computational models are not able
to provide any answer for some input values (e.g., a nonlinear finite element
solver of a structure is not able to converge for a combination of input
parameters). Or, the performance of a system is not a smooth function or
contains discontinuities which pose a problem for most of the gradient-based
algorithms embedded in the reliability analyses. The gradient-based algorithms
are used to locate the most central failure point (or more of these points for
different parts of the design domain), and the shape of the supposedly
continuous performance function is used to estimate the distance to the
boundary between failure and safe regions (the safety margin). However, some
performance functions of a system may be just a discrete number of states or
just a binary function returning either “success” or “failure”. In these
cases, the majority of the above-mentioned algorithms and methods fail
entirely, and thus the motivation behind the present work was to develop a
robust technique that can solve all these problems related to the performance
function.

There is thus a need for a method that is resistant to noise or a non-smooth
shape in the performance function, its jumps or even discrete values, and yet
is able to provide a reasonable failure probability estimation with a small
number of function evaluations. The method should balance between the global
_exploration_ of the design domain in search for new input space territories
leading to failure and the local _exploitation_ of the previously discovered
boundary between safe and failure sets (failure surface) in order to refine
its description. We see it as the problem of “how to divide a territory?": to
detect and geometrically describe the boundaries between the safe and failure
sets. Optimally, the method should sequentially extend the experimental design
(ED) in single steps to maximize the gain from the already obtained
information at each stage. The method should keep adapting its representation
of the true performance via some form of a surrogate model, which is
inexpensive to evaluate and which uses a tailored sampling strategy for quick
on-the-fly reliability estimation. The method should not make unnecessary
calls of the expensive response function in regions surrounded by “safe”
samples (no matter “how” safe they are) as these are almost sure to be also a
part of the safe set. Analogously, we want to avoid unnecessary new samples in
the (almost sure) failure regions. It is evident that the most precious
information is the refinement of the boundary separating safe and failure
regions in proportion to the local density. It is there that the expensive
information needs to be obtained by calling the true performance function,
i.e., the limit state function.

##  2 Problem statement

Let us assume a ${N_{\mathrm{var}}}$-dimensional vector ${\boldsymbol{X}}$ of
continuously distributed basic random variables with known joint probability
density function $f_{{\boldsymbol{X}}}({\boldsymbol{x}})$. Vector
${\boldsymbol{X}}$ is the input to the performance function of a
system/process whose reliability is to be evaluated. Assume that for any
realization (a point from the _design domain_ , $\pazocal{D}$),
${\boldsymbol{x}}\in\pazocal{D}$, we can obtain, albeit at high expense,
information as to whether or not the system fails, or generally about the
system performance. We define the failure set (domain)
${\pazocal{F}}\subset\pazocal{D}$ as a union of all regions within the design
domain in which failure occurs. The probability of failure $p_{\pazocal{F}}$
is then defined as

$\displaystyle{p_{\pazocal{F}}}\equiv\idotsint_{{\pazocal{F}}}f_{{\boldsymbol{X}}}({\boldsymbol{x}})\;\mathrm{\,d}\hskip
0.50003pt{\boldsymbol{x}}=\idotsint_{\pazocal{F}}\mathrm{\,d}\hskip
0.50003ptF_{{\boldsymbol{X}}}({\boldsymbol{x}}),$ (1)

where $F_{{\boldsymbol{X}}}({\boldsymbol{x}})$ is the cumulative density
functions of the random vector ${\boldsymbol{X}}$. The integral over the whole
$\pazocal{D}$ equals one (exhaustion of all possible events).

Let us assume that the performance of a system is evaluated via a performance
function, sometimes referred to as the _limit state function_
$g(\boldsymbol{x})$. From here on, we will consider the crudest case in which
$g(\boldsymbol{x})$ returns either 1 (failure) or zero for a safe state (cases
when the performance function does not provide any answer in some input points
will be discussed separately). Limiting the response to a  _binary function_
does not affect the definition of $p_{\pazocal{F}}$ in Eq. (1) as it never
included any other information. When working with systems for which we can
obtain responses other than a binary response via $g(\boldsymbol{x})$, we use
it to define the _indicator function_
${\boldsymbol{1}_{\pazocal{F}}\left({\boldsymbol{x}}\right)}$, which returns
one for ${\boldsymbol{x}}$ falling to the failure set (typically
$g(\boldsymbol{x})<0$) and zero otherwise

$\displaystyle{\boldsymbol{1}_{\pazocal{F}}\left({\boldsymbol{x}}\right)}$
$\displaystyle=\left\\{\begin{array}[]{ll}1,&{\boldsymbol{x}}\in{\pazocal{F}}\text{
(failure event)},\\\ 0&\text{otherwise}.\end{array}\right.$ (4)

Note that the failure set may be composed of more types of failure or
indications of malfunctioning, and the corresponding failure probabilities may
also be desired. We address this eventuality by generalizing the dichotomous
nature to handle more than two event types. A rare event will be generally
denoted as $\pazocal{T}$. The method described in this paper focuses on the
approximation of the boundaries between various sets. However, as illustrated
below, we assume that the rare-event sets are not scattered over too many
disjoint sets, which are closed and have individually almost zero measure.

By including this indicator function in the integrand in Eq. (1), the integral
over the whole design space gets narrowed down to the failure domain

$\displaystyle{p_{\pazocal{F}}}=\idotsint_{\pazocal{D}}{\boldsymbol{1}_{\pazocal{F}}\left({\boldsymbol{x}}\right)}\,f_{{\boldsymbol{X}}}({\boldsymbol{x}})\,\mathrm{\,d}\hskip
0.50003pt{\boldsymbol{x}}.$ (5)

Without loss of generality, we assume that the joint probability distribution
function (PDF) of ${\boldsymbol{X}}$ is the ${N_{\mathrm{var}}}$-dimensional
standard normal density with independent components. If the original problem
features non-normal marginals and dependencies, we assume that a probability
preserving transformation exists, such as the Nataf model [65, 66, 67, 68]
(sometimes referred to as the Gaussian copula [69]) or the Rosenblatt-
transformation, via a product of conditional distributions [70, 13]. The
Gaussian space is particularly suitable for the proposed method as it is
rationally invariant and open (unbounded). We remark that in many practical
problems, the information about the joint PDF is often not complete. If the
information is limited to univariate marginals and the correlation matrix
only, there exist many possible forms of the joint PDF fulfilling this input,
and it may happen that none of them can be represented via the Nataf
transformation. New ways to transform the original joint PDF to a Gaussian
space which are more flexible [68] exist, and are still appearing. They may
employ polynomial chaos representation [71] or other distributions [72]. We
deliberately refrain from discussing these transformations further as they are
not the main focus of this paper.

This paper proceeds with a quick presentation of tools useful in rotationally
symmetrical standard Gaussian space (Sec. 3), and these are immediately used
to present the proposed sequential adaptive extension of sample size in Sec.
4. The desired failure probability can be estimated by a quick and simple
evaluation of the information obtained at the current ED via binary surrogate,
which we propose should be constructed on-the-fly; see Sec. 5. The byproduct
of this estimation is the information about novel global sensitivity indices
which we propose in Sec. 6.

##  3 Useful tools in the standard Gaussian space

We assume that the probability density in the standard Gaussian space is
jointly Gaussian, i.e., the individual standard normal random variables are
independent. Therefore, the probability density function of any point
${\boldsymbol{x}}=\\{x_{1},\ldots,x_{{N_{\mathrm{var}}}}\\}$ becomes the
product of univariate Gaussian densities of the individual marginals

$f_{{\boldsymbol{X}}}({\boldsymbol{x}})=\prod_{v=1}^{{N_{\mathrm{var}}}}\varphi(x_{i})=\frac{1}{(2\pi)^{{N_{\mathrm{var}}}/2}}\exp{\left(-\frac{1}{2}{\boldsymbol{x}}^{\mkern-1.5mu\mathsf{T}}{\boldsymbol{x}}\right)},$
(6)

where $\varphi(x)={(2\pi)^{-1/2}}\exp{\left(-{x^{2}}/{2}\right)}$ is the
standard univariate Gaussian density.

Due to the rotational symmetry with respect to the origin, the space can be
indexed using (i) ${N_{\mathrm{var}}}-1$ independent directions and (ii) one
univariate Euclidean distance from the origin

$\rho=\lVert{\boldsymbol{x}}\rVert=\sqrt{{\boldsymbol{x}}^{\mkern-1.5mu\mathsf{T}}{\boldsymbol{x}}}=\sqrt{\sum_{v=1}^{{N_{\mathrm{var}}}}x_{v}^{2}},\quad
v=1,\ldots,{N_{\mathrm{var}}}.$ (7)

Therefore, the standard Gaussian density
$f_{{\boldsymbol{X}}}({\boldsymbol{x}})$ at a point ${\boldsymbol{x}}$
simplifies to univariate Gaussian density depending on the radial distance
$\rho$ only

$f_{{\boldsymbol{X}}}({\boldsymbol{x}})=\varphi(\rho)=\frac{1}{(2\pi)^{{N_{\mathrm{var}}}/2}}\exp{\left(-\frac{\rho^{2}}{2}\right)}$
(8)

Figure 1:  Left: probability density and distribution function for a random
radial distance in standard Gaussian space. Middle: radial distance density of
a point selected outside an $N_{\mathrm{var}}$-ball of radius $r$. Right:
radial distance density of a point selected in between two
$N_{\mathrm{var}}$-balls of radii $r$ and $R$.

An important ingredient in the proposed algorithm is the distribution function
for a random distance $\rho$ in the Gaussian space. Assume a randomly selected
point, ${\boldsymbol{x}}$. Its Euclidean distance from the origin has $\chi$
(chi) distribution with ${N_{\mathrm{var}}}$ degrees of freedom. The
probability density function of a distance $\rho>0$ reads (see Fig. 1 bottom
left)

$f_{\rho}(\rho;{N_{\mathrm{var}}})=\frac{2^{1-{N_{\mathrm{var}}}/2}}{\Gamma\left(\frac{{N_{\mathrm{var}}}}{2}\right)}\rho^{{N_{\mathrm{var}}}-1}\exp{\left(-\frac{\rho^{2}}{2}\right)},$
(9)

where $\Gamma(\cdot)$ is the standard (complete) gamma function. The
corresponding cumulative distribution function, defined as
$F_{\rho}(\rho;{N_{\mathrm{var}}})=\int_{0}^{\rho}f_{\rho}(t,{N_{\mathrm{var}}})\mathrm{\,d}\hskip
0.50003ptt$, reads (see Fig. 1 top left)

$\displaystyle
F_{\rho}(\rho;{N_{\mathrm{var}}})=\frac{2^{1-{N_{\mathrm{var}}}/2}}{\Gamma\left(\frac{{N_{\mathrm{var}}}}{2}\right)}\int_{0}^{\rho}t^{{N_{\mathrm{var}}}-1}\exp{\left(-\frac{t^{2}}{2}\right)}\mathrm{\,d}\hskip
0.50003ptt=P\left(\frac{{N_{\mathrm{var}}}}{2},\frac{\rho^{2}}{2}\right),$
(10)

where $P(s,x)=\gamma(s,x)/\Gamma(s)$ is the _regularized lower incomplete
gamma function_. It follows from the lower incomplete gamma function
$\gamma(s,x)=\int_{0}^{x}t^{s-1}\exp(-t)\mathrm{\,d}\hskip 0.50003ptt$ by its
regularization via the complete gamma function $\Gamma(s)$. $P(s,x)$ is one of
the “special functions” and it is a standard part of various mathematical
libraries, such as scipy [73], which is available in Python (function
gammainc($\cdot,\cdot$)).

Therefore, when a point is selected randomly from
${N_{\mathrm{var}}}$-dimensional standard Gaussian space, its distance from
origin has a distribution as if it were selected from a univariate $\chi$
distribution with $N_{\mathrm{var}}$ degrees of freedom. Its mean Euclidean
distance $\mu_{\rho}$, the approximate median (50 % percentile
$\mathrm{me}_{\rho}$), and the modus (the most probable distance
$\mathrm{mo}_{\rho}$) read

$\displaystyle\mu_{\rho}$
$\displaystyle=\sqrt{2}\frac{\Gamma\left(({N_{\mathrm{var}}}+1)/2\right)}{\Gamma({N_{\mathrm{var}}}/2)}\xrightarrow{{N_{\mathrm{var}}}\to\infty}\sqrt{{N_{\mathrm{var}}}-1/2},\quad\mathrm{me}_{\rho}\approx\sqrt{{N_{\mathrm{var}}}}\left(1-\frac{2}{9{N_{\mathrm{var}}}}\right)^{3/2}$
(11) $\displaystyle\mathrm{mo}_{\rho}$
$\displaystyle=\sqrt{{N_{\mathrm{var}}}-1}$ (12)

The variance of a random distance $\rho$ is simply
$\sigma^{2}_{\rho}={N_{\mathrm{var}}}-\mu_{\rho}^{2}$, which asymptotically
equals $1/2$, independently of ${N_{\mathrm{var}}}$.

###  3.1 Interior of an ${N_{\mathrm{var}}}$-ball in standard Gaussian space

An important geometrical entity in the proposed method is an
$N_{\mathrm{var}}$-ball with radius $r$, which is centered at the origin of
the coordinate system; see the green domain in Fig. 1 top middle. The volume
and the surface of the ${N_{\mathrm{var}}}$-ball, denoted as
$B_{r}\in\mathbb{R}^{{N_{\mathrm{var}}}}$, read

$\displaystyle\mathrm{Vol}\left[B_{r}\right]$
$\displaystyle=\frac{\pi^{{N_{\mathrm{var}}}/2}}{\Gamma\left(\frac{{N_{\mathrm{var}}}}{2}+1\right)}r^{{N_{\mathrm{var}}}},$
(13) $\displaystyle\mathrm{Sur}\left[B_{r}\right]$
$\displaystyle=\frac{2\pi^{{N_{\mathrm{var}}}/2}}{\Gamma\left(\frac{{N_{\mathrm{var}}}}{2}\right)}r^{{N_{\mathrm{var}}}-1}.$
(14)

We remark that the $\chi$ density presented in Eq. (9) can be obtained simply
by multiplying the point Gaussian density from Eq. (8) by the ball surface
from Eq. (14), for which we take the radius $r$ equal to the random radial
distance $\rho$, i.e.,
$f_{\rho}=\varphi(\rho)\cdot\mathrm{Sur}\left[B_{\rho}\right]$.

Given this probability density function and the corresponding cumulative
distribution function in Eq. (10), one can easily evaluate the probability
content of the ${N_{\mathrm{var}}}$-ball interior. The probability that a
random Gaussian point falls inside the ball $B_{r}$ reads:

$\displaystyle p_{\mathrm{in}}(r,{N_{\mathrm{var}}})\equiv
F_{\rho}(r;{N_{\mathrm{var}}})=P\left({N_{\mathrm{var}}}/2,r^{2}/2\right).$
(15)

The inverse of the regularized lower incomplete gamma function,
$P^{-1}(\cdot)$, implemented as gammaincinv($\cdot,\cdot$)) in scipy, can be
used to compute the _radius_ of an $N_{\mathrm{var}}$-ball $B_{r}$ that
contains a point with prescribed probability $p_{\mathrm{in}}$:

$r(p_{\mathrm{in}};{N_{\mathrm{var}}})\equiv
F_{\rho}^{-1}(p_{\mathrm{in}};{N_{\mathrm{var}}})=\sqrt{2\,\left[P^{-1}({N_{\mathrm{var}}}/2,p_{\mathrm{in}})\right]}.$
(16)

###  3.2 Exterior of an ${N_{\mathrm{var}}}$-ball in standard Gaussian space

Analogously, one can also define the exterior of $B_{\rho}$ schematically
represented by the red domain in Fig. 1 top middle. Therefore, the probability
corresponding to the region outside the ball, i.e., all points with the radial
distance $\rho\in\left(r,\infty\right)$, must be complementary to the function
$F_{\rho}(\rho;{N_{\mathrm{var}}})$

$\displaystyle p_{\mathrm{out}}(r;{N_{\mathrm{var}}})$ $\displaystyle\equiv
1-p_{\mathrm{in}}(r;{N_{\mathrm{var}}})=1-F_{\rho}(r;{N_{\mathrm{var}}})$
$\displaystyle=1-P\left({N_{\mathrm{var}}}/2,r^{2}/2\right)$ (17)
$\displaystyle=Q\left({N_{\mathrm{var}}}/2,r^{2}/2\right),$

where $Q\left(\cdot,\cdot\right)$ is known as the _regularized upper
incomplete gamma function_. This “special function” is again a standard part
of various mathematical libraries such as scipy available in Python (function
gammaincc($\cdot,\cdot$)). Analogously to Eq. (16), the radius of the ball
_not_ containing a point with given probability $p_{\mathrm{out}}$ can be
obtained by the inversion of Eq. (3.2)

$r(p_{\mathrm{out}};{N_{\mathrm{var}}})=\sqrt{2\,\left[Q^{-1}({N_{\mathrm{var}}}/2,p_{\mathrm{out}})\right]}.$
(18)

Of course, the radii $r$ from Eqs. (16) and (18) match when
$p_{\mathrm{in}}+p_{\mathrm{out}}=1$.

The exterior of an $N_{\mathrm{var}}$-ball is an important geometrical entity
for reliability analysis. Consider a situation in which one samples from the
exterior via importance sampling. The _sampling density_ is nonzero only
outside the $N_{\mathrm{var}}$-ball and is obtained by re-scaling the standard
Gaussian density limited to
$\lVert{\boldsymbol{x}}\rVert\in({\color[rgb]{1,0,0}r},\infty)$ so that the
probability content outside the ball is still unit

$\displaystyle h_{\mathrm{out}}\left({\boldsymbol{x}};r\right)$
$\displaystyle=\displaystyle\frac{\Pi_{v=1}^{{N_{\mathrm{var}}}}\varphi({\boldsymbol{x}}_{v})}{p_{\mathrm{out}}(r,{N_{\mathrm{var}}})}.$
(19)

In practice, it is more convenient to obtain this sampling density for point
${\boldsymbol{x}}$ in terms of its radial distance $\rho$ from the origin.
Such left-bounded $\xi$ density is just the re-scaled density from Eq. (9)

$\displaystyle h_{\mathrm{out}}\left(\rho;r\right)$
$\displaystyle=\displaystyle\frac{f_{\rho}(\rho;{N_{\mathrm{var}}})}{p_{\mathrm{out}}(r,{N_{\mathrm{var}}})}=\displaystyle\frac{f_{\rho}(\rho;{N_{\mathrm{var}}})}{1-F_{\rho}(r;{N_{\mathrm{var}}})}.$
(20)

Indeed, in order to randomly sample a point from the exterior of an
${N_{\mathrm{var}}}$ ball $B_{r}$, we exploit the rotational symmetry of the
density in standard Gaussian space and decompose the task into two sub-tasks:
(i) sampling a random direction with the uniform distribution (which is an
$({N_{\mathrm{var}}}-1)$ dimensional problem) and (ii) sampling a random
Euclidean distance from the origin, which is simply a one-dimensional problem.
A practical method of performing the first step, i.e., sampling a random unit
direction in $N_{\mathrm{var}}$ dimensions, is covered in A, and we propose a
method to improve the spread of multiple random directions in B. What then
remains is to choose a univariate random Euclidean distance from the origin,
which is a one-dimensional problem. Component-wise multiplication of the
previously obtained ${N_{\mathrm{var}}}$ unit direction vectors
$\boldsymbol{s}$ by random distances that obey a left-bounded $\chi$ variable
will move the points along their radii-vectors to the desired distance from
the origin $d\in\left({\color[rgb]{1,0,0}r},\infty\right)$. This random
distance corresponding to a sampling probability $p\in(0,1)$ can be obtained
via inverse transformation as

$d_{\mathrm{out}}({\color[rgb]{1,0,0}r})=F_{\rho}^{-1}\left[p+(1-p)\,F_{\rho}(r;{N_{\mathrm{var}}})\right].$
(21)

Finally, the corresponding point can be obtained simply as

${\boldsymbol{x}}_{\mathrm{out}}=d_{\mathrm{out}}\;\boldsymbol{s}.$ (22)

If a sample of $n$ points all at once is needed, it is advisable to cover the
distances from the origin uniformly with respect to probability, and therefore
a set of even space sampling probabilities
$\boldsymbol{p}:\\{p_{i}=(i-0.5)/n,\;i=1,\ldots,n\\}$ is recommended.

###  3.3 Interior between two ${N_{\mathrm{var}}}$-balls in standard Gaussian
space

A rotationally symmetrical region between two different radii, $r<R$, forms an
annulus (a ring), $B_{\mathrm{ann}}$; see the red domain in Fig. 1 top right.
Its probability density for a point ${\boldsymbol{x}}$ is again just a scaled
standard Gaussian density

$\displaystyle h_{\mathrm{ann}}\left({\boldsymbol{x}};r,R\right)$
$\displaystyle=\displaystyle\frac{\Pi_{v=1}^{{N_{\mathrm{var}}}}\varphi({\boldsymbol{x}}_{v})}{p_{\mathrm{ann}}}=\displaystyle\frac{\Pi_{v=1}^{{N_{\mathrm{var}}}}\varphi({\boldsymbol{x}}_{v})}{p_{R}-p_{r}}$
(23)

and zero otherwise. In this equation, the scaling denominator
$p_{\mathrm{ann}}=p_{R}-p_{r}$ is formed by the difference between two
probabilities, $p_{R}=p_{\mathrm{in}}(R,{N_{\mathrm{var}}})$ and
$r_{r}=p_{\mathrm{in}}(r,{N_{\mathrm{var}}})$. Again, this density can be
rewritten in terms if the radial distance $\rho\equiv\lVert x\rVert$. Whenever
$r<\rho<R$, the probability density becomes nonzero, and it reads

$\displaystyle h_{\mathrm{ann}}\left(\rho;r,R\right)$
$\displaystyle=\displaystyle\frac{f_{\rho}(\rho;{N_{\mathrm{var}}})}{p_{\mathrm{ann}}}.$
(24)

The cumulative distribution for this distance reads
$p(\rho)=\int_{r}^{\rho}f_{\rho}(t;{N_{\mathrm{var}}})/(p_{R}-p_{r})\mathrm{\,d}\hskip
0.50003ptt=\left(F_{\rho}(\rho;{N_{\mathrm{var}}})-p_{r}\right)/(p_{R}-p_{r})$
and its inversion can be used to obtain distances corresponding to selected
sampling probability $p\in(0,1)$

$d_{\mathrm{ann}}=F_{\rho}^{-1}\big{(}p\,(p_{R}-p_{r})+p_{r};{N_{\mathrm{var}}}\big{)}.$
(25)

Sampling random points from the annulus can be performed in the same fashion
as described above, i.e., by using a set of random unit directions from Eq.
(67) and scaling them via Eq. (22) by using the distance $d_{\mathrm{ann}}$
from Eq. (25).

A natural question arises: why would we sample points from the annulus between
two balls? In the proposed algorithm, the Gaussian space can conveniently be
divided into the interior of a ball of radius $r$, and its exterior. If we are
no longer interested in exploring the ball interior (e.g., we are sure no
failure occurs below distance $r$), and we believe that the most central
failure point is roughly at a distance $r$ from the origin, it makes no sense
to consider very remote locations to estimate the failure probability by
integrating the Gaussian density (the exterior of the $N_{\mathrm{var}}$-ball
extends to infinity). It rather makes sense to focus only on regions that
contribute considerably to the exterior probability, i.e., to limit the outer
radius to a value $R$. In this way, the described annulus, or _important ring_
, is formed.

Assume now the worst possible scenario in which the failure region is exactly
the exterior of the ball with radius $r$. In this case we know the
corresponding failure probability $(1-p_{r})$ via knowing the probability
inside the ball $p_{r}=F_{\rho}^{-1}(r;{N_{\mathrm{var}}})$ (safe state). Note
that such a failure domain represents an upper bound of failure probability
because any other failure surface will have a smaller failure probability. It
then makes sense to exclude regions outside a greater ball with radius $R$ by
which only a negligible fraction of $1-p_{r}$ is excluded. For example,
consider excluding the probability of $(1-p_{r})/10^{4}$ only. The outer
diameter $R$ then obeys $1-p_{R}=(1-p_{r})/10^{4}$ meaning that
$p_{R}=1-(1-p_{r})/10^{4}$. Therefore, the outer radius can be computed via
Eq. (16) as

$R=F_{\rho}^{-1}\big{(}1-(1-p_{r})/10^{4};{N_{\mathrm{var}}}\big{)}.$ (26)

##  4 Sequential adaptive _extension_ of the experimental design

The premise behind the development of the presented algorithm is that each
evaluation of the system performance (limit state function) is very expensive.
Crude Monte Carlo type sampling strategies simply throw sampling points ($=$
integration “nodes”) independently of each other, and do not reflect the
structure of the problem being solved at all. Some advanced sampling
techniques learn from previously obtained information but may not be able to
fully exploit its potential.

The set of points at which the true performance function has been evaluated so
far will be referred to as the _experimental design_ (ED), and their number
will be denoted as $N_{\mathrm{sim}}$. The ED is simply a table of
$N_{\mathrm{sim}}$ points each with $N_{\mathrm{var}}$ coordinates and one
additional vector of ${N_{\mathrm{sim}}}$ results obtained from the
performance function. These results can be indices of one of the _discrete
states_ , including failure, success, no result, or a failure type code. The
algorithms proposed in this paper are formulated to operate with such limited
information. In particular, no other information is necessary for either of
the two fundamental steps, namely (i) _extension_ of ED via the selection of
the best candidate point for subsequent evaluation of the expensive
performance function, and (ii) _estimation_ of the current failure probability
using a temporary binary surrogate model (or a surrogate with discrete states
as listed above). These two steps are recognized in Fig. 2, which presents an
overview of the proposed method. It involves the construction of two different
surrogate models: one for the identification of domains for local refinement,
and the other for fast numerical integration via sampling.

Figure 2: Flowchart of the proposed methodology consisting of two basic tasks:
extension and estimation.

The existing methods for the extension of the ED often employ learning
functions [35]. Learning functions are convenient mathematical functions that
weigh the metamodel properties to seek the best candidate to extend the
current ED. They evaluate a set of candidates with criteria that are built on
considerations of uncertainty in the model approximation and the proximity to
the failure region and select the new most promising candidate. Bichon et al.
[74] and Echard et al. [48] introduced two of the most relevant works in this
context: Bichon et al. [74] presented “Efficient Global Reliability Analysis”,
proposing the use of the “Expected Feasibility Function” to extend the ED, and
Echard et al. [48] developed a method combining Kriging with Monte Carlo
Simulation while using what they called a U-function based on the probability
of misclassifying a candidate in order to extend the ED. The literature on
learning functions is quite rich. Lv et al. [75] proposed the H learning
function, which is built on entropy consideration, and combined Kriging with
Line Sampling, Sun et al. [76] defined a learning function named the Least
improvement function (LIF), which combines Kriging-based statistical
information and the joint PDF of basic random variables. Zhou and Lu [77]
proposed a learning function in Polynomial Chaos Expansion that models
uncertainty with a Bayesian approach. Quite recently, Zhang et al. [78]
proposed the Reliability Expected Improvement Function (REIF), which is to be
combined with an adaptive Kriging surrogate, which relates to the expected
improvement (EI) of Jones et al. [79]. Indeed, when the performance function
is smooth and not just binary, its value is deemed to provide an additional
measure of the distance from the failure surface. In such cases, this
information can be used to increase the effectiveness of the algorithm in
various ways. One of the possibilities is to construct a surrogate model based
on the existing ED. The extension of the ED can balance between the
exploration and exploitation of sampling using conditional probability [80].
All of these approaches use the numerical values of the true performance
function $g(\boldsymbol{x})$ and its surrogate.

In the present paper, we consider that the only pieces of information to be
supplied for the selection are: (i) the position of the points in ED, (ii) the
classification of the model output (failure, failure type, success, no
information, etc.) and (iii) the joint probability density function. In order
to discover new failure domains and, at the same time, to refine the existing
approximation of the failure surface (boundary), we prepare two sets of
candidate points: the _exploration set_ (Sec. 4.1) and the _exploitation set_
(Sec. 4.2), combine them together and perform the selection of the best
candidate using the proposed $\psi$ criterion (Sec. 4.3).

###  4.1 Exploration set (global)

We expect that a rare event (failure) region can appear anywhere. The local
probability density of any potential failure region size decreases with
increasing distance from the origin in the standard Gaussian space. Therefore,
in a global search for failure domains, it makes sense to explore the space in
a controlled fashion by checking various directions with an increasing
distance from the origin. Therefore, we propose the construction of a sequence
of nested $N_{\mathrm{var}}$-dimensional balls in the Gaussian space whose
radii increase towards infinity and which occupy the probability in prescribed
levels. In other words, we suggest preparing an _exploration set_ for “onion-
like” layers; see Fig. 3 left. Each layer is numbered by a level number $i$,
and the corresponding $N_{\mathrm{var}}$-ball is to be covered by a prescribed
number of points $n_{i}$ corresponding to this layer. The boundary of each
layer is an $N_{\mathrm{var}}$-ball, and the $n_{i}$ points selected from the
surface of this ball should be uniformly dispersed. Moreover, it is preferable
to avoid directionally collapsible candidates: the candidates corresponding to
level $i$ should have different radii vectors than the points from the
preceding level $i-1$.

Figure 3: The first seven levels of exploration sets (layers) in
${N_{\mathrm{var}}}=2$ and ${N_{\mathrm{var}}}=3$ dimensional standard
Gaussian space. Panels a,c: nested $N_{\mathrm{var}}$-balls. Panels b,d:
realizations of random exploration sets.

The probability content occupied by an $N_{\mathrm{var}}$-ball corresponding
to layer $i$ is selected with regard to what level of failure probability the
ball is roughly associated with. More exactly, the radii of the balls can be
selected using Eq. (16) such that the exterior of that ball has a prescribed
probability, $p_{\mathrm{out}}$, see Eq. (3.2).

There is freedom in the proposed method to select: (i) probabilities
$p_{\mathrm{out},i}$ corresponding to individual levels and (ii) the numbers
of points $n_{i}$ corresponding to these individual levels. The selection of
these parameters is an important step. If the exploration algorithm is too
_aggressive_ , meaning that the ball radii are rapidly increasing while they
are each covered by only a small number of points, there is a risk of
overlooking some localized failure domain $\pazocal{F}$ (regions which are
narrow at small radii; see the illustration in Fig. 3b in which the first four
levels of exploration points miss the tip of the failure domain). In such a
case, it takes several other “onion layers” and expansive $g(\boldsymbol{x})$
evaluations to hit the failure domain for the first time. As will be shown
below, once a distant failure is hit, the algorithm will automatically back-
trace the failure surface down the high probability region; however, many
unnecessary calls to the performance function will be made. If, on the other
hand, the “onion layers” are too densely packed and/or each covered by too
many exploration points, it will take many calls to hit the failure domain for
the first time and start the exploitation phase. In the simplest situation,
the exploration set is pre-generated in advance for a given problem dimension.
We remark, however, that the exploration set can be modified anytime during
the run of the proposed method; see below. For example, suppose the
exploration phase has already reached a rather high radius, and the analyst
suspects that a narrow convex failure domain might have been overlooked. The
exploration set can be enriched with additional layers or candidates at any
time, and the selection criterion described below will automatically select
the most relevant candidate(s) to check for potential failures.

Once the decision regarding the directions/radii coverage is made, it is not
difficult to select $n_{i}$ points for each level approximately uniformly
distributed on the unit sphere and scale them to the desired radius
$\rho_{i}$. We refer to A and B for a simple algorithm for the selection of
the set of $n_{i}$ points from the surface of $N_{\mathrm{var}}$-dimensional
ball.

Based on our experience, it is a good practice to select the
$N_{\mathrm{var}}$-balls so that each level leaves an unexplored exterior with
probability $p_{\mathrm{out},i}=10^{-i}$, where $i=\\{1,2,3,\ldots\\}$. This
means that the first ball consumes 90% and leaves only a probability content
of 0.1, the second ball only leaves 0.01, the third 0.001, etc. Given these
levels, the surfaces of these balls are sufficiently covered by exploration
points when their number is selected as

$n_{i}=-{N_{\mathrm{var}}}\ln\left(\frac{p_{\mathrm{out},i}}{{N_{\mathrm{var}}}}\right).$
(27)

These point counts, e.g., rounded down to integers using the floor function
$n_{i}=\left\lfloor n_{i}\right\rfloor$, will be used in all the numerical
examples presented in this paper. With this rule, the number of evaluations to
cover a certain safety level is known in advance, and so is the number of
function calls for a given failure probability. Of course, any prior knowledge
about the failure probability can be used to adjust the radii. To get an idea
about the point counts proposed in Eq. (27) for dimensions ranging between two
and twenty, see Fig. 23 and Tab. 2 in C.

Until a rare event is encountered in the ED for the first time (be it a
failure or a non-result), the candidate set for selection via the $\psi$
criterion is formed by this _exploration set_ only. Once the desired event is
hit, the surroundings of the point must also be considered in a search for the
most informative candidate. In other words, an _exploitation set_ with
appropriate candidates must be constructed and included when selecting the
next ED point.

###  4.2 Exploitation set (local)

Assume now that the previously evaluated performance function confirmed at
least one rare event (e.g., failure) at any point of the ED. It is clear that
the failure surface has been crossed, and the location of the boundary should
be refined locally. In the case of a binary result, it is reasonable to test a
location somewhere close to existing failure and success points. More exactly,
we propose searching for the boundary at locations (candidates) whose _two
nearest neighbors have different classifications_ (failure, success, no result
or simply _different states_). From all candidates, we propose that only those
having two _different_ classifications of their two nearest neighbors be
retained. The exploitation set is prepared sequentially, i.e., after each
extension of the ED, and so it is based on current information. The
preparation of the candidate set works in two steps: (i) the generation of a
large pool of “dots” in a sufficiently large neighborhood of each existing
rare event point (e.g., a failure point) and (ii) the selection of only those
candidates that potentially refine the boundary between different states (the
two nearest ED neighbors signalling a different state). The censored list of
points is the _exploitation set_.

The pool of candidates can be pre-generated at the beginning and kept constant
throughout the whole analysis. We, however, propose generating them after each
discovery of a new rare event (e.g., evaluation of the limit state function
returning failure). A set of “candidate dots” is simply generated in the
spirit of importance sampling with Gaussian sampling density
$h_{i}({\boldsymbol{x}})$ centered at each rare-event ED point, and the
standard deviation is driven by the problem dimension:

$\displaystyle\mu_{i}$ $\displaystyle={\boldsymbol{x}}_{i},\quad
i=1,\ldots,{N_{\pazocal{T}}}$ (28) $\displaystyle\sigma_{i}$
$\displaystyle=\sqrt{{N_{\mathrm{var}}}-1},$

which spreads the dots sufficiently away from the existing point but at the
same time focuses the dots sufficiently. Since the censoring that extracts
(masks) the candidates for the exploitation set is fast (based merely on the
computation of distances), it is quick to generate and censor $10^{6}$
candidates or more. Moreover, we remark that many operations can be pre-
computed for the potential classifications (failure, success, etc.) while the
true performance function is still being evaluated. Once the performance
function $g(\boldsymbol{x})$ has returned the result, one of the pre-computed
scenarios is used to call $g(\boldsymbol{x})$ with the best candidate
immediately.

Figure 4: Two-dimensional illustration of the combined set candidates and the
selection process (standard Gaussian space) for a binary performance function
dividing the space into a “flower-shaped” safe domain and the surrounding
failure domain. At ${N_{\mathrm{sim}}}=9$ model evaluations, the rare event is
hit for the first time. At ${N_{\mathrm{sim}}}=10$ an _exploitation_ candidate
is selected, while at ${N_{\mathrm{sim}}}=11$ the _exploration_ candidate wins
using the highest value of the $\phi$ criterion expressed by the blue color
saturation. A complete evolution of the extension of the ED using the
$\phi$-based selection process for another run of the same problem is
available as a video.

Fig. 4 presents an illustration of the exploitation sets. The figure shows the
situation at four different stages after ${N_{\mathrm{sim}}}$ evaluations of a
binary performance function. The true failure surface is plotted via the blue
curve (a complicated shape with seven “design points”). The algorithm started
with the exploration candidates (empty circles), and the results obtained from
evaluations of the limit state function are visualized by _filling_ the
circles (either green for “success” or red for “failure” event). Once the
first red point was discovered, the algorithm started proposing _exploitation_
candidates (the results obtained in  these exploitation candidates can be
distinguished from the exploratory ones by the fact that the colored circles
do not have a black ring around them). The boundaries of this region are, in
fact, axes of lines connecting ED points with different outcomes. These lines
are plotted via solid lines of various colors. Note, however, that neither
these lines nor the corresponding boundaries are constructed. The boundaries
are naturally formed by censoring candidates with two different nearest ED
neighbors. The boundaries in higher dimensions become formed by planes in 3D
and hyperplanes in higher dimensions. No computational geometry is needed, and
generalization to high dimensions is straightforward.

Now that the combined exploitation and exploration set has been prepared, the
last step is the selection of the best candidate. In the situation depicted in
the figure, the best candidate (blue cross) is selected from the exploration
set, but this selection depends on the local situation and balances between
the local and global sets to maximize the potential gain from the next
evaluation of the performance function. The selection criterion is described
in the following section.

###  4.3 $\psi$ criterion for candidate selection

Suppose the combined exploration-exploitation set of candidates is available.
The task is to select the candidate in which the (supposedly costly)
evaluation of the performance function delivers the maximum gain in terms of
probability. The new point will support the geometrical interpretation of the
failure domain which is being constructed. The structural information about
the shape and location of the failure surface is important, e.g., for
sensitivity analysis [81], see also Sec. 6, and may also be useful for other
approximation techniques (such as the Taylor expansion at “design points”
etc.).

One of the main motivations behind developing the proposed algorithm is that
no two ED points may be required close to each other unless they deliver a
significant amount of information in terms of probability content (that is,
refine the boundary approximation between territories corresponding to
different event types – such as the failure surface $\partial{\pazocal{F}}$).

We propose a novel $\psi$ criterion for the selection of the best candidate.
Its maximization leads to the (approximate) maximization of the instantaneous
gain in terms of probability content because the meaning of the criterion is
the approximate amount of probability content being classified by evaluating
$g(\boldsymbol{x})$ in the candidate. Indeed, the values of $\psi$ can be
viewed as “bites” of probability with a clear geometrical meaning: each
candidate represents a certain region (its neighborhood) in the design space
and the volume of this region can be multiplied by the average probability
density to obtain the corresponding probability content:

$\psi_{\mathrm{c}}=\underbrace{\sqrt{f_{\mathrm{c}}f_{\mathrm{s}}}}_{\mathrm{ave\;probability}}\,\underbrace{\left(l_{\mathrm{c},\mathrm{s}}\right)^{{N_{\mathrm{var}}}}}_{\propto\mathrm{vol.}},$
(29)

where the term $\sqrt{f_{\mathrm{c}}f_{\mathrm{s}}}$ accounts for the
(Gaussian) probability density of both the candidate, $\mathrm{c}$, and the
density in its nearest ED point, $\mathrm{s}$, (i.e., the existing =
previously evaluated point in the ED). The term represents the geometrical
average probability density of the two. The geometric mean between $n$ numbers
$x_{i}$ is defined as $(\Pi_{i=1}^{n}x_{i})^{1/n}$, and therefore the proposed
term features $\sqrt{f_{\mathrm{c}}f_{\mathrm{s}}}$. Why do we use the
geometric mean instead of the simple arithmetic mean? The arithmetic mean is
not a suitable measure as it would favor infinitely distant candidates. The
problem is that very distant candidates in Gaussian space that have zero
probability density would still form a nonzero average density with an
existing ED point. Therefore, very remote candidates would have a high $\psi$
criterion and become selected, which is not preferable – the expansion to
distant locations is already controlled by the exploration set.

Figure 5: Geometrical meaning of the proposed $\psi$ criterion for a candidate
“$\mathrm{c}$” in ${N_{\mathrm{var}}}=1$ and ${N_{\mathrm{var}}}=2$
dimensions. The black solid circles are the existing ED points and the blue
circle “$\mathrm{s}$” is the nearest neighbor to candidate “$\mathrm{c}$”.
$\psi_{c}$ is therefore proportional to the probability content of the red
domain close to the candidate.

The term $l_{\mathrm{c},\mathrm{s}}$ is the distance between a candidate,
$\mathrm{c}$, and its nearest existing ED point, $\mathrm{s}$. Simply, the
distances between all candidates and all $N_{\mathrm{sim}}$ currently existing
ED points are evaluated, and each candidate is paired with its nearest point.
To compute the $N_{\mathrm{sim}}$ distances from each candidate, a suitable
metric must be selected. It this paper, we consider Euclidean distance between
the candidate and its nearest neighbor as

$l_{\mathrm{c},\mathrm{s}}=\sqrt{\sum_{v=1}^{{N_{\mathrm{var}}}}\big{|}{\boldsymbol{x}}_{v}^{(\mathrm{c})}-{\boldsymbol{x}}_{v}^{(\mathrm{s})}\big{|}^{2}}.$
(30)

When this distance is raised to the problem dimension:
$\left(l_{\mathrm{c},\mathrm{s}}\right)^{{N_{\mathrm{var}}}}$, it becomes
proportional to the volume of a geometrical object selected from the design
domain. This geometrical object, be it an $N_{\mathrm{var}}$-ball,
$N_{\mathrm{var}}$-hypercube, or some other shape, has a volume that
quantifies the size extent of the neighborhood of the candidate $c$. We do not
need to specify the exact shape of the volume as all the geometrically similar
volumes for various candidates under comparison would share the same positive
multiplier related to the specific geometric shape. This multiplier can be
dropped from the definition as it does not alter the ranking of compared
candidates. The geometrical interpretation of the criterion is illustrated in
Fig. 5.

The criterion is constructed as a product of two independent terms:

  * •

the volume term, which favors rapid _exploration_ of previously unexplored
regions by expansion, and

  * •

the probability term, which favors the _exploitation_ of high probability
contents occupied by a candidate near the origin.

The balance between exploration and exploitation is thus maintained
automatically: a candidate located close to an existing ED point can only be
selected if the corresponding probability density is sufficiently high so that
the gain is greater compared to other more distant points. And vice-versa,
overly distant candidates are not preferred because their densities are low.
The criterion favors candidates leading to a fast classification of the space
because, the $\psi$ criterion at a given stage of the computation is the
maximum value of all $\psi_{\mathrm{c}}$ selected among all candidates under
comparison

$\psi=\max_{\mathrm{c}}\left(\psi_{\mathrm{c}}\right).$ (31)

The illustrative example in Fig. 4 reveals that the criterion automatically
switches between the exploration and exploitation sets. Points selected from
the predefined exploration set have a black circumference and are filled with
either red or green color. Filled circles without the black circumference have
been selected from the exploitation sets, and these points could have appeared
only after the discovery of the first rare event point (red). All candidates
having two different types of two nearest neighbors are considered to be
exploitation candidates and ranked together with the exploration candidates
using the proposed $\psi$ criterion. The maximum $\psi$ criterion values
always correspond to the blue crosses. During the exploitation steps, the
algorithm tends to alternate between safe and failure domains by refining the
approximation of the failure domain.

It is useful to record the history of the $\psi$ criterion throughout the
sequential extension of the ED. The criterion provides a rough orientation as
to how much the next candidate can modify the current estimation of
$p_{\pazocal{F}}$, and thus these two pieces of information may serve for
formulation of a  _stopping criterion_. If the value of $\psi$ drops below a
certain limit (say ${p_{\pazocal{F}}}/100$), it might be reasonable to quit
the algorithm and/or switch to another type of analysis. Experience has shown
that $\psi$ will keep decreasing if the computation continues, and typically
the rate of decrease slows down as further refinements of the boundaries cause
only minor changes to the estimations.

This formulation of the stopping criterion does not guarantee that the
algorithm has already discovered all important regions. It can happen that
further computation will lead to the discovery of new, previously overlooked
regions. In such cases, the estimation of underestimated ${p_{\pazocal{F}}}$
may suddenly increase. Analogously, it may also happen that the
${p_{\pazocal{F}}}$ estimation is exaggerated, especially in high dimensions.
If the ED is very small for a problem dimension, it is possible that the
discovery of a failure point is not followed by refinement of the boundary,
and an overly rough approximation of the true failure surface via an
oversimplified surrogate model may then estimate the extent of the failure
domain as being too large, leading to the overestimation of the true
${p_{\pazocal{F}}}$.

###  4.4 Discussion on the $\psi$ criterion and the assumptions behind the
proposed algorithm

#### High dimensions and the hyperspherical method

As pointed out by Katafygiotis and Zuev [82], the geometry of the reliability
problem in high dimensions is challenging: the volume of the probability space
grows with the number of random variables. In such a large volume, most of the
contributions to the failure probability come from failure regions with small
probability densities but large volumes. Katafygiotis and Zuev [82] discuss
the significance of the “design point” (the most central failure point) in the
case of strongly nonlinear problems, and they show that the selection of an
appropriate importance sampling density is practically impossible in these
cases. Indeed, high-dimensional reliability problems are hard to treat using
importance sampling (no matter if the sampling density is centered in the
origin or multiple densities are used centered in the identified design
points). When the form of the importance sampling density is not appropriate,
the variance of importance sampling estimation is known to explode with domain
dimension [25]. Katafygiotis and Zuev [82] introduced a concept of an
Important Ring, which is a rotationally symmetrical region between radii
${N_{\mathrm{var}}}\pm\epsilon$ centered at the origin of the design space,
and explained why IS with density centered around the “design point” becomes
impractical. These geometric insights were exploited by Wang and Song [83],
who presented a hyper-spherical extrapolation method for high-dimensional
problems. Indeed, the shape of an annulus, Fig. 1 left, and Eq. (12) reveal
the fact that most of the probability is associated with a thin layer around
an $N_{\mathrm{var}}$-ball with a radius of approximately
$\sqrt{{N_{\mathrm{var}}}-1}$. The traditional assumption that the region in
the vicinity of the most central point (a.k.a. the design point) has a
prevalent contribution to the failure probability may not be valid, especially
for very high dimensions ${N_{\mathrm{var}}}$. The present method solves this
problem automatically by employing the proposed $\psi$ criterion, which
automatically favors the largest probability contributions of candidates, no
matter where they are.

#### Assumptions and limits

As Rackwitz [1] pointed out, an important step in the development of methods
is to _show_ where they do _not_ work, i.e., to find the limits of the
applicability of a concept and to construct counterexamples. The assumption
behind the developed algorithm in its simplest version is that the
corresponding set (territory) forms an _open domain_. As has now become clear,
cases with scattered localized _closed_ rare event domains may not be treated
well; it can happen that localized domains will be missed by the exploration
set and encapsulated by the application of the proposed $\psi$ criterion
during the refinement. A separate numerical example featuring the “Modified
Rastrigin function” (see Fig. 9) is devoted to this weakness and is presented
in Sec. 7.6.

#### Incorporation of a priori knowledge

Sometimes, the performance function represents a problem for which the outcome
can be known _a priori_ , i.e., without the need to run an expensive computer
simulation. For example, suppose a specific region is known to be a safe
region. This information can be passed to the algorithm easily without
modifying the proposition of candidates or modifying the selection criterion.
Whenever a candidate is selected to become an ED point and a
$g(\boldsymbol{x})$ call is required, the point location is checked first, and
if the point belongs to the described region, the call is simply bypassed by
associating the ED point with the known classification.

##  5 Failure probability _estimation_ via sampling analysis

At any stage of the ED extension process, the desired probabilities can be
estimated based on point-wise information, i.e., the current ED with known
$g(\boldsymbol{x})$ outcomes; see the right-hand part of Fig. 2. Indeed, a
true computationally expensive model (a simulator) can be substituted by a
computationally cheaper model. The surrogate model is constructed in a solely
non-intrusive way with respect to the original simulator, i.e., it is purely
data-driven. When the performance function $g(\boldsymbol{x})$ returns
continuously distributed, trustworthy, and well-behaved output, traditional
surrogate models such as PCE, Kriging, radial basis functions, etc., can be
employed for fast sampling analysis. Such a surrogate model has the potential
to improve the accuracy of the importance sampling estimation presented in
this section because the failure surface (the boundary
$\partial{\pazocal{F}}$) may be approximated more accurately. It is guaranteed
that the above-described extension of the ED was performed in such a way that
the surrogate model was well supported, especially when close to the failure
surface.

In this paper, however, we focus on cases in which the original simulator is
trusted to provide categorical information only, i.e., one of a finite set of
classifications. Similarly, the surrogate will simply be a finite-state
classifier.

###  5.1 Rare event surrogate via the nearest neighbor classifier

We propose the construction of an intermediate discrete-valued surrogate at
any location ${\boldsymbol{x}}$ based simply on its nearest ED neighbor at
which the performance of the system has already been evaluated. For example,
an integration node ${\boldsymbol{x}}$ will only be considered as a “failure
location” if the nearest neighbor from the ED signaled failure. In this way,
an effective approximation indicator function
$\boldsymbol{1}_{\pazocal{F}}^{{N_{\mathrm{sim}}}}\left({\boldsymbol{x}}\right)$
of the true indicator function
$\boldsymbol{1}_{\pazocal{F}}\left({\boldsymbol{x}}\right)$ is obtained at any
stage with $N_{\mathrm{sim}}$ evaluated points. Analogously, a general event
$\pazocal{T}$ is estimated to occur if the known model response in the nearest
neighbor is $\pazocal{T}$ and a collection of such points becomes associated
with a nonzero indicator value
${\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{X}}\right)}$. This
corresponds to the division of the design domain $\pazocal{D}$ into a finite
number of territory types (see the two top right parts of Fig. 6) with just a
binary classification. The surrogate must decide about the previously
“undecided” regions visible in the top left part of the figure, i.e., region
considered in the ED extension as potential boundary regions between two
different classes of model outputs.

Figure 6: Top row: on-the-fly classification of territories. Top left panels:
two nearest neighbors for the classification of boundary regions – red and
white candidates are censored out, and the retained (blue) candidates are
considered for the _extension_ of ED. Top right panels: classification of dots
and nodes via the nearest neighbor; used for screening and cubature. Bottom:
history of the rare event probability _estimation_ and the $\psi$ criterion
quantifying the approximate “probability bites” by selected candidates. The
complete evolution of all panels is demonstrated in the Wavy circle video.

Making a response prediction based on its nearest neighbor is equivalent to
constructing a Voronoï diagram, i.e., Dirichlet tessellation, which is a
standard way to partition a domain into disjoint cells called Dirichlet
regions ($=$Thiessen polytopes $=$ Voronoï polygons). The $N_{\mathrm{sim}}$
points from the current ED with known classifications play the role of the
_seeds_ (sometimes called sites or generators). This Voronoï-type
approximation may seem to be too crude, but as will be shown below, when the
classification is built using points obtained via the proposed sequential
adaptive algorithm based on the $\psi$ criterion,
$\boldsymbol{1}_{\pazocal{T}}^{({N_{\mathrm{sim}}})}\left({\boldsymbol{x}}\right)$
quickly tends to the true
${\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{x}}\right)}$ as
$N_{\mathrm{sim}}$ increases. An important aspect is that no computational
geometry is involved; the only operation needed is the comparison of
distances.

We remark that the described distance-based surrogate model is not restricted
to the Gaussian space only. The surrogate model and also the estimation of
$p_{\pazocal{F}}$ can be made either in the standard Gaussian space or, it can
be computed directly in the original space, standardized to zero mean and unit
variance. If, however, any of the variables are bounded, then experience has
shown that the proposed sequential refinement and classification becomes less
efficient.

###  5.2 Importance sampling estimation

Assume the surrogate model for event $\pazocal{T}$ is ready in the form of
indicator function
${\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{X}}\right)}$. The related
event probability defined by Eq. (5) can today be estimated very efficiently,
e.g., via importance sampling (IS). In IS the original joint density
$f_{{\boldsymbol{X}}}({\boldsymbol{x}})$ for sample selection is replaced by
the importance sampling density $h({\boldsymbol{x}})$, which we consider to be
jointly Gaussian with independent marginals. Generating millions of
realizations of Gaussian random variables takes less than a second on
contemporary computers, and implementations are provided ready to use in
various free numerical packages (for example the NumPy package [84] for
powerful operations over multidimensional dimensional arrays in combination
with vectorized functions and routines [85] from scipy available in Python).

The probability of an event type $\pazocal{T}$ occurring when the indicator
function ${\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{X}}\right)}$ signals
it is defined as the expectation:
${p_{\pazocal{T}}}=\mathbb{E}\left[{\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{X}}\right)}\right]=\idotsint_{\pazocal{D}}{\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{x}}\right)}\,f_{{\boldsymbol{X}}}({\boldsymbol{x}})\,\mathrm{\,d}\hskip
0.50003pt{\boldsymbol{x}}$. Let $h({\boldsymbol{x}})$ be the IS density which
is positive wherever event $\pazocal{T}$ occurs. The probability of event
$\pazocal{T}$ can be rewritten as

$\displaystyle{p_{\pazocal{T}}}=\int_{\pazocal{D}}\frac{{\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{x}}\right)}\,f_{{\boldsymbol{X}}}({\boldsymbol{x}})}{h({\boldsymbol{x}})}h({\boldsymbol{x}})\mathrm{\,d}\hskip
0.50003pt{\boldsymbol{x}}=\mathbb{E}_{h}\left[\frac{{\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{X}}\right)}\,f_{{\boldsymbol{X}}}({\boldsymbol{X}})}{h({\boldsymbol{X}})}\right],$
(32)

where $\mathbb{E}_{h}\left[\cdot\right]$ denotes the expectation for
${\boldsymbol{X}}$ being distributed according to $h$: ${\boldsymbol{X}}\sim
h$. The IS estimation of ${p_{\pazocal{T}}}$ based on the current
approximation of
${\boldsymbol{1}_{\pazocal{T}}^{({N_{\mathrm{sim}}})}\left({\boldsymbol{x}}\right)}{}$
of the true indicator function
${\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{x}}\right)}$ is made with
$n_{\textsf{IS}}$ integration nodes via the arithmetic average

$\displaystyle{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}$
$\displaystyle\approx\frac{1}{n_{\textsf{IS}}}\sum_{i=1}^{n_{\textsf{IS}}}{\boldsymbol{1}_{\pazocal{T}}^{({N_{\mathrm{sim}}})}\left({\boldsymbol{x}}_{i}\right)}{}\frac{f_{{\boldsymbol{X}}}({\boldsymbol{x}}_{i})}{h({\boldsymbol{x}}_{i})},\quad{\boldsymbol{X}}_{i}\sim
h.$ (33)

The IS estimator is unbiased by construction. The variance of the IS estimator
follows the standard definition

$\displaystyle\text{Var}_{h}[{p_{\pazocal{T}}}]$
$\displaystyle=\mathbb{E}_{h}\left[\left(\frac{{\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{x}}\right)}\,f_{{\boldsymbol{X}}}({\boldsymbol{x}})}{h({\boldsymbol{x}})}\right)^{2}\right]-{p_{\pazocal{T}}}^{2}=\int_{\pazocal{D}}\left(\frac{{\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{x}}\right)}\,f_{{\boldsymbol{X}}}({\boldsymbol{x}})}{h({\boldsymbol{x}})}\right)^{2}h({\boldsymbol{x}})\mathrm{\,d}\hskip
0.50003pt{\boldsymbol{x}}-{p_{\pazocal{T}}}^{2}$ (34)

Its estimation using $n_{\textsf{IS}}$ integration nodes
${\boldsymbol{X}}_{i}$ sampled from $h$, $i=1,\ldots,n_{\textsf{IS}}$, can
again be performed using an arithmetic average

$\displaystyle\text{Var}_{h}[{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}]$
$\displaystyle\approx\frac{1}{n_{\textsf{IS}}}\Bigg{\\{}\left[\frac{1}{n_{\textsf{IS}}}\sum_{i=1}^{n_{\textsf{IS}}}{\boldsymbol{1}_{\pazocal{T}}^{({N_{\mathrm{sim}}})}\left({\boldsymbol{x}}_{i}\right)}{}\frac{f^{2}_{{\boldsymbol{X}}}({\boldsymbol{x}}_{i})}{h^{2}({\boldsymbol{x}}_{i})}\right]-\left({p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}\right)^{2}\Bigg{\\}}$
(35)

and finally, the coefficient of variation of the estimation reads

$\displaystyle\text{CoV}_{h}[{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}]$
$\displaystyle\approx\frac{\sqrt{\text{Var}[{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}]}}{{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}}$
(36)

A straightforward application of IS estimation would be to use Gaussian
sampling density centered at the origin (zero mean values) and a magnified
standard deviation which is somehow based on the smallest distance between the
failed ED point and the center point. Though this is a robust strategy, it is
not very efficient, as many integration nodes fall outside the territories
corresponding to the desired event type. Another disadvantage is the need to
compute the weight ratio for each integration node.

It is known that the optimal choice of the importance sampling density for
event $\pazocal{T}$ (e.g. a failure) is proportional to the original density,
but defined only over the importance (failure) region in which the event
$\pazocal{T}$ occurs:
$h_{\text{opt}}({\boldsymbol{x}})={\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{x}}\right)}\,f_{{\boldsymbol{X}}}({\boldsymbol{x}})/{p_{\pazocal{T}}}$.
Such an optimal density leads to zero variance of the IS estimator Eq. (33).
Generally, using such an optimal density is not feasible for two reasons: the
probability ${p_{\pazocal{T}}}$ being estimated is not known, and even if the
density $h_{\text{opt}}({\boldsymbol{x}})$ can be constructed, a method to
efficiently sample from such a density is generally not available.

We argue that the crucial component in formulating the best IS density is to
determine territories for exclusion from the sampling because the sampling of
points from the important region is then performed proportionally to the
original density. Therefore, we argue that it is worthwhile to obtain the
precious description of the important region, i.e., the territory indicating
the target event $\pazocal{T}$ via
${\boldsymbol{1}_{\pazocal{T}}\left({\boldsymbol{x}}\right)}$. As described in
the preceding subsection, its point-wise approximation can be obtained as
${\boldsymbol{1}_{\pazocal{T}}^{({N_{\mathrm{sim}}})}\left({\boldsymbol{x}}\right)}$.
Therefore, we propose starting with a simple importance sampling run to obtain
“screening dots” providing good geometrical information about the territory of
event $\pazocal{T}$.

#### Screening via a set of localized importance sampling runs

Given the point-wise information from ED, we expect that the target territory
is in the vicinity of points associated with $\pazocal{T}$. From all the
$N_{\mathrm{sim}}$ points obtained so far, we select only those
$n_{\pazocal{T}}$ points which were classified as $\pazocal{T}$-events. In
order to map their surroundings and describe the corresponding Voronoï cells,
we propose the performance of such a set of $N_{\pazocal{T}}$ localized IS
runs around each such ED point that the Gaussian sampling densities $h_{i}$
are always centered at one of the $\pazocal{T}$-classified points, and
$\sigma_{i}=\sqrt{{N_{\mathrm{var}}}-1}$. These are the same kind of dots as
were used as candidates in the _extension_ task, see Eq. (28). It suffices to
throw roughly hundreds or thousands of screening dots around each of the
${n_{\pazocal{T}}}$ points if one wishes to sample the territory fairly
accurately. Based on the nearest neighbor classification, the territory
belonging to $\pazocal{T}$ will now be described by a high number of retained
screening dots and not only the $n_{\pazocal{T}}$ points. Even though the
Voronoï cells are not constructed explicitly via computational geometry (such
as the QuickHull algorithm), their descriptors, such as the locations of
vertices or volumes, can be estimated very accurately [86].

#### Cubature via a set of _local_ IS runs

One can use the retained screening dots as cubature nodes and estimate the
probability of event $\pazocal{T}$ directly in the fashion of IS. Each such
point has its likelihood ratio computed as the fraction
$f_{X}({\boldsymbol{x}})/h_{i}({\boldsymbol{x}})$, where
$h_{i}({\boldsymbol{x}})$ can be different for dots obtained from sampling
around a different $\pazocal{T}$ point. This direct approach is referred to as
“localized IS” from here on and will be reported in the convergence diagrams;
see the magenta points and line in Fig. 6. Such a systematic cubature can be
quite accurate (provided the most important failure regions are well covered
by ED). The disadvantage is the need to compute the weight ratio for each
integration node (the ratio between the original density and the local
sampling density) and also the number of nodes that are not classified as
$\pazocal{T}$-event tends to be high, which degrades the results.

#### Cubature via a  _global_ IS via rotationally invariant density in an
important ring

Suppose the screening dots were not localized in a very small region for which
the localized importance sampling would be efficient. We propose the use of
use a large pool of “integration nodes” selected from sampling density
$h_{\mathrm{ann}}({\boldsymbol{x}})$ which excludes the useless
$N_{\mathrm{var}}$-ball of the radius $r$ ($=$ the distance of the most
central “screening dot” from the origin; see the magenta points in Fig. 6).
Additionally, we propose that the density also excludes the exterior of the
$N_{\mathrm{var}}$-ball with the radius $R>r$. The outer radius can initially
be calculated using Eq. (26), which excludes only a negligible fraction of the
contents of the exterior of the ${N_{\mathrm{var}}}$-ball of radius $r$; see
also Fig. 1 top. Once a previous estimation of the rare event probability
exists, such as the preceding
${p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}-1\right)}}$, the outer radius
should be based on this information, i.e., the excluded probability should
form a negligible portion of it, say
${p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}-1\right)}}/10^{4}$. Using Eq.
(16), the outer radius is obtained as

$R({p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}{\color[rgb]{1,0,0}-1}\right)}};{N_{\mathrm{var}}})=F_{\rho}^{-1}\left(1-\frac{{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}{\color[rgb]{1,0,0}-1}\right)}}}{10^{4}};{N_{\mathrm{var}}}\right).$
(37)

Based on the information from the screening dots, we presume that the interior
of the $N_{\mathrm{var}}$-ball of radius $r$ does not contain any failure
event.

All numerical examples in this paper are analyzed using this strategy. It is
our experience that using IS with density
$h_{\mathrm{ann}}({\boldsymbol{x}})$, hereinafter called “global IS”, provides
a good balance between robustness and efficiency. We remark that in a
hypothetical case, in which the rare event occurs exactly outside a
circle/ball/hyperball with radius $r$ in the standard Gaussian space, the
sampling density $h_{\mathrm{ann}}(x)$ becomes an _optimal_ IS sampling
density (with zero variance of the estimator). The reason is that
$h_{\mathrm{ann}}(x)$ is proportional to the original standard Gaussian
density and yet it removes completely the contribution from the safe region
($N_{\mathrm{var}}$-ball with radius $r$).

Let us now consider the sampling density
$h_{\mathrm{ann}}\left({\boldsymbol{x}};r,R\right)$ introduced in Eq. (23) and
substitute it into Eq. (33)

$\displaystyle{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}$
$\displaystyle\approx\frac{1}{n_{\textsf{IS}}}\sum_{i=1}^{n_{\textsf{IS}}}{\boldsymbol{1}_{\pazocal{T}}^{({N_{\mathrm{sim}}})}\left({\boldsymbol{x}}_{i}\right)}\frac{f_{{\boldsymbol{X}}}({\boldsymbol{x}}_{i})}{h_{\mathrm{ann}}({\boldsymbol{x}}_{i})}=\frac{p_{\mathrm{ann}}}{n_{\textsf{IS}}}\underbrace{\sum_{i=1}^{n_{\textsf{IS}}}{\boldsymbol{1}_{\pazocal{T}}^{({N_{\mathrm{sim}}})}\left({\boldsymbol{x}}_{i}\right)}}_{n_{\textsf{IS},\pazocal{T}}}=p_{\mathrm{ann}}\frac{{n_{\textsf{IS},\pazocal{T}}}}{n_{\textsf{IS}}}.$
(38)

In other words, computation of the likelihood ratio ($f_{{\boldsymbol{X}}}/h$)
at each node is _not_ needed as the Gaussian densities cancel out. This is
because the samples generated via the simple procedure from Sec. 3.3 have
their density proportional to $f_{{\boldsymbol{X}}}$. It suffices to simply
compute the proportion of IS nodes that signaled event $\pazocal{T}$ (e.g.,
failure) and multiply it with $p_{\mathrm{ann}}$. The closer the true event
domain is to the annulus, the closer $h_{\mathrm{ann}}$ is to the optimal IS
density, and so the estimation variance vanishes. With this strategy, the
annulus between the two radii $r$ and $R$ is effectively examined. Whatever
event occurs outside the outer radius $R$ is associated with a negligible
probability.

The variance of such an IS estimation is also simple to obtain. By
substituting $n_{\textsf{IS},\pazocal{T}}$ as the number of samples that lead
to the event $\pazocal{T}$ and $p_{\mathrm{ann}}$ as the constant likelihood
ratio into Eq. (35), the estimation variance becomes

$\displaystyle\text{Var}_{h}[{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}]\approx\frac{1}{n_{\textsf{IS}}}\Bigg{\\{}\left[\frac{n_{\textsf{IS},\pazocal{T}}}{n_{\textsf{IS}}}p^{2}_{\mathrm{ann}}\right]-\left({p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}\right)^{2}\Bigg{\\}}=\frac{1}{n_{\textsf{IS}}}\Bigg{\\{}\left[{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}p_{\mathrm{ann}}\right]-\left({p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}\right)^{2}\Bigg{\\}}=\frac{{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}}{n_{\textsf{IS}}}\left(p_{\mathrm{ann}}-{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}\right),$
(39)

and therefore, by using Eq. (38), the general coefficient of variation in Eq.
(36) simplifies to

$\displaystyle\text{CoV}_{h}[{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}]$
$\displaystyle\approx\frac{1}{\sqrt{n_{\textsf{IS}}}}\sqrt{\frac{p_{\mathrm{ann}}}{{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}}-1}=\frac{1}{\sqrt{n_{\textsf{IS}}}}\sqrt{\frac{n_{\textsf{IS}}}{n_{\textsf{IS},\pazocal{T}}}-1}$
(40)

There can be two reasons for this coefficient of variation of the estimator
potentially being unacceptably high: (i) the number $n_{\textsf{IS}}$ of
integration nodes is small (the square root in the first fraction), and (ii)
the failure region is highly localized due to which which the important ring
contains only a small proportion of nodes associated with the rare event (the
fraction in the second square root is much greater than one). In the first
case, the relevant measure is to increase $n_{\textsf{IS}}$, while in the
second case, it pays to switch to the set of local IS runs around various
$\pazocal{T}$ points (see above), or generally around clusters of rare-event
clouds which can be localized using the “k-means clustering” algorithm.

Fig. 6 presents all the important information about the process of _extension_
of ED (top left) and the _estimation_ of probabilities which can be performed
at any time during the process. The evolution of both is shown in the Wavy
circle video covering the history from the very first limit state function
evaluation up to ${N_{\mathrm{sim}}}=200$. The video frames, which can be
displayed one by one, correspond to individual stages of the process, thus
enabling a detailed inspection of the process. The accuracy of the estimation
is excellent already at ${N_{\mathrm{sim}}}\approx 70$ when all of the seven
“failure regions” become discovered. The blue line plotted in the bottom
diagram shows the amount of the “probability bite $\psi$” occupied by the
neighborhood of the selected candidate. It can be seen that once about
${N_{\mathrm{sim}}}=30$ points have been evaluated, the $\psi$ contributions
become smaller than the rare event probability itself. The decrease in $\psi$
with an increasing ${N_{\mathrm{sim}}}$ can be used, along with the
stabilization of probability estimations, to decide the profitability of
further $g(\boldsymbol{x})$ function evaluations (stopping criterion).

The red line, such as the one visible in Fig. 6, is always accompanied in this
paper by a scatterband of $\pm$ one standard deviation obtained from Eq. (40).
In particular, the two lines parallel the estimate
${{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}}$ are computed as
${{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}}\left(1\pm\text{CoV}_{h}[{p_{\pazocal{T}}^{\left({N_{\mathrm{sim}}}\right)}}]\right)$.
We remark that such a scatterband accounts for the variance in the IS
estimator for a given classification surrogate only. Since we used large
numbers of integration dots (many thousands), the scatterband is very narrow
(almost invisible). However, it must be stressed that Eq. (39) does _not_
account for the variability associated with the varying surrogate classifier
$\boldsymbol{1}_{\pazocal{T}}^{({N_{\mathrm{sim}}})}\left({\boldsymbol{x}}\right)$,
which evolves differently in various runs of the ED extension algorithm. This
means that another run of the proposed extension algorithm may result in
different evolution of the rare event probability estimations. However, our
experience is that in almost all tested cases, the differences in the
evolution of the classifier are negligible and they quickly diminish as the
number of points ${N_{\mathrm{sim}}}$ increases. The only exception is the
“Modified Rastrigin” counterexample (Sec. 7.6), in which the surprise
discoveries of various closed parts of the failure domains are very much
randomly variable.

##  6 Global sensitivities and importance measures of individual variables

One of the by-products of the design point search for FORM/SORM is the simple
evaluation of $\alpha$-sensitivities [87, 88]. These values have a clear
geometrical meaning, and they have immediate application in the partial safety
factor method in design codes. Suppose we have a linear performance function
$g(\boldsymbol{x})$ with one distinct design point
${\boldsymbol{x}^{\star}}=\arg\,\min\\{\lVert{\boldsymbol{x}}\rVert\;|\;g(\boldsymbol{x})=0\\}$.
The Euclidean distance of this most central failure point to the origin of the
standard Gaussian space is the safety index
$\beta=\lVert{\boldsymbol{x}^{\star}}\rVert=\sqrt{{\boldsymbol{x}^{\star}}^{\mkern-1.5mu\mathsf{T}}\cdot{\boldsymbol{x}^{\star}}}=\sqrt{\sum_{v=1}^{{N_{\mathrm{var}}}}(x_{v}^{\star})^{2}}$.
The distance can also be written using the vector of $\alpha$-sensitivities,
${\boldsymbol{\alpha}}$, which has the unit size:
$\lVert{\boldsymbol{\alpha}}\rVert=\sqrt{\sum_{v=1}^{{N_{\mathrm{var}}}}\alpha_{v}^{2}}=1$.
The linear and normalized approximation $M({\boldsymbol{x}})$ to the safety
margin can thus be expressed as
$M({\boldsymbol{x}})={\boldsymbol{\alpha}}^{\mkern-1.5mu\mathsf{T}}({\boldsymbol{x}^{\star}}-{\boldsymbol{x}})=\beta-{\boldsymbol{\alpha}}^{\mkern-1.5mu\mathsf{T}}{\boldsymbol{x}}$.
The variance of $M$ is
$\sqrt{\sum_{v=1}^{{N_{\mathrm{var}}}}\alpha_{v}^{2}}=1$ [88]. Vector
${\boldsymbol{\alpha}}$ is the negative of the gradient of the performance
function at the origin of the standard normal space, i.e.,
${\boldsymbol{\alpha}}$ points in the _important direction_ (see Fig. 7a), and
the coordinates of the design point can be obtained as
${\boldsymbol{x}^{\star}}=\beta{\boldsymbol{\alpha}}$. The components of
${\boldsymbol{\alpha}}$, i.e., the direction cosines of $\boldsymbol{\alpha}$
are called the $\alpha$-sensitivities (or $\alpha$-factors), and they are
regarded as measures of the sensitivity of the reliability index to
inaccuracies in the value of ${\boldsymbol{x}}$ at the design point

$\displaystyle\alpha_{v}=\left.\frac{\partial\beta}{\partial
x_{v}}\right\rvert_{{\boldsymbol{x}^{\star}}}.$ (41)

The squared components of vector ${\boldsymbol{\alpha}}$ can be obtained
directly from the coordinates of the design point:

$\displaystyle\alpha_{v}^{2}=\frac{(x_{v}^{\star})^{2}}{\beta^{2}}=\frac{(x_{v}^{\star})^{2}}{\sum_{v=1}^{{N_{\mathrm{var}}}}(x_{v}^{\star})^{2}},\quad
v=1,\ldots,{N_{\mathrm{var}}}.$ (42)

Figure 7: Measures of importance with respect to failure probability. a)
illustration of FORM $\alpha$-sensitivities constructed using linearization of
a failure surface in the most central failure point
${\boldsymbol{x}^{\star}}$. b)-c) illustration of the proposed global
sensitivity measures $s_{v}^{2}$ by decomposition of ${p_{\pazocal{F}}}$ into
additive contributions by individual variables.

The relative importance measures by FORM/SORM are highly affected by the
accuracy of FORM/SORM solutions, and they may be inadequate for component or
system reliability problems featuring multiple design points or surfaces with
high curvatures. Yet, good importance measures are needed to identify critical
random variables concerning the reliability of a product or a process and
enable focusing on critical uncertainties (disregarding unimportant ones) in
the optimization of the product reliability.

In this paper, we introduce a new generalized importance measure that
considers all (known) points contributing to failure, weighted by the original
density. This measure may use the the performance function’s gradient when
available, but can be defined also for discrete-state performance function for
which the gradient is locally zero. The proposed measure comes as a by-product
of the above-proposed technique for reliability estimations. Consider that the
total unit contribution of _any_ point
${\boldsymbol{x}}=\\{x_{1},\ldots,x_{v},\ldots,{\boldsymbol{x}}_{{N_{\mathrm{var}}}}\\}$
to an event of interest can be decomposed into individual coordinates in the
spirit of the $\alpha$-sensitivities

$\displaystyle\alpha_{v}^{2}({\boldsymbol{x}})=\left(\frac{\color[rgb]{1,0,0}\nabla
g_{v}({\boldsymbol{x}})}{\color[rgb]{1,0,0}\lVert\nabla\boldsymbol{g}({\boldsymbol{x}})\lVert}\right)^{2}=\frac{\color[rgb]{1,0,0}g_{v}^{2}({\boldsymbol{x}})}{\sum_{v=1}^{{N_{\mathrm{var}}}}\color[rgb]{1,0,0}g_{v}^{2}({\boldsymbol{x}})},\quad\text{i.e.}\;\sum_{v=1}^{{N_{\mathrm{var}}}}\alpha_{v}^{2}({\boldsymbol{x}})=1\;\;\text{for
any}\;{\boldsymbol{x}},$ (43)

where $g_{v}({\boldsymbol{x}})$ is the projection of the local gradient
$\lVert\nabla\boldsymbol{g}({\boldsymbol{x}})\lVert$ onto direction $v$. This
gradient is perpendicular to contour lines of the performance function
$g(\boldsymbol{x})$. In FORM approximation, the assumption of linearity of
$g(\boldsymbol{x})$ guarantees that the direction of the gradient is identical
for all points in the design domain and can be taken as the radius vector
connecting the origin with the design point. When $g(\boldsymbol{x})$ is a
general smooth function, the gradient and its components can be computed by
evaluating the partial derivatives (either analytically or numerically). For
discrete-state performance functions treated in this paper, such gradient does
not provide any meaningful information about the steepest descent and must be
replaced by another information available from the previously evaluated points
in the ED. A reasonable way to replace the gradient in an existing
(supposedly) failure point is the direction to the _nearest safe point_. These
directions are perpendicular to the failure surface, see the magenta arrows in
Fig. 7b which illustrate the local gradients at various points in the failure
set $\pazocal{F}$. Technically, these directions are not difficult to compute.
Assume the extension algorithm has already discovered at least one failure
point. The nearest-neighbor classification therefore divides a number of
integration nodes into safe and failure groups. In each failure integration
node, one can now approximate the direction perpendicular to the failure
surface by taking $K$ nearest neighbors from the group of safe nodes. The line
connecting the failure node with the centroid of $K$ nearest safe nodes is a
good approximation of the “gradient” direction, see the empty circles Fig. 7b
representing such centroids. An efficient way to perform the computations is
to employ K-d tree space-partitioning data structure for organizing points.
This data structure enables classification used for fast identification of the
$K$ nearest neighbors from the group of safe nodes. For Python, the
functionality is readily available via sklearn.neighbors.KNeighborsClassifier
in scikit-learn library [89]. The function allows for vectorized call which
quickly finds the $K$ nearest neighbors for a large set of nodes from the
failure set at once.

The resulting squared projections $\alpha_{v}^{2}({\boldsymbol{x}})$ can be
viewed as importance measures of the individual dimensions and can be used as
additive shares of probability density at any point ${\boldsymbol{x}}$. This
choice is natural as $\alpha_{v}$s weigh locally the individual increments
needed to change the system state:
$\alpha_{v}(x)={\color[rgb]{1,0,0}\frac{\partial g({\boldsymbol{x}})}{\partial
x_{v}}}$. Using these shares, the total failure of probability defined in Eq.
(1) can be rewritten as

$\displaystyle{p_{\pazocal{F}}}\equiv\idotsint_{{\pazocal{F}}}\sum_{v=1}^{{N_{\mathrm{var}}}}\alpha_{v}^{2}({\boldsymbol{x}})f_{{\boldsymbol{X}}}({\boldsymbol{x}})\;\mathrm{\,d}\hskip
0.50003pt{\boldsymbol{x}}=\sum_{v=1}^{{N_{\mathrm{var}}}}\idotsint_{{\pazocal{F}}}\alpha_{v}^{2}({\boldsymbol{x}})f_{{\boldsymbol{X}}}({\boldsymbol{x}})\;\mathrm{\,d}\hskip
0.50003pt{\boldsymbol{x}}=\sum_{v=1}^{{N_{\mathrm{var}}}}{p_{\pazocal{F},v}}.$
(44)

In this way, the probability of failure (or analogously any event type) is
obtained as a sum of contributions ${p_{\pazocal{F},v}}$ of individual
variables. Each variable contributes to $p_{\pazocal{F}}$ by

$\displaystyle{p_{\pazocal{F},v}}=\idotsint_{{\pazocal{F}}}\alpha_{v}^{2}({\boldsymbol{x}})f_{{\boldsymbol{X}}}({\boldsymbol{x}})\;\mathrm{\,d}\hskip
0.50003pt{\boldsymbol{x}}=\idotsint_{{\pazocal{F}}}f_{v}({\boldsymbol{x}})\;\mathrm{\,d}\hskip
0.50003pt{\boldsymbol{x}},$ (45)

where we define a part of the standard Gaussian density ascribed to a single
variable $v$ as (see also the illustrations in Figs. 7c and d)

$\displaystyle
f_{v}({\boldsymbol{x}})={\color[rgb]{1,0,0}\alpha^{2}_{v}({\boldsymbol{x}})}f_{{\boldsymbol{X}}}({\boldsymbol{x}}).$
(46)

This local contribution to ${p_{\pazocal{F},v}}({\boldsymbol{x}})$ is easy to
evaluate as it is dependent only on the local gradients at point
${\boldsymbol{x}}$ and its standard Gaussian density. Finally, the shares
${p_{\pazocal{F},v}}$ can be standardized by ${p_{\pazocal{F}}}$ to form the
proposed global importance measures $s_{v}^{2}$ associated with individual
variables

$\displaystyle
s_{v}^{2}=\frac{{p_{\pazocal{F},v}}}{{p_{\pazocal{F}}}},\quad\text{i.e.}\;\sum_{v=1}^{{N_{\mathrm{var}}}}s_{v}^{2}({\boldsymbol{x}})=1.$
(47)

An important aspect is that the proposed $s_{v}^{2}$s are not based on values
of the performance function, as only the binary information indicating an
_event_ is needed. This is an important property because the importance
measures should _not_ be dependent on the way a performance function is
defined if it provides the same failure boundary. A robust importance measure
for sensitivity to an event should be invariant under reformulations or
reparametrizations of the underlying problem.

The numerical values of $s_{v}^{2}$ generally do not match with the classical
factors $\alpha_{v}^{2}$ defined in Eq. (42). The reason is that the proposed
measures $s_{v}^{2}$ weigh the contributions to $p_{\pazocal{F}}$ from the
whole failure set $\pazocal{F}$, while the classical $\alpha$-factors, in
fact, consider only a single design point. The dark red line in Fig. 8 right
presents the classical $\alpha$-sensitivities for a linear failure surface in
two dimensions in dependence on its rotation around the origin (see the dashed
lines in Fig. 8 middle). The meaning of the proposed measure is different from
the $\alpha$-sensitivities. We argue that it provides meaningful results also
for complicated failure domains (nonlinear, non-smooth, disconnected failure
sets, etc.).

To show an example in which the results are quite different from
$\alpha$-sensitivities, we consider the following nonlinear bivariate
performance function: $g(\boldsymbol{x})=\beta-x_{1}^{4}/c-x_{2}$, where
$\beta=3$. In order to have a unique most central failure point, the constant
$c$ must be greater than 32; we take $c=33$. Since $g(\boldsymbol{x})$ is
linear in $x_{2}$ and constant in $x_{1}$ when $x_{1}=0$, the gradient point
search initiated at the origin easily finds the unique design point
${\boldsymbol{x}^{\star}}=\\{0,3\\}$; see the point trajectory in Fig. 8 left.
The classical $\alpha$-sensitivities are thus
${\boldsymbol{\alpha}}=\\{0,1\\}$. These are the sensitivities related to the
unique $\beta$ point: the failure surface $\partial{\pazocal{F}}$ is linear at
${\boldsymbol{x}^{\star}}$, and the SORM correction fails because the
curvature in the $\beta$-point vanishes. However, it is evident that there are
many highly probable regions with nonzero $x_{1}$ which considerably
contribute to the failure probability. In particular, there are two
symmetrically distributed regions around points
$x_{1}=\pm{\color[rgb]{1,0,0}2.84}$, $x_{2}=1.02$; see the red regions in Fig.
8 left). Indeed, the proposed global importance measures are not much
different from each other: $s_{1}^{2}\approx{\color[rgb]{1,0,0}0.57}$ and
$s_{2}^{2}\approx{\color[rgb]{1,0,0}0.43}$. The fact is that
$s_{1}^{2}>s_{2}^{2}$ may be seen as contradicting the FORM sensitivity
$\alpha^{2}_{1}=0$, i.e., the zero importance of variable $x_{1}$.

Figure 8: Left: the nonlinear performance function and the design point
search. Middle: failure surfaces for various rotation angles $\gamma$. Right:
Comparison of the classical $\alpha^{2}$-sensitivities with the proposed
importance measures $s_{v}^{2}$ for various rotations of the problem by angle
$\gamma$.

Consider now a clockwise rotation of the coordinate system by an angle
$\gamma$. The failure surface revolves as a rigid body and moves the
$\beta$-points along a circle of radius $\beta=3$; see Fig. 8 middle, and also
the case of $\gamma=30\degree$, which is visualized in Fig. 7. The classical
$\alpha$-sensitivities become strongly dependent on the angle:
${\boldsymbol{\alpha}}=\\{\sin^{2}(\gamma),\cos^{2}(\gamma)\\}$; see Fig. 8
right. The same figure documents that the proposed global importance measures
are almost insensitive to the angle of rotation, because it considers all
three failure regions.

An extreme case documenting the fundamental differences between the local and
the proposed global importance measure is the case of a rotationally symmetric
failure domain in standard Gaussian space (the exterior of the
$N_{\mathrm{var}}$-ball) has all importance measures identical:
${p_{\pazocal{F},v}}=1/{N_{\mathrm{var}}}$, which expresses the fact that no
variable is more important than any other for the achievement of failure. The
classical $\alpha$-sensitivities would, in this case, take arbitrary values in
order to make their squared sum equal to one, i.e., any point on the failure
surface $\partial{\pazocal{F}}$ of the $N_{\mathrm{var}}$-ball can be
selected.

The evaluation of individual ${p_{\pazocal{F},v}}$ (and therefore also
$s_{v}^{2}$) is very cheap as it can be seen as a by-product of the sampling
analysis employed to deliver an estimation of ${p_{\pazocal{F}}}$. Suppose we
have an existing set of points that were sampled proportionally to the
standard Gaussian density $f_{{\boldsymbol{X}}}({\boldsymbol{x}})$. It can be
the global importance sample with $N_{\mathrm{sim}}$ points that was obtained
outside the $N_{\mathrm{var}}$-ball or in the annuloidal
$h_{\mathrm{ann}}({\boldsymbol{x}})\propto
f_{{\boldsymbol{X}}}({\boldsymbol{x}})$. From this sample, we only select a
vector ${\boldsymbol{x}}$ containing ${n{F}}$ points marked as “failure” (or
any other event type, based on the nearest neighbor surrogate). All these
points are equally probable, and they each represent the same share of the
estimated failure probability
${p_{\pazocal{F}}}({\boldsymbol{x}}_{i})={p_{\pazocal{F}}}/{n{F}}$,
$i=1,\ldots,{n{F}}$. This share can be further split into individual
directions in the spirit of Eq. (43) by selecting the nearest safe point to
the failure point ${\boldsymbol{x}}_{i}$, and computing the $v$ squared
standardized projections of such a pair at point ${\boldsymbol{x}}_{i}$:
${\color[rgb]{1,0,0}\alpha^{2}_{v}({\boldsymbol{x}}_{i})}$. In this way, the
contribution of the $i$th point in the $v$th direction reads

$\displaystyle{p_{\pazocal{F},v}}({\boldsymbol{x}}_{i})={p_{\pazocal{F}}}({\boldsymbol{x}}_{i}){\color[rgb]{1,0,0}\alpha^{2}_{v}({\boldsymbol{x}}_{i})}=\frac{{p_{\pazocal{F}}}}{{n{F}}}{\color[rgb]{1,0,0}\alpha^{2}_{v}({\boldsymbol{x}}_{i})},\quad\text{i.e.,
each point contributes
}\;{p_{\pazocal{F}}}({\boldsymbol{x}}_{i})=\sum_{v=1}^{{N_{\mathrm{var}}}}{p_{\pazocal{F},v}}({\boldsymbol{x}}_{i})=\frac{{p_{\pazocal{F}}}}{{n{F}}}\;\text{for
any}\;{\boldsymbol{x}}_{i},$ (48)

where $x_{i,v}$ is the $v$th coordinate of point ${\boldsymbol{x}}_{i}$. The
desired estimation of importance measures $s_{v}^{2}$ of individual directions
can be obtained by summation over ${n{F}}$ failing nodes with a fixed
direction index $v$ and dividing by the failure probability

$\displaystyle
s_{\pazocal{F},v}^{2}\approx\sum_{i=1}^{{n{F}}}\frac{{p_{\pazocal{F},v}}({\boldsymbol{x}}_{i})}{{p_{\pazocal{F}}}}=\frac{1}{{n{F}}}\sum_{i=1}^{{n{F}}}{\color[rgb]{1,0,0}\alpha^{2}_{v}({\boldsymbol{x}}_{i})}.$
(49)

Sensitivity to any other event type $\pazocal{T}$ is obtained analogously by
retaining only ${n_{\pazocal{T}}}$ points corresponding to that event and
computing the average:
$s_{\pazocal{T},v}^{2}\approx\sum_{i=1}^{{n_{\pazocal{T}}}}\left[{\color[rgb]{1,0,0}\alpha^{2}_{v}({\boldsymbol{x}}_{i})}\right]/{{n_{\pazocal{T}}}}$.
Therefore, the proposed sensitivity measure is just a cheap by-product of the
proposed method. In cases when the sampling probability is not proportional to
the standard Gaussian density, it is no longer true that all points have the
same contribution of $1/{{n_{\pazocal{T}}}}$, and therefore straightforward
re-scaling in a similar manner to importance sampling must be employed

$\displaystyle
s_{\pazocal{T},v}^{2}\approx\frac{1}{p{T}}\frac{1}{n_{\textsf{IS}}}\sum_{i=1}^{n_{\textsf{IS}}}{\boldsymbol{1}_{\pazocal{T}}^{({N_{\mathrm{sim}}})}\left({\boldsymbol{x}}_{i}\right)}{}\frac{f_{{\boldsymbol{X}}}({\boldsymbol{x}}_{i})}{h({\boldsymbol{x}}_{i})}{\color[rgb]{1,0,0}\alpha^{2}_{v}({\boldsymbol{x}}_{i})}$
(50)

##  7 Numerical examples

We present a variety of numerical examples which have been selected to explore
different classes of problems posing unique challenges. In the first seven
examples, two-dimensional problems are defined in the space of independent
standard Gaussian random variables; see Fig. 9 for a quick overview of the
selected functions. After that, an engineering example with Gaussian inputs
shows the applicability of the proposed method to real nonlinear computational
mechanics problem. The next example remains bivariate, however, it documents
the straightforward applicability of the approach to problems with correlated
non-Gaussian variables. Then, problems in higher dimensions are analyzed. In
all definitions of the functions, we present the expressions that return
smoothly or non-smoothly varying output variables. However, the proposed
extension algorithm receives categorical information only (such as binary
“failure-success” codes). The same holds for the estimation, which uses only
the indicator functions signaling an event. The only exceptions are “Four
Branch” and “Metaballs” examples for which we also examine the degree of
improvement in probability estimation when a smooth interpolation of the
point-wise information in the ED via the Radial Basis Function is employed as
the classifier.

Figure 9: Overview of seven two-dimensional examples used for demonstration of
the proposed method. Design points are visualized using solid black balls.

###  7.1 Wavy circle

The failure surface used to demonstrate the _extension_ and _estimation_ steps
in Fig. 6 is a sine wave of amplitude $a=1$ superposed onto a circle of
average radius $r_{\mathrm{ave}}=4$. The corresponding limit state function
can be formulated as, e.g.,

$\displaystyle g\left(x_{1},x_{2}\right)$
$\displaystyle=r_{\mathrm{ave}}+a\sin(k\varphi_{{\boldsymbol{x}}})-\lVert{\boldsymbol{x}}\rVert,$
(51)

where $\lVert{\boldsymbol{x}}\rVert=\sqrt{x_{1}^{2}+x_{2}^{2}}$ is the
Euclidean distance from the origin,
$\varphi_{{\boldsymbol{x}}}=\operatorname{atan_{2}}(x_{2},x_{1})=\arctan(x_{2}/x_{1})$
is the angle of vector ${\boldsymbol{x}}$ from axis $x_{1}$, and $k=7$ is the
number of design points. The sum of the first two terms represents the maximum
safe distance from the origin, and the rare event (failure) occurs “behind”
the wavy boundary, i.e., when $g\leq 0$. The presented algorithm uses only
binary information (failure or success). The seven most central failures are
located at an identical distance of $r_{\mathrm{ave}}-a=3$. Fig. 6 presents
the situation after ${N_{\mathrm{sim}}}=100$ evaluations of the limit state
function, and the full history starting with the initial evaluation at the
origin is included in the Wavy circle video. The video captures the
progressive (i) discovery of all the design points and (ii) the refinement of
the failure surface, both alternating based on ranking via the $\psi$
criterion to occur proportionally to the Gaussian density
$f_{{\boldsymbol{X}}}$. The estimation of failure probability converges
towards the exact failure probability result ${p_{\pazocal{F}}}\approx
2.582\cdot 10^{-3}$.

The _global_ IS performed on the nearest-neighbor surrogate provides stable
estimations in the ring between $r=3$ and $R\approx 5.2$, see Fig. 6. The
proposed global sensitivities estimated via Eq. (47) are all identical
$s_{\pazocal{F},1}^{2}=s_{\pazocal{F},2}^{2}=0.5$, which corresponds to the
distribution of seven failure regions, while the local FORM-style values of
${\boldsymbol{\alpha}}$ sensitivities are very different and depend on which
design point is considered (e.g., ${\boldsymbol{\alpha}}=\\{0,1\\}$ for the
point ${\boldsymbol{x}^{\star}}=\\{0,3\\}$).

We used OpenTURNS software [90] to run various techniques for failure
probability estimation. The gradient-based design point search for FORM
analysis initiated at a random location needs, after some help, 33 model
evaluations to discover one of the design points, and thus FORM approximates
the probability incorrectly as $\Phi(-3)\approx 1.35\cdot 10^{-3}$. The SORM
analysis increases the number of limit state function calls because it
computes the failure surface curvature about the design point and decreases
the estimation by about $3.2\cdot 10^{-4}$ (the values slightly differ
depending on the method used: Tvedt, Hohenbichler, and Breitung).

The results obtained with SuS implemented in OpenTURNS very much depend on the
computational budget. If the number of limit state function evaluations
reserved for each probability level is sufficiently high (a few thousand), the
stochastic gradient optimization is able to correctly locate all seven failure
regions and estimate the failure probability accurately. When, however, the
total number of function calls drops below roughly one thousand, the estimate
becomes incorrect. Similar statements are true about results obtained from the
Adaptive Directional Stratification Algorithm in OpenTURNS: the number of
function evaluations must be in the thousands for the method to provide
sufficiently good results. For this function, the proposed technique provides
better efficiency than the above methods because it provides more accurate
estimations with considerably fewer function calls.

Figure 10: The “Wavy line” problem in Eq. (52). The complete evolution of all
panels is available as Wavy line video.

The picture changes completely when the methods are supplied with binary
information about $g(\boldsymbol{x})$ only. Design point search, SuS,
directional sampling, etc., break down completely, and no efficiency/accuracy
comparison is possible with the proposed technique. The same is true also for
various tail approximation methods, which are based on fitting various
distribution functions to the supposedly smooth output variable. This is true
for all functions used in the remainder of this paper, and therefore, a
comparison of the presented algorithm with SuS or directional sampling will
not be presented in most of the examples.

The only method class from the ones available from the wide selection in
OpenTURNS is the class of pure sampling methods, such as crude Monte Carlo, or
IS around the mean values. These methods are able to operate given the
dichotomous nature of information: they count the number of failed samples. Of
course, the efficiency is much worse compared to both proposed IS techniques
(global and local), which involve sampling such that the fraction of failure
samples of all function evaluations is much higher. For example, a single run
with ${N_{\mathrm{sim}}}=1\,000$ of the standard IS around mean values, and
with a Gaussian sampling density whole standard deviations are equal to three,
the ${p_{\pazocal{F}}}$ estimation is around 0.003, and the cov of the
estimator is 0.15, which is quite high.

###  7.2 Wavy line

The previous example revealed that the algorithm distributes attention evenly
among all seven failure regions, which share the same distance from the
origin. Another test featuring multiple design points documents the ability to
treat multiple “design points” with different contributions to the failure
probability. Inspired by Sundar and Shields [34], who presented a similar
problem to illustrate their construction of a surrogate model based on an
artificial neural network, we define the performance of the system as

$\displaystyle g\left(x_{1},x_{2}\right)$
$\displaystyle=-\frac{1}{4}x_{1}-x_{2}+\sin(5x_{1})+5.5,$ (52)

where the constants were selected such that the failure surface is a sine
function (with its frequency controlled by the multiplier 5) superposed onto
linearly decreasing function (with a slope of $-1/4$), which is shifted in the
direction of $x_{2}$ by a constant 5.5. The example is a low-dimensional
problem with a very strongly nonlinear limit state possessing multiple
effectively disjoint failure domains, see Fig. 10.

FORM analysis is based on the single $\beta$-point
${\boldsymbol{x}^{\star}}=\\{0.943626,4.26411\\}$, whose distance from the
origin is $\beta=4.3672719$. Therefore, the FORM solution estimates the
failure probability as $\Phi(-\beta)=6.29\cdot 10^{-6}$. The proposed method
quickly converges to the correct solution ${p_{\pazocal{F}}}=1.217\cdot
10^{-6}$, see the wavy line video associated with Fig. 10.

The estimated global sensitivities highlight the role of variable $x_{1}$:
$s_{\pazocal{F},1}^{2}\approx 0.829$ and $s_{\pazocal{F},2}^{2}\approx 0.171$.
This is counter-intuitive and very much different from the reasonable  FORM
sensitivities based on ${\boldsymbol{x}^{\star}}$ and thus highlighting the
role of variable $x_{2}$: ${\boldsymbol{\alpha}}=\\{0.0467,0.9533\\}$.  The
reason for this discrepancy is that the performance function is taken as
binary and therefore the information about the linear gradient of
$g(\boldsymbol{x})$ in Eq. (52) is not available. The algorithm approximates
the local gradient at any point as the direction to the nearest safe point,
i.e., perpendicularly to the wavy failure surface. Most of these directions
are close to horizontal (almost parallel with $x_{1}$) and due to this
prevailing horizontal projection, the sensitivity to the first variable
appears to be dominant.

###  7.3 Metaball function – complicated topology

Breitung [54], inspired by [91], used the Metaball example to document the
risk that SuS does not move its stochastic gradient optimization towards the
region with the highest contribution to failure probability. The particular
definition of the Metaball function in [54] reads

$\displaystyle
g\left(x_{1},x_{2}\right)=\frac{30}{\left(\frac{4\left(x_{1}+2\right)^{2}}{9}+\frac{x_{2}^{2}}{25}\right)^{2}+1}+\frac{20}{\left(\frac{\left(x_{1}-2.5\right)^{2}}{4}+\frac{\left(x_{2}-0.5\right)^{2}}{25}\right)^{2}+1}-5$
(53)

Failure occurs when $g(\boldsymbol{x})<0$ and the associated failure
probability reads ${p_{\pazocal{F}}}\approx 1.12857\cdot 10^{-5}$. The
proposed method quickly converges to this result once the first failure is
hit, see Fig. 11. When the ED contains about 40 sequentially added points via
the $\psi$ criterion, the estimation of failure probability using the crude
nearest neighbor surrogate is sufficient. The “probability bites” approximated
by the $\psi$ criterion become more than one order of magnitude less than
${p_{\pazocal{F}}}$ when ${N_{\mathrm{sim}}}>100$, see the blue line. Further
extension of the ED leads to even better refinement of the true boundary and,
therefore, also more accurate probability estimation; see Fig. 11 and its
evolution captured by the MetaBalls video.

Figure 11: The “MetaBalls” problem in Eq. (53). The complete evolution of all
panels is available in the MetaBalls video.

The estimation method presented in this paper uses the simple nearest neighbor
classification surrogate. However, other classifiers such as Support Vector
Machine, Radial Basis Function, Gaussian process, Random Forest, Neural
Networks, Naive Bayes, AdaBoost, Quadratic Discriminant Analysis, Polynomial
Chaos, etc., can be employed to increase the accuracy of the estimation. We
now examine how the estimation accuracy increases by employing radial basis
function (RBF) interpolation built using the true $g(\boldsymbol{x})$ values,
not just the binary output. We have tested the simplest RBF with default
parameters (class Rbf from the scipy.interpolate sub-package [73]) to perform
classification into failure and safe domains based on the sign of the RBF
surrogate. For each sample size, the RBF interpolator was set up based on the
current ED. The global IS estimation was then applied for each sample size to
estimate $p_{\pazocal{F}}$. It leads to a stable almost exact estimate already
for ${N_{\mathrm{sim}}}\geq 60$; see the brown line in Fig. 11. The third
panel from the left in the top row of the figure and also in the associated
MetaBalls video captures the evolution of classification based on smooth RBF
interpolation. The involvement of the smooth interpolator indeed slightly
improves the accuracy and stability of the IS estimate. However, in this
particular example, the estimate based on the simple nearest neighbor
surrogate provides almost equally good results, despite the fact that only the
binary information about $g(\boldsymbol{x})$ is being used.

We now discuss the role of the complicated contour of this particular function
(see Fig. 9) for gradient-based optimization methods such as FORM and SuS.
There are four distinct “design points” for this function (see the white
diamonds in the third panel from the left in Fig. 11): the most central
failure point is
${\boldsymbol{x}^{\star}_{[}}1]\approx\\{-4.26376,-0.00067\\}$ with the
distance from the origin $\beta_{1}\approx 4.26376$, and the other three read:
${\boldsymbol{x}^{\star}_{[}}2]\approx\\{5.1265,-0.2256\\}$ with
$\beta_{2}\approx 5.131$,
${\boldsymbol{x}^{\star}_{[}}3]\approx\\{0.399,-5.241\\}$ with
$\beta_{1}\approx 5.257$, and
${\boldsymbol{x}^{\star}_{[}}4]\approx\\{0.201,5.895\\}$ with
$\beta_{4}\approx 5.898$. The corresponding four FORM approximations are
$p_{1}=\Phi(-\beta_{1})=1.005\cdot 10^{-5}$, $p_{2}=1.437\cdot 10^{-7}$,
$p_{3}=7.336\cdot 10^{-8}$, and $p_{4}=1.838\cdot 10^{-9}$, which reveal that
the first design point is the only decisive one, and its discovery and usage
in FORM approximation leads to very accurate probability estimate (12% error
only).

However, the most central failure point ${\boldsymbol{x}^{\star}_{[}}1]$ is
hard to discover because the path to it from around the origin leads over a
“high peak”. Therefore, the standard gradient-based search algorithms do not
discover it, and instead slide towards ${\boldsymbol{x}^{\star}_{[}}3]$. They
then estimate the $p_{\pazocal{F}}$ three orders of magnitude wrong as
$p_{3}$. As argued by Breitung [54], SuS, which is in fact a stochastic
gradient method, also struggles to move the point cloud towards
${\boldsymbol{x}^{\star}_{[}}1]$. Indeed, we confirm that using excessively
small sample sizes in individual levels of the SuS algorithm implemented in
OpenTURNS [90] makes the point cloud descend from the saddle point
corresponding to the coordinate origin down the depression towards
${\boldsymbol{x}^{\star}_{[}}3]$. The probability estimate is wrong because
the SuS algorithm was unable to navigate the cloud towards
${\boldsymbol{x}^{\star}_{[}}1]$ when the number of $g(\boldsymbol{x})$
evaluations dropped below several thousand. This behavior was also visually
reported by Breitung [54], who used it to show that decisive regions, located
relatively close to the origin, i.e., with high failure probability content,
can be overshadowed by other, less important regions in SuS. The reason is
that for small sample sizes, the sequences of point clouds orient the downhill
search according to basically local information only.

The Directional Sampling implemented in OpenTURNS was also unable to provide a
correct result with a small sample size. The method provides good and stable
results, but the number of limit state evaluations must be in the thousands.

The proposed global sensitivities equal $s_{\pazocal{F},1}^{2}=0.994$ and
$s_{\pazocal{F},2}^{2}=0.006$ underlying the fact that all the “gradients” in
high density regions of $\pazocal{F}$ (i.e. in the vicinity of points
${\boldsymbol{x}^{\star}_{[}}1]$ and ${\boldsymbol{x}^{\star}_{[}}2]$) are
almost perfectly aligned with $x_{1}$. The FORM sensitivities based on
${\boldsymbol{x}^{\star}_{[}}3]$ are
${\boldsymbol{\alpha}}\approx\\{0.0058,0.9942\\}$ incorrectly putting
significance to variable $x_{2}$. This is a consequence of not discovering the
much more important point ${\boldsymbol{x}^{\star}_{[}}1]$ for which the FORM-
style sensitivities would be ${\boldsymbol{\alpha}}\approx\\{1,0\\}$.

###  7.4 Four Branches problem

Figure 12: The “Four Branches” problem in Eq. (58). The complete evolution of
all panels is shown in the Four Branches video.

The following 2-dimensional “Four Branch function” [92] is a common benchmark
problem in reliability analyses; see, e.g., studies with various parameter
settings [93, 48, 26, 94]. The function describes the failure of a  _series
system_ with four distinct limit state components: two linear and two
nonlinear branches of the failure surface. The limit state function reads

$\displaystyle g\left(x_{1},x_{2}\right)$
$\displaystyle=\min\left\\{\begin{array}[]{ll}3+0.1\left(x_{1}-x_{2}\right)^{2}-{\left(x_{1}+x_{2}\right)}/{\sqrt{2}}\\\
3+0.1\left(x_{1}-x_{2}\right)^{2}+{\left(x_{1}+x_{2}\right)}/{\sqrt{2}}\\\
x_{1}-x_{2}+7/\sqrt{2}\\\ x_{2}-x_{1}+7/\sqrt{2}\end{array}\right..$ (58)

The failure event is defined as $g(\boldsymbol{x})\leq 0$. Various authors
have used various parameters replacing the number $7$ (originally 3.5,
sometimes 6, or 11). In our definition, there are two pairs of design points:
two points when $x_{1}=x_{2}=\pm 3\sqrt{2}/2$ at a distance $\beta_{1,2}=3$
and another two when $x_{1}=-x_{2}=\pm 7\sqrt{2}/4$ at a distance
$\beta_{1,2}=3.5$. The exact result failure probability is
${p_{\pazocal{F}}}\approx 2.222\cdot 10^{-3}$, and the proposed global
sensitivities are identical
($s_{\pazocal{F},1}^{2}=s_{\pazocal{F},2}^{2}=0.5$) due to symmetry reasons.

Fig. 12 shows the four blue lines forming the failure surface. The associated
Four Branches video demonstrates the behavior of the proposed method by adding
points one by one. The extension algorithm refines the boundary proportionally
to the probability density featured in the $\psi$ criterion. The consequence
of this is that the classification close to the four remote intersections of
failure surfaces is not performed correctly for small designs. While this is
no problem for the global IS probability estimation, which focuses on the
high-density regions, it is a source of erroneous estimation for the local IS.
As can be seen, the magenta dash-dot line in Fig. 12 represents wrong and
unstable results, which are degraded due to the accentuation of the
inaccurately classified regions (heavily covered by the local IS integration
nodes). The global IS applied to the same binary surrogate classifier provides
accurate and stable results for EDs as small as ${N_{\mathrm{sim}}}\approx
80$. The top right panel in Fig. 12 shows that the corner regions are not
important for the global IS.

We have examined the simplest Gaussian RBF-based classifier in the same way as
we did in the preceding example. Using the global IS estimation leads to
stable, almost exact results which are already in the range of
${N_{\mathrm{sim}}}\in(20,80)$, see the brown line in Fig. 12. The improvement
is due to the more accurate classification of the failure surface in the
vicinity of the four design points; the “corners” were not classified
correctly neither by RBF nor by the binary surrogate. The efficiency of the
RBF classifier is as good as the efficiency of the best methods known in the
literature. Tab. 1 compares many advanced techniques applied to the Four
Branches example, along with the numbers of limit state evaluations
($N_{\mathrm{sim}}$) and the rare event probability estimates. Many sampling
techniques, such as the sequential importance sampling employed in [26], use
too many limit state function calls, and good results are obtained only when
building a surrogate model. All the presented techniques, however, use the
numerical values for the $g(\boldsymbol{x})$. We repeat again that the
proposed algorithm achieves almost the same efficiency; however, it does so
using only the categorical information about $g(\boldsymbol{x})$, which makes
the proposed method very robust.

Table 1: Summary of results for the Four Branches problem obtained with other methods from the literature. Year | Method | ${N_{\mathrm{sim}}}$ | ${p_{\pazocal{F}}^{\left({N_{\mathrm{sim}}}\right)}}\,(\cdot 10^{3})$ | $\mathrm{cov}\left({p_{\pazocal{F}}^{\left({N_{\mathrm{sim}}}\right)}}\right)\left[\%\right]$
---|---|---|---|---
2000 | DS [95] | 227 | $2.2556$ | 37
2011 | AK-MCS+U [48] | 96 | $2.233$ | -
2011 | AK-MCS+EFF [48] | 101 | $2.232$ | -
2011 |  2SMART [39] | 1035 | $2.21$ | 1.7
2013 | CE-AIS-GM [96] | 3 943 | $2.15$ | 3
2014 | MetaAK-IS2 [97] | 48+90 | $2.22$ | 1.7
2016 | AK–SS [98] | 45 | $2.233$ | 4.94
2017 | KRA [99] | 116 | $2.22$ | 4.7
2017 | ASVM-MCS [41] | 89 | $2.13$ | 2.2
2018 | iRS [100] | 33 | $2.24$ | -
2019 | AKEE-SS [101] | 41.7 | $2.20$ | 3.07
2020 | BSC+RLCB [102] | 37 | $2.240$ | 4.89
2020 | DRL [103] | 2597 | $2.314$ | -
2021 | ABSVR1 [104] | 30 | $2.214$ | -
2021 | ABSVR2 [104] | 43 | $2.222$ | -
2021 | SuS+K [105] | 24 | $2.234$ | 2.05
2021 | RVM [106] | 73 | $2.218$ | 2.14
2021 | SVM [107] | 92 | $2.13$ | 10
2022 | RVM-MIS [108] | 161 | $2.252$ | 2.759
2022 | CE-DIS [109] | 275 | $2.20$ | 4.79
2022 | APCK-PDEM [110] | 35 | $2.226$ | 0.32

###  7.5 Black Swan

Consider now a highly localized failure region characterized by the
_simultaneous_ violation of two thresholds: failure $\pazocal{F}$ occurs when
$x_{1}>2$ and $x_{2}>5$. This failure region emerges for the following simple
function studied in chapter 5 of a book [32], where the authors discuss the
difficulties of SuS reaching the failure domain:

$\displaystyle g\left(x_{1},x_{2}\right)$
$\displaystyle=\left\\{\begin{array}[]{ll}5-x_{1},&\text{for }x_{1}\leq 2\\\
5-x_{2},&\text{for }x_{1}>2\end{array}\right.$ (61)

With this definition of the limit state function, the failure region described
above corresponds to $g(\boldsymbol{x})<0$. The exact solution is the product
of the probabilities obtained from two simple FORM-like solutions:
${p_{\pazocal{F}}}=\Phi(-5)\Phi(-2)\approx 6.52136\cdot 10^{-9}$. In this
example, we show that it can be very hard to hit the “black swan” event. The
exploration set must have at least one point located in the rare event domain
(failure region); without it, the discovery of the event is not possible.
However, once a failure event is localized, the proposed extension algorithm
is very effective in the refinement of the small part of the failure surface
which is associated with a high probability.

Figure 13: The “Black Swan” problem in Eq. (61). The complete evolution of all
panels is shown in the Black Swan video.

We generally recommend generating many layers of the predefined exploration
set, covering remote territories in the input space. Unnecessarily distant
points will not be selected while refinement is in progress exploiting more
important regions. However, without offering remote exploration points, the
discovery of new disjoint failure regions would not be possible.

The estimated global sensitivities
$s_{\pazocal{F},1}^{2}\approx{\color[rgb]{1,0,0}0.3189}$ and
$s_{\pazocal{F},2}^{2}\approx{\color[rgb]{1,0,0}0.6811}$ proposed in this
paper are somewhat different from the FORM-style values
${\boldsymbol{\alpha}}$: $2^{2}/(2^{2}+5^{2})\approx 0.1379$ and
$5^{2}/(2^{2}+5^{2})\approx 0.862$, which are based on the single most central
failure point ${\boldsymbol{x}^{\star}}=\\{2,5\\}$. The proposed sensitivities
are based on the whole failure region in which the directions to the nearest
failure surface approximating the missing gradient information are taken as
either parallel or perpendicular to $x_{1}$. In this way, the failure region
${\pazocal{F}}:\\{x_{1}>2\text{ and }x_{2}>5\\}$ becomes split into two
triangles separated by a straight line $x_{2}=x_{1}+3$. These two triangles
contribute separately to ${p_{\pazocal{F}}}=p_{1}+p_{2}$, where
$p_{1}=\int_{x_{1}=2}^{\infty}\left(\int_{x_{2}=x_{1}+3}^{\infty}f_{{\boldsymbol{x}}}(x_{1},x_{2})\mathrm{\,d}\hskip
0.50003ptx_{2}\right)\mathrm{\,d}\hskip 0.50003ptx_{1}\approx 2.0794\cdot
10^{-9}$, leading to $s_{\pazocal{F},1}^{2}=p_{1}/{p_{\pazocal{F}}}$ and
$s_{\pazocal{F},2}^{2}=({p_{\pazocal{F}}}-p_{1})/{p_{\pazocal{F}}}$.

###  7.6 Modified Rastrigin – scattered small closed failure regions

As discussed already in Sec. 4.4, the proposed algorithm may struggle when the