# Active Learning-based Domain Adaptive Localized Polynomial Chaos Expansion

Lukáš Novák<EMAIL_ADDRESS>Brno University of Technology, Brno, Czech
Republic  Michael D. Shields<EMAIL_ADDRESS>Johns Hopkins
University, Baltimore, USA  Václav Sadílek<EMAIL_ADDRESS>Miroslav
Vořechovský<EMAIL_ADDRESS>Brno University of Technology, Brno, Czech
Republic

###### Abstract

The paper presents a novel methodology to build surrogate models of
complicated functions by an active learning-based sequential decomposition of
the input random space and construction of localized polynomial chaos
expansions, referred to as domain adaptive localized polynomial chaos
expansion (DAL-PCE). The approach utilizes sequential decomposition of the
input random space into smaller sub-domains approximated by low-order
polynomial expansions. This allows approximation of functions with strong
nonlinearties, discontinuities, and/or singularities. Decomposition of the
input random space and local approximations alleviates the Gibbs phenomenon
for these types of problems and confines error to a very small vicinity near
the non-linearity. The global behavior of the surrogate model is therefore
significantly better than existing methods as shown in numerical examples. The
whole process is driven by an active learning routine that uses the recently
proposed $\Theta$ criterion to assess local variance contributions [1]. The
proposed approach balances both _exploitation_ of the surrogate model and
_exploration_ of the input random space and thus leads to efficient and
accurate approximation of the original mathematical model. The numerical
results show the superiority of the DAL-PCE in comparison to (i) a single
global polynomial chaos expansion and (ii) the recently proposed stochastic
spectral embedding (SSE) method [2] developed as an accurate surrogate model
and which is based on a similar domain decomposition process. This method
represents general framework upon which further extensions and refinements can
be based, and which can be combined with any technique for non-intrusive
polynomial chaos expansion construction.

###### keywords:

Polynomial Chaos Expansion , Adaptive Sampling , Sequential Sampling , Local
Approximations , Active Learning , Stochastic Spectral Embedding

††journal: Computer Methods in Applied Mechanics and Engineering

Effective construction of a general purpose surrogate model based on
polynomial chaos expansion.

Novel method for sequential decomposition of the input random space and
construction of local approximations.

Sequential domain decomposition and sample size extension based on an active
learning methodology.

Active learning is represented by variance-based $\Theta$ criterion developed
for polynomial chaos expansion.

##  1 Introduction

The Polynomial Chaos Expansion (PCE), originally proposed by Norbert Wiener
[3] and further investigated in the context of engineering problems by many
researchers, e.g. [4, 5], is a preferred method for uncertainty quantification
(UQ) and surrogate modeling in industrial applications [6, 7] thanks to its
efficiency and powerful post-processing. Once a PCE is available for a given
problem, the constructed explicit function can be exploited to directly
estimate important properties of the original problem including its
statistical moments, response probability distribution or sensitivity indices
(without additional sampling [8]), which brings significant efficiency for
surrogate modeling, sensitivity analysis, uncertainty quantification and
reliability analysis [9].

The PCE, in its non-intrusive form, offers a convenient way to perform
probabilistic analysis of any black-box model, e.g. finite element models
representing complex physical systems in engineering. There are generally two
types of non-intrusive methods to calculate the deterministic PCE
coefficients: spectral projection and linear regression. The spectral
projection approach utilizes the orthogonality of the multivariate polynomials
and calculates the coefficients using inner products. The spectral projection
leads to an explosion of computational complexity referred to as the _curse of
dimensionality_. Therefore, the non-intrusive approach based on linear
regression is often preferred. Although it is typically less expensive than
the spectral projection (the number of samples should be at least
$\mathcal{O}(P\,\ln(P))$, where $P$ is the number of terms in the PCE [10,
11]), it suffers from the _curse of dimensionality_ as well, since the number
of PCE terms grows rapidly with both dimension and maximum polynomial order.
Therefore, it becomes necessary to employ advanced adaptive techniques to
construct sparse PCEs that yield efficient solutions for real-world physical
systems.

Regression-based PCE can be significantly affected by the selected sampling
scheme, as was recently shown in an extensive review paper [12] comparing
several general statistical sampling techniques. However, PCE construction as
a linear regression model is a very problem specific task and it can be highly
beneficial to use methods that exploit information from the given mathematical
model and sequentially update the surrogate model – referred to as _active
learning_. Active learning is a common approach for surrogate-based
reliability analysis, wherein an initial experimental design is iteratively
updated based on the current estimate of the limit-state surface [13, 14, 15].
Active learning for reliability analysis with PCE was used e.g. in [16, 17,
18]. For general UQ studies, some recent studies have focused on general
sequential sampling for PCE based on space-filling criteria or alphabetical
optimality [19, 20]. However, it is beneficial to use both _exploitation_
(leveraging model behavior) criteria and _exploration_ (space filling)
criteria to define an optimally balanced criterion [21]. Such sequential
sampling for sparse Bayesian learning PCE combining both aspects – epistemic
uncertainty of the statistical inference (exploration) together with quadratic
loss function (local exploitation) – was recently proposed in [22]. However,
its application is limited to PCE built by sparse Bayesian learning only.

The authors of this paper recently proposed a general active learning method
based on sequential adaptive variance-based sampling [1], which is an
efficient tool for accurate surrogate modeling that is sufficiently general
for further extension [23]. Although this approach leads to superior results
in comparison to standard approaches without active learning, it is limited by
the inherently smooth nature of the PCE. More specifically, polynomial basis
functions are not able to approximate functions with discontinuities or
singularities. Moreover, it is necessary to use high-order polynomials to
approximate functions with local non-linearities, even when the rest of the
input random space could be easily approximated by a low-order PCE. This can
lead to spurious oscillations in the approximation and over-fitting. To
overcome this limitation, we propose a method to construct localized PCEs
based on the concept of _divide-and-conquer_ , i.e. decomposition of the input
random space to sub-domains approximated by many low-order PCEs instead of a
single high-order global PCE. Although this concept is not entirely new in
stochastic finite elements [24] and stochastic collocation [25, 26], there is
no such approach for non-intrusive PCE. However there are two primary
techniques based on similar concepts as described in the following section.

###  1.1 Related Developments

Stochastic Spectral Embedding (SSE) [2] is a general approximation technique
based on a decomposition of the input random space and the construction of
embedded local approximations. Although it is generally possible to use any
spectral approximation technique, it is beneficially coupled with PCE. SSE is
based on a novel idea of _embedding_ – instead of constructing local
approximations of the original mathematical model, local surrogates are
constructed to approximate the _residuals_ between the model and approximation
from the previous level of the decomposed space. Although such an approach can
lead to significant improvement in comparison to a single global approximation
[2], it is not a sequential approach based on active learning and thus it does
not iteratively reflect new information obtained from the previous steps of
the algorithm. Active learning is crucial in analysis of functions with
discontinuity or singularity because it allows for the aforementioned
exploration and exploitation necessary to find and resolve these features. For
the sake of completeness, active learning for SSE has been proposed for
reliability analysis [27], but it does not lead to an accurate approximation
over the entire input random space. Its accuracy is limited to regions around
the limit surface, which are important for an estimation of failure
probability.

The second related technique is Multi-element generalized Polynomial Chaos
Expansion (ME-gPC) [28]. ME-gPC was developed as an extension of generalized
PCE based on Wiener-Askey scheme [29] allowing analysis of models with
arbitrary distribution of input random vector. The ME-gPC method consists of
three main parts: decomposition of the input random space, numerical
construction of locally orthogonal polynomials and an adaptive procedure based
on the decay rate of local error in estimated variance derived from local PCE.
ME-PCE applies an $h$-type mesh refinement procedure akin to mesh refinement
in finite element methods. By doing so, they introduce a structured grid of
uniform points in each new element and solve for the PCE coefficients. This
can be cumbersome and does not afford the flexibility to adaptively select
sparse and near-optimal training points. Moreover, we note that the ME-gPC was
created mainly for uncertainty propagation in models with arbitrary input
distributions, and thus in contrast to SSE, its objective is not necessarily
to construct the best possible surrogate model using adaptive algorithms, but
rather to minimize errors in response statistics. This is a subtle, but
important difference that distinguishes its use as a predictive tool from that
of a tool for statistical estimation.

###  1.2 Contributions of this paper

This paper describes a novel method, termed Domain Adaptive Localized PCE
(DAL-PCE) that applies adaptive sequential decomposition of the input random
space and adaptive sequential sampling within the sub-domains. Both of these
features are based on recently a proposed criterion for variance-based
sequential statistical sampling, developed specifically for PCE in [30]. In
the context of previously described methods SSE and ME-gPC, the proposed novel
approach can be though to lie between them. Like SSE, it is developed
specifically for the construction of accurate surrogate models, especially for
functions with high non-linearity or discontinuity. But the decomposition of
the input random space is rather similar to ME-gPC. The uniqueness of our
proposal lies in the combination of active learning, sequential sampling,
sequential decomposition of the input space and regression-based PCE using
sparse solvers such as Least Angle Regression (LARS) allowing adaptivity and
learning in each iteration of the proposed algorithm.

##  2 Polynomial Chaos Expansion

Assume a probability space ($\Omega,\pazocal{F},\pazocal{P}$), where $\Omega$
is an event space, $\pazocal{F}$ is a $\sigma$-algebra on $\Omega$ and
$\pazocal{P}$ is a probability measure on $\pazocal{F}$. If the input variable
of a mathematical model, $Y=f(X)$, is a random variable
$X(\omega),\omega\in\Omega$, the model response $Y$($\omega$) is also a random
variable. Assuming that $Y$ has a finite variance, PCE represents the output
variable $Y$ as a function of an another random variable $\xi$ called the
_germ_ with a known distribution

$Y=f(X)=f^{\textsf{PCE}}(\xi),$ (1)

and represents the function $f(X)$ via infinite polynomial expansion. A set of
polynomials, orthogonal with respect to the distribution of the germ, are used
as a basis of the Hilbert space $L^{2}$ ($\Omega,\pazocal{F},\pazocal{P}$) of
all real-valued random variables of finite variance, where $\pazocal{P}$ takes
over the meaning of the probability distribution. The orthogonality condition
is given by the inner product of $L^{2}$ ($\Omega,\pazocal{F},\pazocal{P}$)
defined for any two functions $\psi_{j}$ and $\psi_{k}$ for all $j\neq k$ with
respect to the weight function $p_{\xi}$ (probability density function of
$\xi$) as:

$\langle\psi_{j},\psi_{k}\rangle=\int\psi_{j}(\xi)\psi_{k}(\xi)p_{\xi}(\xi)\;\mathrm{d}\xi=0.$
(2)

This means that there are specific orthogonal polynomials associated with the
corresponding distribution of the germ via its weighting function. For
example, Hermite polynomials orthogonal to the Gaussian measure are associated
with normally distributed germs. Orthogonal polynomials corresponding to other
distributions can be chosen according to Wiener-Askey scheme [29] or
constructed numerically [31]. For further processing, it is beneficial to use
normalized polynomials (orthonormal), where the inner product of $i$th and
$j$th polynomials is equal to the Kronecker delta $\delta_{jk}$, i.e.
$\delta_{jk}=1$ if and only if $j=k$, and $\delta_{jk}=0$ otherwise.

In the case of $\boldsymbol{X}$ and ${\bm{\xi}}$ being vectors containing $M$
independent random variables, the polynomial $\Psi({\bm{\xi}})$ is
multivariate and it is built up as a tensor product of univariate orthonormal
polynomials, i.e.

$\Psi_{{\boldsymbol{\alpha}}}({\bm{\xi}})=\prod_{i=1}^{M}\psi_{\alpha_{i}}(\xi_{i}),$
(3)

where ${{\boldsymbol{\alpha}}}\in\mathbb{N}^{M}$ is a set of integers called
the _multi-index_ reflecting polynomial degrees associated to each $\xi_{i}$.
The quantity of interest (QoI), i.e. the response of the mathematical model
$Y=f(\boldsymbol{X})$, can then be represented as [5]

$Y=f(\boldsymbol{X})=\sum_{{\boldsymbol{\alpha}}\in\mathbb{N}^{M}}\beta_{{\boldsymbol{\alpha}}}\Psi_{{\boldsymbol{\alpha}}}({\bm{\xi}}),$
(4)

where $\beta_{{\boldsymbol{\alpha}}}$ are deterministic coefficients and
$\Psi_{{\boldsymbol{\alpha}}}$ are multivariate orthonormal polynomials.

###  2.1 Non-intrusive computation of PCE coefficients

For practical computation, the PCE expressed in Eq. (4) must be truncated to a
finite number of terms $P$. One can generally choose any truncation rule (e.g.
tensor product of polynomials up to the selected order $p$), but the most
common truncation is achieved by retaining only terms whose total degree
$|{\boldsymbol{\alpha}}|$ is less than or equal to a given $p$, in which case
the truncated set of PCE terms is then defined as

$\pazocal
A^{M,p}=\left\\{{{\boldsymbol{\alpha}}}\in\mathbb{N}^{M}:\left|{{\boldsymbol{\alpha}}}\right|=\sum_{i=1}^{M}\alpha_{i}\leq
p\right\\}.$ (5)

The cardinality of the truncated _index set_ $\pazocal A^{M,p}$ is given by

$\mathrm{card}\>\pazocal A^{M,p}=\frac{\left(M+p\right)!}{M!\>p!}\equiv P\,.$
(6)

When the PCE is truncated to a finite number of terms, there is an error
$\varepsilon$ in the approximation such that

$Y=\displaystyle f{(\boldsymbol{X})}=\sum_{{\boldsymbol{\alpha}}\in\pazocal
A}\beta_{{\boldsymbol{\alpha}}}\Psi_{{\boldsymbol{\alpha}}}({\bm{\xi}})+\varepsilon\,.$

From a statistical point of view, PCE is a simple linear regression model with
intercept. Therefore, it is possible to use _ordinary least squares_ (OLS)
regression to minimize the error $\varepsilon$.

Knowledge of vector ${\boldsymbol{\beta}}$ fully characterizes the
approximation via PCE. To solve for ${\boldsymbol{\beta}}$, first it is
necessary to create ${N_{\mathrm{sim}}}$ realizations of the input random
vector $\boldsymbol{X}$ and the corresponding results of the original
mathematical model $\pazocal Y$, together called the experimental design (ED).
Then, the vector of $P$ deterministic coefficients ${\boldsymbol{\beta}}$ can
be determined by OLS as

${\boldsymbol{\beta}}=({\bm{\Psi}}^{T}{\bm{\Psi}})^{-1}\
{\bm{\Psi}}^{T}\pazocal Y,$ (7)

where ${\bm{\Psi}}$ is the data matrix

${\bm{\Psi}}=\left\\{\Psi_{ij}=\Psi_{j}({\bm{\xi}}^{(i)}),\;i=1,\ldots,{N_{\mathrm{sim}}},\;j=0,\ldots,P-1\right\\}.$
(8)

A well-known problem, the _curse of dimensionality_ , states that $P$ is
highly dependent on the number of input random variables $M$ and the maximum
total degree of polynomials $p$, which is clear from Eq. (6). Considering that
estimation of ${\boldsymbol{\beta}}$ by regression requires at least
$\mathcal{O}(P\,\ln(P))$ number of samples for stable solution [10, 11], the
problem can become computationally highly demanding in case of a large or
strongly non-linear stochastic models. Although one can use advanced model
selection algorithms such as Least Angle Regression (LAR) [32, 4], orthogonal
matching pursuit [33] or Bayesian compressive sensing [34] to find an optimal
set of PCE terms, and thus reduce the number of samples needed to compute the
unknown coefficients, the benefit of these techniques is significant only if
the true coefficient vector is sparse or compressible. The sparse set of basis
functions obtained by any adaptive algorithm is further denoted by
$\pazocal{A}$ for the sake of clarity.

###  2.2 Approximation Error Estimation

Once the PCE is constructed, it is crucial to estimate its accuracy. Further,
the PCE accuracy can be used to directly compare several PCEs to choose the
best surrogate model. Ideally the ED should be divided into validation and
training sets, but this might be extremely computationally demanding in
engineering applications with complex numerical models. Therefore in the field
of uncertainty quantification (UQ) of engineering models, it is preferred to
estimate the approximation error directly from the training set, without any
additional sampling of the original model. A common choice is the coefficient
of determination $R^{2}$, which is well-known from machine learning or
statistics. However, $R^{2}$ may lead to over-fitting and thus advanced
methods should be used. One of the most widely-used methods is the leave-one-
out cross-validation (LOO-CV) error $Q^{2}$. The LOO-CV is based on residuals
between the original surrogate model and the surrogate model built with the ED
while excluding one realization. This approach is repeated for all
realizations in the ED and the average error is estimated. Although the
calculation of $Q^{2}$ is typically highly time-consuming, it is possible to
obtain results analytically from a single PCE as follows [35]:

$Q^{2}=\frac{\displaystyle{\frac{1}{{N_{\mathrm{sim}}}}\sum_{i=1}^{{N_{\mathrm{sim}}}}\displaystyle{{{\left[{\frac{g\left(\mbox{\boldmath$x$}^{(i)}\right)-g^{\textsf{PCE}}\left(\mbox{\boldmath$x$}^{(i)}\right)}{1-{h_{i}}}}\right]}^{2}}}}}{\sigma^{2}_{Y,\mathrm{ED}}},$
(9)

where ${\sigma^{2}_{Y,\mathrm{ED}}}$ is the variance of the ED calculated
using the original mathematical model and ${h_{i}}$ represents the $i$th
diagonal term of matrix
$\mathbf{H}={\bm{\Psi}}\left({\bm{\Psi}}^{T}{\bm{\Psi}}\right)^{-1}{\bm{\Psi}}^{T}$.

###  2.3 Statistical Moments Derived from PCE

The form of PCE as a linear summation over orthonormal polynomials allows for
powerful and efficient post-processing. In particular, once a PCE
approximation is created, it is possible to directly estimate statistical
moments of the output from the expansion.

The first statistical moment (the mean value) is simply the first
deterministic coefficient of the expansion
$\mu_{Y}=\big{<}Y^{1}\big{>}=\beta_{{\boldsymbol{0}}}$. The second raw
statistical moment, $\big{<}Y^{2}\big{>}$, can be estimated by

$\displaystyle\left\langle{{Y^{2}}}\right\rangle$
$\displaystyle=\int{\left[{\sum\limits_{{\boldsymbol{\alpha}}\in\pazocal{A}}{{\beta_{{\boldsymbol{\alpha}}}}{\Psi_{{\boldsymbol{\alpha}}}}\left({\bm{\xi}}\right)}}\right]}^{2}p_{{\bm{\xi}}}\left({\bm{\xi}}\right)\;\mathrm{d}{\bm{\xi}}=\sum\limits_{{\boldsymbol{\alpha}}_{1}\in\pazocal{A}}\sum\limits_{{\boldsymbol{\alpha}}_{2}\in\pazocal{A}}\beta_{{{\boldsymbol{\alpha}}}_{1}}\beta_{{{\boldsymbol{\alpha}}}_{2}}\int\Psi_{{{{\boldsymbol{\alpha}}}_{1}}}\left({\bm{\xi}}\right)\Psi_{{{{\boldsymbol{\alpha}}}_{2}}}\left({\bm{\xi}}\right)p_{\bm{\xi}}\left({\bm{\xi}}\right)\;\mathrm{d}{\bm{\xi}}$
(10)
$\displaystyle=\sum\limits_{{{\boldsymbol{\alpha}}}\in{\pazocal{A}}}{\beta_{{\boldsymbol{\alpha}}}^{2}}{\int{{\Psi_{{\boldsymbol{\alpha}}}}\left({\bm{\xi}}\right)}^{2}}{p_{\bm{\xi}}}\left({\bm{\xi}}\right)\;\mathrm{d}{\bm{\xi}}=\sum\limits_{{{\boldsymbol{\alpha}}}\in{\pazocal{A}}}{\beta_{{\boldsymbol{\alpha}}}^{2}}\left\langle{{\Psi_{{\boldsymbol{\alpha}}}},{\Psi_{{\boldsymbol{\alpha}}}}}\right\rangle.$

Considering the orthonormality of the polynomials, it is possible to obtain
the variance $\sigma_{Y}^{2}=\big{<}Y^{2}\big{>}-\mu_{Y}^{2}$ as the sum of
all squared deterministic coefficients except the intercept (which represents
the mean value), i.e.

$\sigma_{Y}^{2}=\sum_{\begin{subarray}{c}{\boldsymbol{\alpha}}\in\pazocal{A}\\\
{\boldsymbol{\alpha}}\neq{\boldsymbol{0}}\end{subarray}}\beta_{{\boldsymbol{\alpha}}}^{2}.$
(11)

Note that the computation of higher statistical central moments, specifically
skewness $\gamma_{Y}$ ($3^{\text{rd}}$ moment) and kurtosis $\kappa_{Y}$
($4^{\text{th}}$ moment), are more complicated since they require triple and
quad products. These can be obtained analytically only for certain polynomial
families, e.g. formulas for Hermite and Legendre polynomials (and their
combination) can be found in [30].

##  3 Active Learning-based Domain Adaptive Localized PCE (DAL-PCE)

In this section, we propose a novel methodology to constructed localized PCEs
designed for highly non-linear functions, termed Domain Adaptive Localized PCE
(DAL-PCE). Instead of increasing the maximum polynomial order $p$
($p$-adaptivity), which brings high computational requirements due to the
_curse of dimensionality_ , we propose to decompose the input random space
into several sub-domains approximated by low-order PCEs ($h$-adaptivity).
Although this idea is not entirely new, we use this approach in combination
with novel active learning methods to identify domains for refinement and for
sequential sample selection and regression-based PCEs. This allows us to use
any sparse adaptive solver (e.g. LAR) and thus it can be easily implemented
into the existing software packages [36, 37]. In the following sections, we
define the requisite components of the proposed method and provide an
algorithm (Algorithm 1) for its implementation.

###  3.1 Variance-based Adaptive Sequential Sampling

The decomposition of the input random space is a sequential process coupled
with adaptive sampling assuring optimal coverage of the sub-domains of
interest. The whole process thus consists of two steps: (i) identification of
an important sub-domain, that is, a domain that is either large compared to
other sub-domains or that is associated with a high local variance; and (ii)
identification of the best positions for additional samples extending the
current ED in the selected sub-domain. Each of these steps must be based on a
criterion that balances _exploration_ of the input random space with
_exploitation_ of the surrogate model, which in our case is in the form of a
PCE. The $\Theta$-criterion for adaptive sequential sampling, which is driven
by the output variance and its approximation via local variance using PCE [1],
is employed for both steps. We will first discuss the process for adaptive
sequential sampling within a specified sub-domain in this section. This will
be followed by the process for refinement of the domain in the subsequent
sections.

Consider a pool of candidate samples containing realizations of the random
vector ${\bm{\xi}}$ generated by an arbitrary sampling technique, e.g., Latin
Hypercube Sampling (LHS) [38, 39] or Coherence sampling [40, 41, 10]. From
this pool of candidates, we select the best sample using a method inspired by
the sequential sampling proposed in [21] and based on Koksma-Hlawka inequality
[42]. The $\Theta$-criterion for PCE, which accounts for both variation of the
function and discrepancy of the samples, was proposed as follows [1]:

$\displaystyle\Theta({\bm{\xi}}^{(\mathrm{c})})\equiv\Theta^{c}$
$\displaystyle=\underbracket{\sqrt{\sigma_{\\!\\!\pazocal{A}}^{2}({\bm{\xi}}^{(\mathrm{c})})\cdot\sigma_{\\!\\!\pazocal{A}}^{2}({\bm{\xi}}^{(\mathrm{s})})}}_{\mathrm{ave\;variance\;density}}\>\underbracket{l_{\mathrm{c},\mathrm{s}}^{M}}_{\mathrm{vol.}}\equiv\sqrt{\sigma_{\mathrm{c}}^{2}\cdot\sigma_{\mathrm{s}}^{2}}\>l_{\mathrm{c},\mathrm{s}}^{M}.$
(12)

The criterion is a product of two terms – the _exploitation_ term (denoted as
“ave variance density”) and the _exploration_ part (the distance term
$l_{\mathrm{c},\mathrm{s}}$ raised to the domain dimension) – which are
multiplied to maintain an optimal balance between exploration and exploitation
[1].

The _exploration_ aspect is maintained by accounting for the distance
$l_{\mathrm{c},\mathrm{s}}$ between a candidate ${\bm{\xi}}^{(\mathrm{c})}$
and its nearest neighboring realization from the existing ED,
${\bm{\xi}}^{(\mathrm{s})}$ as

$l_{\mathrm{c},\mathrm{s}}=\sqrt{\sum_{i=1}^{M}|\xi_{i}^{(\mathrm{c})}-\xi_{i}^{(\mathrm{s})}|^{2}}.$
(13)

If the criterion was reduced to this term only, sequential filling of the
greatest empty regions would occur, converging to uniform space coverage in
the spirit of the space-filling “miniMax criterion” [43, 44, 45].

The _exploitation_ component is motivated by the desire to sample points in
regions with the greatest contributions to the total variance of the QoI
$\sigma_{Y}^{2}$, i.e. at points with the highest _variance density_. Once the
PCE has been established at any given stage of the algorithm, the _variance
density_ is computationally cheap to evaluate for any location ${\bm{\xi}}$ as

$\sigma_{\\!\\!\pazocal{A}}^{2}({\bm{\xi}})=\big{[}\sum_{\begin{subarray}{c}{\boldsymbol{\alpha}}\in\pazocal{A}\\\
{\boldsymbol{\alpha}}\neq{\boldsymbol{0}}\end{subarray}}\beta_{{\boldsymbol{\alpha}}}{\Psi_{{\boldsymbol{\alpha}}}\left({\bm{\xi}}\right)}\big{]}^{2}{p_{\xi}}\left({\bm{\xi}}\right).$
(14)

The local variance is therefore estimated directly using the basis functions
and coefficients $\beta$ of the PCE. When considering a candidate
“$\mathrm{c}$”, an estimate of the variance contribution of the region between
the candidate and its nearest neighbor “$\mathrm{s}$” may be obtained by
averaging the local variance densities between the two. Therefore, we can say
that the candidate with the greatest $\Theta^{c}$ criterion is the one that
represents the largest amount of total variance to be refined by its
selection.

A significant advantage of this method is the ability to add candidates into
an existing ED one-by-one. Thus, it can be employed at any moment of the PCE
construction process. Moreover, this learning function can be combined with
any sampling algorithm for the construction of the initial ED and candidates
for extension. The ideas behind the $\Theta$ criterion will now be used in the
proposed domain decomposition and ED extension algorithm.

###  3.2 Decomposition of Input Random Space

The core of the proposed approach is a sequential decomposition of the input
random space $\mathcal{D}$ for the construction of local approximations. This
approach assumes that the original mathematical model can be approximated by
piecewise low-order PCEs that are valid only in individual sub-domains of
$\mathcal{D}$. Therefore, in the proposed approach, the input random space is
sequentially decomposed into $n_{\mathcal{D}}$ smaller non-overlapping sub-
domains $\mathcal{D}_{i}\subset\mathcal{D}$ that collectively fill the full
input random space $\mathcal{D}$, i.e.

$\bigcup_{i=1}^{n_{\mathcal{D}}}\mathcal{D}_{i}=\mathcal{D}\quad\text{such
that}\quad\mathcal{D}_{i}\cap\mathcal{D}_{j}=\emptyset\quad\forall i,j$ (15)

In each iteration of the algorithm, a single sub-domain $\bm{\mathcal{D}_{i}}$
(referred to as the parent) is identified for refinement and divided by a
plane perpendicular to the direction of one selected input random variable.
Specifically, $\bm{\mathcal{D}_{i}}$ is divided into a refinement-child
$\mathcal{D}_{i}$, which is further processed, and an inheriting-child
$\mathcal{D}_{i}^{\star}$ adopting the PCE from the parent as illustrated for
a one-dimensional function in Fig. 1. In this case, we see that the space is
divided into two subdomains. In the left (refinement child) a new PCE is
constructed. In the right (inheriting child), the original PCE is retained.
Such process assures an exhaustive decomposition into disjoint subsets i.e.
$\bm{\mathcal{D}_{i}}=\mathcal{D}_{i}\oplus\mathcal{D}_{i}^{\star}$. This
sequential domain decomposition is illustrated in Fig. 2, which depicts the
original input random space and the first four iterations of the decomposition
process.

Figure 1: The first iteration of the algorithm: the original sub-domain is
split and the new local PCE is constructed in $\mathcal{D}_{i}$ (red
background), while the second part in $\mathcal{D}_{i}^{\star}$ inherits the
PCE approximation from the original domain.

In contrast to SSE [2], the selection of a single sub-domain for refinement in
each iteration is based on an active learning approach, the details of which
are provided in subsequent sections. Importantly, actively integrating
information from the original mathematical model leads to a significantly more
effective decomposition of the space and thus assures accurate approximations,
even for small-size EDs. On the other hand, the identified decomposition and
the associated ED are directly connected to the given mathematical model and
therefore might be inefficient for general statistical analysis.

The complete surrogate model is assembled from the $n_{\mathcal{D}}$ local
PCEs associated with all sub-domains $\mathcal{D}_{i}$ as:

$Y\approx\sum_{i=0}^{n_{\mathcal{D}}}\sum_{{\boldsymbol{\alpha}}_{i}\in\pazocal
A_{i}}\beta_{{\boldsymbol{\alpha}}_{i}}\Psi_{{\boldsymbol{\alpha}}_{i}}({\bm{\xi}})\mathbb{1}_{\mathcal{D}_{i}}({\bm{\xi}}),$
(16)

where $\mathbb{1}_{\mathcal{D}_{i}}({\bm{\xi}})$ represents indicator
function, i.e. $\mathbb{1}_{\mathcal{D}_{i}}({\bm{\xi}})=1$ only if
${\bm{\xi}}\in\mathcal{D}_{i}$ and
$\mathbb{1}_{\mathcal{D}_{i}}({\bm{\xi}})=0$ otherwise. In other words, to
approximate the original model at any point, it suffices to determine the one
relevant sub-domain and use the corresponding local PCE. Each such local PCE
has its own set of basis functions $\pazocal A_{i}$ and corresponding
coefficients $\beta_{{\boldsymbol{\alpha}}_{i}}$, which can be obtained by any
model-selection algorithm. In this paper the OLS and LAR algorithms are
employed, but generally any non-intrusive technique can be used.

Figure 2: The first four steps of the decomposition of a 3D space of input
random variables. The thick black lines outline the parent domain selected for
division. The red and green boxes inside it represent the two newly created
refinement-child $\mathcal{D}_{i}$ (red) and inheriting-child
$\mathcal{D}_{i}^{\star}$ (green) sub-domains created by splitting the parent
domain $\bm{\mathcal{D}_{i}}$ (bold boundaries), selected via Eq. (17), by the
cutting plane (blue). The cutting plane is perpendicular to the variable
selected for splitting (blue arrow).

###  3.3 Domain Selection via Modified Variance-based Criterion

The selection process to identify the “best” subdomain for possible division
is governed by extending the $\Theta$-criterion from Eq. (12) as follows:

$\Theta_{i}=\underbracket{\pazocal{W}_{i}\cdot\exp{(Q_{i}^{2})}}_{\mathrm{weight\>of\>subdomain}}\cdot\underbracket{\sqrt{\sigma_{\\!\\!\pazocal{A}_{i}}^{2}(\bm{\xi}^{(c)})\cdot\sigma_{\\!\\!\pazocal{A}_{i}}^{2}(\bm{\xi}^{(s)})}\>l_{c,s}^{M}}_{\Theta^{c}\>\text{in
\>}i\text{th subdomain}}.$ (17)

This extended criterion aims to identify sub-domains of the input random space
associated with the maximum value of $\Theta^{c}$, while simultaneously
accounting for the size of each subdomain and the accuracy of the existing
local PCE. The former is calculated using Eq. (12) calculated for a rich pool
of screening global candidates, while the latter are measured by incorporating
the volume of each sub-domain $\pazocal{W}_{i}$ and the LOO-CV error
$Q_{i}^{2}$, respectively. The LOO-CV term, $\exp{(Q_{i}^{2})}$, can be
thought to artificially inflate the domain volume as a penalization for
inaccurate approximation. When the approximation is perfect $(Q_{i}^{2}=0)$
the true volume of the sub-domain is used. Meanwhile, a poor approximation
with $Q_{i}^{2}=1$ leads to roughly 2.72 times increased volume.

The three terms featured in Eq. (17) aim at different aspects affecting the
accuracy of the final surrogate model: large sub-domains are preferred by
$\pazocal{W}_{i}$, sub-domains containing poor PCE approximation are promoted
via $\exp{(Q_{i}^{2})}$ and finally, $\Theta^{c}$ prefers sub-domains with
high concentration of variance. Note that $\Theta^{c}$ is calculated for a
rich pool of screening candidates, and $\pazocal{W}_{i}$ and
$\exp{(Q_{i}^{2})}$ are calculated directly from the geometry of existing sub-
domain and the local PCE model, respectively. The product of all three terms
in the extended criterion therefore maintains the desired balance and assures
the selection of the sub-domain, $\mathcal{D}_{i}$, that currently seems to be
the most important for increasing the accuracy of the PCE surrogate model.

Sub-domain $\mathcal{D}$ with the greatest $\Theta_{i}$ is selected and one of
the operations described in detail in Sec. 3.6 is performed, depending on
whether $\mathcal{D}_{i}$ contains a critical number of ED points. Two
scenarios can occur:

  * •

$\mathcal{D}_{i}$ contains a sufficient number of ED points
($n_{i}\geq{n_{\mathrm{sim}}}$) to ensure accuracy of a PCE on the domain.
Therefore, it becomes a parent $\bm{\mathcal{D}_{i}}$ (bold boundaries in Fig.
2) and is divided into two parts by a selected rule. The child domain
containing the decisive candidate with the greatest $\Theta^{c}$ becomes the
refinement-child $\mathcal{D}_{i}$ (see the red subdomains in steps $1-4$ in
Fig. 2). The remaining volume becomes an inheriting-child denoted
$\mathcal{D}_{i}^{\star}$ (see the green subdomains in Fig. 2), which retains
the PCE from the parent. Division occurs by a cutting plane, oriented
perpendicular to the selected direction (blue arrows in Fig. 2) and naturally,
the coordinates of the cutting plane are restricted to the bounding box of the
selected parent $\bm{\mathcal{D}_{i}}$, see Sec. 3.6. If needed, the
refinement-child domain $\mathcal{D}_{i}$ is sequentially filled with
additional ED points (according to $\Theta^{c}$) to reach
$n_{i}={n_{\mathrm{sim}}}$ needed to construct a new PCE approximation.

  * •

$\mathcal{D}_{i}$ does _not_ contain a sufficient number of ED points
($n_{i}<{n_{\mathrm{sim}}}$). The domain is not divided because the suggestion
for division is based on insufficient information. Instead, new ED points are
sequentially added to $\mathcal{D}_{i}$, again using the $\Theta^{c}$
criterion. Note that this scenario practically arises when the selected domain
was an inheriting-child in the previous iteration. In this case, the selected
domain has inherited a PCE model that was constructed over a larger domain.
When that domain was divided, it was left with an insufficient number of
points from which to construct a new PCE.

###  3.4 PCE Basis Functions

Without loss of generality, the proposed method operates on the
$M$-dimensional unit hypercube with uniform distributions of input random
variables, i.e. $\boldsymbol{X}\sim\pazocal{U}[0,1]^{M}$. In the case of a
general joint probability distribution of $\boldsymbol{X}$, it is always
possible to transform input random vector to the unit hypercube by Rosenblatt
transformation [46], Nataf transformation [47] or various methods based on
copulas [48]. Standard normalized Legendre polynomials, orthonormal to the
uniform distribution, can thus be used as basis functions for the PCE.
However, due to the decomposition of the input random space to smaller sub-
domains, each with lower bound $a_{i}$ and upper bound $b_{i}$, it is
necessary to use univariate scaled orthonormal Legendre polynomials of $n$th
order $\tilde{\psi}_{n}(\xi)$ defined as follows:

$\tilde{\psi}_{n}(\xi)=\psi_{n}\left(\frac{2\xi-
a_{i}-b_{i}}{b_{i}-a_{i}}\right),$ (18)

where $\psi_{n}$ represents standard orthonormal Legendre polynomials.
Naturally, the transformation of the original input random vector to the unit
hypercube might bring additional non-linearity, and thus one might prefer the
direct construction of polynomials locally orthonormal to the given original
probability measure as proposed in the Me-gPC [28]. While certainly possible,
this brings additional computational demands and thus it is not employed here.

###  3.5 Local and Global Statistical Estimates from DAL-PCE

The significant advantage of PCE is that analytically post-processing of the
expansion yields highly efficient estimates of statistical moments [30],
sensitivity indices [8] and LOO-CV [4]. In the proposed DAL-PCE, since the
original domain $\mathcal{D}$ is decomposed into a set of sub-domains (see Eq.
(15)), standard analytical post-processing can be applied locally and global
characteristics can be obtained by simple weighted summations that converge to
the true values as $n_{\mathcal{D}}$ increases. Specifically, the global mean
value and variance of a QoI are obtained from localized PCEs (denoted by
subscript $\mathcal{D}_{i}$) as follows:

$\mu_{Y}=\sum_{i=1}^{n_{\mathcal{D}}}\pazocal{W}_{i}\beta_{0_{i}}=\sum_{i=1}^{n_{\mathcal{D}}}\pazocal{W}_{i}\mu_{\mathcal{D}_{i}},\
$ (19)
$\sigma_{Y}^{2}=\sum_{i=1}^{n_{\mathcal{D}}}\pazocal{W}_{i}\sum_{\begin{subarray}{c}{\boldsymbol{\alpha}}_{i}\in\pazocal{A}_{i}\\\
{\boldsymbol{\alpha}}_{i}\neq{\boldsymbol{0}}\end{subarray}}\beta_{{\boldsymbol{\alpha}}_{i}}^{2}=\sum_{i=1}^{n_{\mathcal{D}}}\pazocal{W}_{i}\sigma_{\mathcal{D}_{i}}^{2}.$
(20)

where the local mean $\mu_{\mathcal{D}_{i}}$ and variance
$\sigma_{\mathcal{D}_{i}}^{2}$ are obtained as described in Section 2.3.

Local Sobol’ indices, $S_{\mathcal{D}_{i}}$, of any order can be derived
directly from localized PCEs and their first-order (main effect) estimates are
given by

$S_{\mathcal{D}_{i}}^{X_{j}}=\frac{1}{{\sigma_{\mathcal{D}_{i}}^{2}}}{{\sum\limits_{{\boldsymbol{\alpha}}_{i}\in\pazocal{A}^{X_{j}}_{i}}{\beta_{\bf{{\boldsymbol{\alpha}}_{i}}}^{2}}}}{}\;\quad{\pazocal{A}^{X_{j}}_{i}}=\left\\{{{{\boldsymbol{\alpha}}_{i}}\in{\pazocal{A}_{i}}:{\alpha_{i}^{j}}>0,{\alpha_{i}^{k\neq
j}}=0}\right\\}.$ (21)

These local Sobol’ indices are used in the DAL-PCE to determine the cut
direction (see Section 3.6). Likewise, global Sobol’ indices can be obtained
easily from weighted summation of local contributions to partial variances
normalized by $\sigma_{Y}^{2}$ as follows:

$S_{X_{j}}=\frac{\sum_{i=1}^{n_{\mathcal{D}}}\pazocal{W}_{i}{{\sum\limits_{{\boldsymbol{\alpha}}_{i}\in\pazocal{A}^{X_{j}}_{i}}{\beta_{{{\boldsymbol{\alpha}}_{i}}}^{2}}}}}{\sigma_{Y}^{2}}.$
(22)

Similarly, global LOO-CV, $Q^{2}$, of a QoI can be approximated by the
weighted summation of the local contributions as

$Q^{2}=\sum_{i=1}^{n_{\mathcal{D}}}\pazocal{W}_{i}Q^{2}_{\mathcal{D}_{i}},$
(23)

where $Q^{2}_{\mathcal{D}_{i}}$ are obtained from each local PCE using Eq.
(9).

These estimates are used throughout the proposed DAL-PCE, as described in
detail next.

###  3.6 Numerical Algorithm

Based on the presented theoretical background, we now present the numerical
algorithm for the domain adaptive localized PCE. As mentioned above, the whole
process can be divided to two iterative tasks: (i) decomposition of the input
random space and (ii) construction of localized PCEs. Both of these tasks are
described in the following paragraphs with specific reference to the steps in
Algorithm 1.

Algorithm 1 DAL-PCE: Active Domain Decomposition and Construction of Localized
PCEs

1:Input: maximum local polynomial order $p$, number of screening global
candidates $n_{c,g}$, number of local candidates $n_{c,l}$, number of
iterations $n_{\mathrm{iter}}$

2:set the minimum number of realizations for local PCE construction
${n_{\mathrm{sim}}}\in\langle P,2P\rangle$

3:generate a rich pool of $n_{c,g}$ screening candidates

4:generate the initial ED (size ${n_{\mathrm{sim}}}$) and construct the
initial global PCE

5:for $1$ to $n_{\mathrm{iter}}$ do

6: identify the sub-domain $\mathcal{D}_{i}$ with the highest $\Theta_{i}$
based on screening candidates

7: $n_{i}\leftarrow$ number of ED samples existing in $\mathcal{D}_{i}$

8: if $n_{i}\geq{n_{\mathrm{sim}}}$ then

9: the identified sub-domain $\mathcal{D}_{i}$ becomes a parent
$\bm{\mathcal{D}_{i}}$

10: identify the direction of the highest first-order Sobol’ index
$S_{\mathcal{D}_{i}}$ of the parent $\bm{\mathcal{D}_{i}}$

11: restrict coordinates of $\bm{\mathcal{D}_{i}}\rightarrow\mathcal{D}_{i}$
and create $\mathcal{D}_{i}^{\star}$

12: $n_{i}\leftarrow$ number of ED samples existing in $\mathcal{D}_{i}$

13: end if

14: generate $n_{c,l}$ local candidates in $\mathcal{D}_{i}$

15: while $n_{i}<{n_{\mathrm{sim}}}$ do

16: extend size of local ED $n_{i}$ using the local $\Theta^{c}$ criterion

17: end while

18: reconstruct local PCEs in the $\mathcal{D}_{i}$

19:end for

20:Output: list of subdomains and corresponding PCEs

The first task identifies the important sub-domain $\mathcal{D}_{i}$ that
should be divided and over which low-order local PCE should be constructed.
The sub-domain $\mathcal{D}_{i}$ is specifically identified using the
$\Theta_{i}$ criterion from Eq. (17), which again incorporates three important
characteristics for accurate surrogate modeling – the size of the sub-domain
$\pazocal{W}_{i}$, the accuracy of the existing local PCE measured by
$Q^{2}_{\mathcal{D}_{i}}$, and the original $\Theta^{c}$ criterion measuring
the variance contribution in $\mathcal{D}_{i}$. While $\pazocal{W}_{i}$ and
$Q^{2}_{\mathcal{D}_{i}}$ are computed for the whole sub-domain, $\Theta^{c}$
is computed at specific realizations of input random vector. Therefore, it is
necessary to cover the sub-domains by a sufficiently large number of screening
candidates, such that the total global number of screening candidates is given
by $n_{c,g}$. Based on numerical experiments, we recommend $n_{c,g}\geq
1000\,M$ to ensure that each sub-domain contains a sufficient number of
screening candidates. Note that the screening candidates are used only to
identify $\mathcal{D}_{i}$ [_step 5_]. They are not used for the ED, and thus
even high $n_{c,g}$ does not bring any additional computational demand.

Once $\mathcal{D}_{i}$ is identified, it is necessary to check whether there
are enough samples to construct a PCE inside the sub-domain. We start with
finding out how many points belong to the selected domain $\mathcal{D}_{i}$
[_step 6_]. If the number of samples in the identified sub-domain, $n_{i}$, is
greater than (or equal to) ${n_{\mathrm{sim}}}$ [_step 7_], a local PCE
already exists for $\mathcal{D}_{i}$. The subdomain is then assigned as a
parent $\bm{\mathcal{D}_{i}}$ for division [_step 8_] and the first-order
Sobol’ indices are estimated by Eq. (22) [_step 9_]. This identified parent
$\bm{\mathcal{D}_{i}}$ is divided in the direction of the highest first-order
Sobol’ index $S_{\mathcal{D}_{i}}^{X_{j}}$. The new restricted coordinates of
refinement-child $\mathcal{D}_{i}$ are identified and the inheriting-child
$\mathcal{D}_{i}^{\star}$ is created [_step 10_]. Further, the number of ED
samples $n_{i}$ in the refinement-child $\mathcal{D}_{i}$ is determined [_step
11_]. On the other hand, if the identified sub-domain $\mathcal{D}_{i}$ does
not contain enough samples (i.e. $n_{i}<{n_{\mathrm{sim}}}$), the inherited
PCE from the previous iteration is not sufficiently local (it was trained over
a domain that has since been divided) and it is necessary to add new samples
to $\mathcal{D}_{i}$ before constructing a new local PCE.

The second task of the proposed algorithm is sequential sampling and adaptive
PCE construction in sub-domain $\mathcal{D}_{i}$. Recall that this domain may
be either

  1. (i)

a refinement-child that was just divided but does not contain a sufficient
number of points ($n_{i}<{n_{\mathrm{sim}}}$) or,

  2. (ii)

an inheriting-child that now does not contain at least ${n_{\mathrm{sim}}}$ ED
samples.

Next, a set of local candidates is generated in region $\mathcal{D}_{i}$
[_step 13_]. To ensure sufficient assessment of the coverage of the domain,
the number of local candidates is empirically recommended as
$n_{c,l}\in\langle 3P,5P\rangle$ [1]. From these candidates, the standard
$\Theta^{c}$ criterion in Eq. (12) is used to iteratively select the best
candidates until there are ${n_{\mathrm{sim}}}$ samples in $\mathcal{D}_{i}$
[_step 14-16_]. This sequential extension of the sample in $\mathcal{D}_{i}$
is adaptive in the sense that the pairwise distances in Eq. (12) between
candidates and existing ED points are updated after the addition of each new
point. However, because $n_{i}<{n_{\mathrm{sim}}}$ the local variance
densities are estimated from the previously existing PCE, which cannot be
updated until a sufficient number of samples are available in
$\mathcal{D}_{i}$.

The last step of each iteration is to construct the local PCE using scaled
Legendre polynomials as basis functions (see Eq. (18)) [_step 17_]. Any non-
intrusive technique can be used to estimate the coefficients
${\boldsymbol{\beta}}$; we use LARS and OLS for an adaptive construction of
the local PCEs in this paper. At the end of the iteration, all sub-domains are
re-numbered and a list of sub-domains with corresponding PCEs can be exported
or the next iteration can be started.

###  3.7 Adaptivity in PCE Construction and Domain Decomposition

Adaptivity is central to the proposed DAL-PCE. In the proposed algorithm,
there are two types of adaptivity employed:

  1. (i)

adaptivity in PCE construction (selection of the optimal set of basis
functions), and

  2. (ii)

adaptivity in domain decomposition

Since the PCE can be constructed by any regression technique in each sub-
domain, PCE adaptivity is incorporated by sparse solvers and best model
selection algorithms, e.g. Least Angle Regression [32], orthogonal matching
pursuit [33] or Bayesian compressive sensing [34]. Although sparse solvers are
often used for PCE with high $p$, this adaptivity is also important for
reducing the number of basis functions (and thus the minimum number of ED
samples) for high-dimensional examples or, in our case, for very low-size ED
in each $\mathcal{D}_{i}$ approximated by low-$p$ local PCE.

The second type of adaptivity is the proposed adaptivity in the domain
decomposition. At any point in the iterative process, the existing ED samples
can be used to construct local PCEs or a single global PCE. The DAL-PCE is not
guaranteed to provide a better approximation than the global PCE. This can be
measured via $Q^{2}$, specifically by computing $Q^{2}_{\mathrm{local}}$ from
Eq. (23) and $Q^{2}_{\mathrm{global}}$ from a single global PCE according to
Eq. (9). If $Q^{2}_{\mathrm{local}}>Q^{2}_{\mathrm{global}}$ at a given
iteration, the domain decomposition is deemed to be poor and the whole
decomposition process is _re-started_. That is, the complete geometrical
decomposition is forgotten and all existing ED points are taken as an initial
ED for a brand new run of the algorithm. This is illustrated in Fig. 3 which
shows the decomposition (top) and the associated error (bottom) right before
the restart a) at ${N_{\mathrm{sim}}}=181$, b) the new decomposition and error
right after the restart, and c) the final decomposition/error which shows
significant improvement over the global PCE. These histories show the standard
$R^{2}$ error defined in Eq. (24). It is not necessary to check this criterion
at every iteration, but it is suggested to check it periodically, every
$n_{r}$ steps, to ensure adequate local refinement.

Figure 3: Illustration of domain decomposition restart. a) decomposition and
error evolution prior to restart, b) rebuilt decomposition and error drop
right after the restart, c) final decomposition and error showing that the
restart unlocks a dramatic decrease in approximation error.

###  3.8 Stopping Criteria

The proposed DAL-PCE algorithm can be fully automated by adding an adequate
stopping criterion. A simple but practical stopping criterion is based on
computational budget, i.e. once the total number of model evaluations
${N_{\mathrm{sim}}}$ or number of iterations $n_{\mathrm{iter}}$ have reached
a critical level/budget. One may also use a stopping criterion based on
decomposition pattern, e.g. the smallest or the largest volumes of any
subdomain, to ensure a desired resolution. Valuable stopping criterion can be
also obtained directly from $Q^{2}$, corresponding to a target/threshold level
of achieved accuracy. Regardless of the selected stopping criteria, it can
easily be applied before step 5 of the proposed algorithm (start of each
iteration).

##  4 Numerical Experiments

The proposed DAL-PCE is presented on four numerical examples of increasing
complexity and which illustrated different aspects of the approach. The
obtained results are compared (a) to the standard global PCE approach with
adaptive maximum order $p\in[5,25]$ and (b) to SSE [2], as current state-of-
the-art non-intrusive surrogate modeling technique based on the domain
decomposition. The PCE is constructed using the UQPy package [36] and the
original implementation of SSE is used from the UQLab package [37]. To compare
methods, the relative mean squared errors $\epsilon$ are calculated for all
three approximations $\tilde{f}$ on a validation set containing a large pool
of $10^{6}$ integration points generated by crude Monte Carlo according to:

$\epsilon(\boldsymbol{X})\coloneqq\frac{\mathbb{E}\Big{[}\big{(}f(\boldsymbol{X})-\tilde{f}(\boldsymbol{X})\big{)}^{2}\Big{]}}{\mathbb{D}\Big{[}f(\boldsymbol{X})\Big{]}},$
(24)

where $\mathbb{E}[]$ and $\mathbb{D}[]$ are the mean value and variance
operators, respectively.

To show representative results of the proposed DAL-PCE algorithm, the
calculations were repeated 100 times, and the same settings of the algorithm
for all examples were selected as follows: maximum local polynomial degree
$p=2$, number of global candidates $n_{c,g}=1000\ M$, number of local
candidates $n_{c,l}=5P$, minimum number of samples for local PCE construction
${n_{\mathrm{sim}}}=1.5P$, minimum number of iterations before checking for
restart $n_{r}=20$, and ${\boldsymbol{\beta}}$ are obtained by LARS and OLS
algorithm. Minimum number of samples in sub-domains required to justify an
expansions for SSE was set identically to DAL-PCE and polynomial order is
adaptively selected in the range $p\in[2,6]$. Since the SSE is not a
sequential approach, the presented results were obtained for 10 discrete
sample sets of increasing size to compare convergence of the method. Note that
all samples and candidates are generated by LHS for all compared approaches,
though it was shown [1] that for the variance-based sequential sampling, it is
significantly better to use advanced techniques such as Coherence D-optimal
sampling [41].

###  4.1 One-dimensional Toy Example

The first example involves a simple 1D function [2] that is extremely
difficult to approximate with PCE due to the third, highly nonlinear “$\exp$”
term:

$f(X)=-X+0.1\sin(30X)+\exp(-(50(X-0.65))^{2}),\quad X\sim\pazocal{U}[0,1].$
(25)

The poor performance of a single global PCE learned from 200 samples is
depicted by the blue line in Fig. 4c where it is clear that a single global
PCE is not able to accurately approximate the function even for a high number
of samples and high maximum polynomial order $p\in[5,25]$. This function was
originally developed to demonstrate the efficiency of SSE based on domain
decomposition and thus it is a natural choice for comparison of the proposed
DAL-PCE and SSE.

Fig. 4a-b show a typical realization of the DAL-PCE where the algorithm
sequentially decomposes the domain and adds additional samples to the ED.
Specifically shown are the 4th and 11th iterations. The boundaries of sub-
domains are represented by blue vertical lines and red dots show the positions
of samples in the ED. Once the algorithm discovers the highly nonlinear region
(the steep peak caused by $\exp$), it progressively refines this region and
adds more samples there as a result of the high variance density. Of course,
these figures show only one realization of the algorithm and the decomposition
is dependent on the initial ED. Therefore, it is necessary to repeat the
algorithm many times with random initial ED to assess convergence.

Figure 4: (a), (b) The adapted domain and ED before (iteration 4) and after
(iteration 11) exploration and discovery of the exponential part of the
mathematical model. (c) Final surrogate models from global PCE and DAL-PCE.
(d) Convergence plot comparing the mean square error for global PCE SSE, and
DAL-PCE. The convergence plots for Global PCE and DAL-PCE show continuous mean
value $\pm\sigma$ intervals from 100 repeated trials, while those for SSE are
plotted for several discrete ED sizes.

Fig. 4d shows convergence of the error $\epsilon$ from 100 repeated trials.
The single global PCE is unable to accurately approximate the original
function even when using high $p$ and thus the $\epsilon$ does not converge,
as expected. Both methods based on domain decomposition (DAL-PCE and SSE)
achieve great accuracy already for $200$ samples. However, the DAL-PCE
consistently has $1$–$2$ orders of magnitude higher accuracy than SSE for the
given number of samples. Moreover, increase in variance of $\epsilon$ is, in
general, slower in DAL-PCE than in SSE. Fast increment in variance of SSE can
be seen also in the original paper [2]. Finally, we again observe that
convergence is continuous with DAL-PCE, where convergence can only be assessed
at discrete sample sizes with SSE through a new analysis. All of these
advantages of the DAL-PCE can be attributed to the active learning, which both
explores the space and exploits the behavior of the function to decompose the
domain and add samples. Although active learning might lead to lower accuracy
(higher $\epsilon$) initially (for small ${n_{\mathrm{sim}}}=10$–$20$) as it
is dominated by exploration, it rapidly improves once it identifies important
features and begins to favor exploitation.

###  4.2 Two-dimensional Singularity

Figure 5: Results for the $2$-dimensional Singularity function: a) original
mathematical model, b) approximation via DAL-PCE (background color), current
domain division and the corresponding ED, c) local LOO-CV
$Q^{2}_{\mathcal{D}_{i}}$ and $\Theta_{i}$ value for each sub-domain, d)
convergence plots for DAL-PCE, Global PCE, and SSE showing the mean value and
$\pm\sigma$ interval. Convergence plots for SSE show the mean $\pm\sigma$ at
discrete sample sizes.

The second example involves a $2$D function with mirrored quarter-circle arc
line singularities [1]. The form of the function is give by:

$\centering f(\boldsymbol{X})=\frac{1}{\lvert
0.3-X_{1}^{2}-X_{2}^{2}\rvert+\delta}-\frac{1}{\lvert
0.3-(1-X_{1})^{2}-(1-X_{2})^{2}\rvert+\delta},\quad\boldsymbol{X}\sim\pazocal{U}[0,1]^{2},\@add@centering$
(26)

where the strength of the singularities is controlled by the parameter
$\delta$, which we set as $\delta=0.1$. The singularities in this example
represent a challenging task for a global PCE even with high order, due to the
well-known Gibbs phenomenon [49]. It is thus beneficial to identify the
location of the singularity, locally decompose the domain, and construct low-
order local PCEs.

Fig. 5 illustrates the decomposition and DAL-PCE approximation at a given
stage of the computation. Panel a) visualizes the true values of the function
via a background color. The same coloring scheme is used in panel b) for the
pointwise information available in the current ED (small circles) and for the
function approximation via DAL-PCE by the background color. Panels b) and c)
show also the final domain decomposition. The symmetry in the decomposition
documents the great convergence of the DAL-PCE thanks to an adaptive
decomposition described in the previous section. Plot c) shows the local
$Q^{2}_{\mathcal{D}_{i}}$ error in each individual sub-domain (darker color
corresponds to higher local error). These local errors clearly show
localization of the prediction error to very small areas near singularities,
which are continually being refined. The color of the small solid squares in
the center of each sub-domains shows the $\Theta_{i}$ value for that sub-
domain.

Finally, the convergence plot in Fig. 5d) shows that both DAL-PCE and SSE
outperform the global PCE, as expected. The SSE performs comparable to or
slightly better than DAL-PCE for small $N_{\mathrm{sim}}$, but the DAL-PCE
begins to outperform SSE as $N_{\mathrm{sim}}$ grows thanks to the active
learning approach that targets samples in the vicinity of the singularities.
Note that the error converges for both SSE and DAL-PCE as we approach 1000
samples and does not seem to substantially reduce after this. This is due to
the fundamental limitation of trying to approximate this singularity, even
locally, with low-order polynomials.

###  4.3 $M$-dimensional Discontinuity

Figure 6: Results for the $2$-dimensional discontinuiy function: a) original
mathematical model, b) approximation via DAL-PCE and ED, c) local LOO-CV
$Q^{2}_{\mathcal{D}_{i}}$ and $\Theta_{i}$ value for each sub-domain, d)
convergence plots for DAL-PCE, Global PCE, and SEE showing the mean value and
$\pm\sigma$ interval. Convergence plots for SSE show the mean $\pm\sigma$ at
discrete sample sizes.

The third example investigates the role of dimensionality on the performance
of the proposed DAL-PCE. The following discontinuous function is defined for
an arbitrary number of input random variables $M$ [26]:

$f(\boldsymbol{X})=\begin{cases}\sin{(X_{1}\pi)}\sin{(X_{2}\pi)}&\text{if
$x_{1}\leq 0.5$ and $x_{2}\leq 0.5$}\\\
\sum_{i=3}^{M}X_{i}&\text{otherwise}\end{cases},\quad\boldsymbol{X}\sim\pazocal{U}[0,1]^{M}.$
(27)

This function has a discontinuity in the first two input random variables,
which can be seen in Fig. 6a. A single global PCE cannot accurately
approximate the function because of the discontinuity, although the function
$f(\boldsymbol{X})$ can be easily approximated by two separate PCEs in the two
regions for which the definitions differ. But, this requires a priori
knowledge of the discontinuity location. Since the location of the
discontinuity is assumed to be unknown, this function is a good example for
domain adaptation using DAL-PCE.

The detailed results for a 2D version of this problem are depicted in Fig. 6
in identical form as in the previous example. Note that the local $Q^{2}_{i}$
errors Fig. 6c show perfect accuracy in the part of the input random space
where $f(\boldsymbol{X})=0$ and thus the associated sub-domains are not
preferred for further decomposition. The convergence plot in Fig. 6d confirms
that a single global PCE is not able to create an accurate approximation and
adding more points to ED does not lead to significant improvements in the
approximation. The mean values of errors $\epsilon$ associated to the proposed
DAL-PCE approach are significantly lower in comparison to SSE ($1$–$2$ orders
of magnitude) similarly as in the first example, though the convergence trend
is similar for both methods. SSE, however, uses a random splitting routine.
This can lead to very high variance of results, since the accuracy is highly
dependent on the pattern of the decomposed input random space. This clearly
shows the advantage of an active learning approach.

Figure 7: Convergence plots for the $M$-dimensional function: a)
$3$-dimensional version, b) $5$-dimensional version, c) $6$-dimensional
version, and d) $8$-dimensional version. Convergence plots for the DAL-PCE and
global PCE show the mean value $\pm\sigma$ interval. Convergence plots for SSE
also show the mean $\pm\sigma$, but at discrete sample sizes.

The influence of dimensionality $M$ on convergence of the DAL-PCE, SSE, and
global PCE is studied in Fig. 7 for a) 3, b) 5, c) 6, and d) 8 input random
variables. As the domain dimension increases, the linear part of the function
$f(\boldsymbol{X})$ occupies an increasing proportion of the domain while the
discontinuity remain low-dimensional. The proposed DAL-PCE greatly improves
the convergence because it is able to identify an ideal decomposition and
local samples to resolve the discontinuity. For low-dimensions ($M=2,3$), SSE
error $\epsilon$ shows a decreasing trend that is better than global PCE but
has an extremely high variance. This is caused by a lack of control in sample
placement. The domain decomposition in SSE is a product of sample location and
without active learning to guide sample placement, SSE will sometimes produce
a very good decomposition and sometimes a very poor decomposition. Meanwhile,
the proposed DAL-PCE errors have comparably low variance for low-dimensions
and consistently have accuracy comparable to, or better than, the best SSE
realizations.

As the dimension, $M$, increases the DAL-PCE is able to maintain a very high
level of accuracy, while the accuracy degrades completely for the SSE such
that it is comparable to the global PCE. The DAL-PCE is able to maintain its
low error because the discontinuity remains low-dimensional and the active
learning process is able to target this region for domain refinement and
sampling. This means that the DAL-PCE remains largely independent of the
problem dimension, and instead depends predominantly on the intrinsic
dimension of the discontinuous/nonlinear features of the model. The
performance of SSE, on the other hand, degrades with dimension because its
domain decomposition depends only on a set of a priori specified points that
are not selected in a way that is aware of the important features of the
model. Consequently, as the dimension increases the algorithm becomes less
likely to refine the domain appropriately around an embedded low-dimensional
feature. We remark that this desirable scalable convergence trend of the DAL-
PCE is not likely a universal property, as the trend may break down in
problems where the intrinsic dimension of the discontinuity/nonlinearity is
high or where the discontinuity occupies a very small proportion of the domain
– in which case exploration of the space to find the important feature may
take a very large number of samples.

In the present example, the discontinuity in the function given in Eq. (27)
lies at $x_{1}=0.5$ and $x_{2}=0.5$, which corresponds to the exact location
where the domain will be split for both SSE and during the early iterations of
the DAL-PCE. One might argue that this presents an unreasonable advantage for
the proposed algorithm. We therefore modified the function such that the
discontinuity lies at $x_{1}=0.61$ and $x_{2}=0.61$. Fig. 8 shows the
convergence for the DAL-PCE and SSE for this modified function with varying
dimension, $M$. The absolute errors $\epsilon$ exhibit slower decrease,
especially for dimensions $M=3$ and $M=5$. However, the proposed active
learning still leads to superior results (especially for higher dimensions as
in the previous case). Note that there are visible spikes in the DAL-PCE
convergence graph for the 3-dimensional example. Although the results were
statistically processed, these spikes are caused by the restart adaptivity
occurring at the same ${N_{\mathrm{sim}}}$ in each replication. In this case,
the optimal decomposition pattern is very complicated and therefore the
algorithm activates the restart adaptivity frequently (after multiples of
$n_{r}$ steps), until it finds a suitable pattern to continue convergence. SSE
in the 3- and 5-dimensional cases has higher mean error and significantly
lower variance in comparison to the previous example. This is caused by the
fact that the modified discontinuity location no longer lies along the
boundary of the domain decomposition. In the previous example, some SSE
realizations achieved near-perfect accuracy because the domain was
coincidentally divided along the discontinuity.

Figure 8: Convergence plots for the modified $M$-dimensional function: a)
$3$-dimensional version, b) $5$-dimensional version, c) $6$-dimensional
version, and d) $8$-dimensional version. Convergence plots for the DAL-PCE and
global PCE show the mean value $\pm\sigma$ interval. Convergence plots for SSE
also show the mean $\pm\sigma$, but at discrete sample sizes. Figure 9:
Convergence plots for DAL-PCE and SSE with additional boxplots for SSE showing
the median, lower and upper quartiles and outliers for: a) the 3D example with
discontinuity at $x_{1}=0.5$ and $x_{2}=0.5$, b) the 3D example with
discontinuity at $x_{1}=0.61$ and $x_{2}=0.61$.

This phenomenon is investigated more closely in Fig. 9, which compares number
of outliers in both versions of 3D examples. In addition to the mean
$\pm\sigma$ seen previously, the figure also shows standard boxplots for SSE
(median along with lower and upper quartiles) and the corresponding number of
“extreme” realizations producing very high accuracy (top axis) for a) the
original position of discontinuity; and b) discontinuity at $x_{1}=0.61$ and
$x_{2}=0.61$. As can be seen, in panel a) there are many outliers producing
$\epsilon<-7$, which effectively decreases $\mu$ relative to the median while
also significantly increasing the variance. In contrast DAL-PCE has no
outliers and it leads to very consistent results. In panel b), there are no
outliers for either SSE or DAL-PCE and the results are thus consistent with
low variance for both methods.

###  4.4 Asymmetric shallow von Mises truss

In this section, we demonstrate the relevance of the proposed method for a
representative engineering example exhibiting discontinuous response. Consider
the shallow two-bar planar truss subjected to a vertical load at its top
joint, as presented in [50] and illustrated in Fig. 10a.

The truss is formed by two prismatic bars made of a hard wood (density $800$
kg/m3, modulus of elasticity $E=12$ GPa). There are two variables in the
studied von Mises truss: (i) the loading vertical force $F$, and (ii) a half
sine-wave imperfection of the left bar having magnitude $\delta$, see the
sketch in Fig. 10a. The load is applied dynamically as a step function at time
zero for an unlimited duration. The structure is modeled, as illustrated in
Fig. 10b. In particular, the mass of the bar is concentrated in 21 mass
points, including the supports and the loading point. These mass points are
connected via $10+10$ translational springs representing the normal stiffness
of the true bars. The pairs of the axial members are connected via rotational
spring having zero moment for a zero angle between adjacent bars. The only
exceptions are the loading ans support points where there are no rotational
springs attached (hinges). The damping is associated with the mass points via
linear viscous damping coefficient set to $11\leavevmode\nobreak\
\text{N}\cdot\text{s}/(\text{kg}\cdot\text{m})$ approximating the relative
damping of about 3%. Explicit dynamics solver FyDiK [51, 52] was used to solve
the equations of equilibrium at the mass points. The numerical solution lasts
to up to two seconds, which is the time needed for almost complete
stabilization of the solution (kinetic energy drops below a negligible
threshold).

Since the structure is very shallow, sudden application of the vertical force
can cause snap-through buckling, wherein the loading point drops down between
the supports and the members switch from a state of compression to tensile
stresses in the final stable state. We specifically study the horizontal
coordinate $y_{F}$ of the loading point after the dynamic response stabilizes
to the final deformed shape. The force $F\in(31.6,772.6)$ kN and initial
imperfection $\delta\in(-0.4,0.4)$ m are treated as uniform random variables
mapped to the unit square such that the model input
$\mathbf{X}\sim\pazocal{U}[0,1]^{2}$. Because of the potential snap-through
buckling, the solution is discontinuous as illustrated in Fig. 10c. On each
side of the discontinuity, the solution $y_{F}$ is smooth and slowly-varying
having values near +1 m and -1 m, respectively. Note that the output is _not
symmetric_ with respect to $\delta=0$ because the dynamical response evolves
differently for concave and convex initial displacements.

The sharp boundary between the buckled and unbuckled regions, shown in Fig.
11a cause global PCE to produce poor approximations that are vulnerable to the
Gibbs phenomenon, similar to the example in subsection 5. This is shown by the
convergence plots in Fig. 11d comparing global PCE, DAL-PCE, and SSE. Clearly,
the complexity of this example and the complicated shape of the discontinuity
limits the accuracy of all the surrogate models. The proposed DAL-PCE achieves
low accuracy for small sample sizes because the corresponding small number of
sub-domains and low-order PCEs are unable to sufficiently approximate the
boundary. Therefore, the global PCE and SSE (with a low number of embedding
levels) are initially better. With increasing number of samples, the proposed
DAL-PCE approach leads to superior results because the active learning is able
to resolve the discontinuity as illustrated in Fig. 11b, which shows the
domain decomposition and approximation after 2000 samples. Fig. 11c shows the
corresponding LOO-CV errors for each subdomain, demonstrating the errors are
confined to small, localized regions near the boundary.

Figure 10:  Asymmetric shallow von Mises truss. a) Initial geometry with two
random variables $F$ and $\delta$; b) illustrative sketch of the discrete
dynamical model and the meaning of output variable $y_{F}$, c) illustration of
the discontinuous response function of the two input variables.  Figure 11:
Results for the von Misses truss example: a) original mathematical model
(numerical solution), b) approximation via DAL-PCE and ED, c) local LOO-CV
$Q^{2}_{\mathcal{D}_{i}}$ and $\Theta_{i}$ value for each sub-domain, d)
convergence plots for DAL-PCE, Global PCE, and SSE showing the mean value and
$\pm\sigma$ interval; convergence plots for SSE show the mean $\pm\sigma$ at
discrete sample sizes.

##  5 Discussion & Future Work

The proposed DAL-PCE approach is a general methodology for the decomposition
of the input random space and construction of localized PCEs using active
learning. The proposed active learning is based on a novel $\Theta$ criterion
that optimally balances global _exploration_ with local _exploitation_ of the
model. Although this paper presents one specific learning algorithm, the
methodology is general and amenable to modifications to reflect the specific
user’s needs. The whole process can be divided into two tasks: A)
decomposition of the input random space and B) construction of localized PCEs;
and both can be easily modified as discussed further:

  1. A)

The most important sub-domain $\mathcal{D}_{i}$ is identified by extended
$\Theta$ according to Eq. (17) evaluated for a large number of global
candidates. In this paper, we use standard LHS for candidate generation, but
it may be beneficial to use different sampling methods that produce more
uniform coverage of the whole input random space (see e.g. [53, 54, 45]).
Although it is generally possible to generate a large number of candidates, it
might be challenging to uniformly cover the entire input random space,
especially in high dimensions. Thus, one can use any sampling technique
suitable for a specific example, e.g. [55].

Once the $\mathcal{D}_{i}$ is identified via Eq. (17), it is either divided
(providing it contains enough ED points) or the sample is extended inside it,
to achieve a better PCE approximation. The simplest division occurs by
splitting the volume into two parts of identical hypervolume in the direction
of the highest first-order Sobol’ index. However, the algorithm can
accommodate various different approaches. For example, it is possible to
divide the $\mathcal{D}_{i}$ into a higher number of sub-domains, not just
two. Moreover, instead of splitting the domain into parts of equal
hypervolume, other criteria can be used. For example, the cutting plane can be
positioned so to split the domain variance into equal parts.

  2. B)

The user can choose to employ any existing method to construct the non-
intrusive PCEs, including various sparse solvers or adaptive algorithms, which
may be preferable for certain applications [12]. For example, we use LARS with
OLS. However, it is generally more efficient to use active learning based on
the $\Theta$ criterion for PCE as shown in [1], which employs variance-based
sequential sampling. This improvement can be integrated within the DAL-PCE to
make local PCE more efficient in each subdomain, and thereby improving the
overall convergence. The can be compounded by the use of advanced sampling
techniques within the subdomains such as Coherence D-optimal sampling [40,
41].

As seen from the previous paragraphs, the whole algorithm can be adapted for
specific needs reflecting the characteristics of a given mathematical model,
such as dimensionality, sparsity, non-linearity etc., by simply exchanging
components of the proposed algorithm for suitable existing (or new)
techniques. Note that even after the modification, the whole methodology based
on $\Theta$ criterion is still valid and can be used for uncertainty
quantification and surrogate modelling as described in this paper. Moreover,
in comparison to SSE, the DAL-PCE sequentially adds points and divides the
sub-domains one-by-one based on information obtained from the previous
iteration.

Another significant advantage of the DAL-PCE is that it provides estimates of
the local errors, $Q_{\mathcal{D}_{i}}$, associated with each sub-domain.
Since localized PCEs are constructed independently, local errors estimate the
local accuracy of the surrogate model directly, and can be assembled to
provide global error measures. Naturally, local accuracy is very important
information that can be used for further probabilistic analysis and active
learning. Although this paper does not propose any specific approach for
further processing of this information, it could serve as a main ingredient
for various active learning algorithms. For example, it could be directly used
to predict uncertainty in industrial applications and possibly extend the ED
in a sub-domain of interest.

Finally, an important topic of further research is to study the behavior of
the proposed criterion in higher dimensions. In particular, the geometrical
terms $l_{c,s}^{M}$ and $\mathcal{W}_{i}$ likely cause poor convergence in
high dimensions. Although some preliminary results focused on investigating of
$l_{c,s}^{M}$ in high dimensions was previously performed in the paper [1]
proposing the original $\Theta$ criterion, it is still necessary to perform an
extensive study of its behavior as well as investigating the influence of
$\mathcal{W}_{i}$, which may need to be reformulated for high dimensions.

##  6 Conclusion

The paper presented a novel approach, domain adaptively localzed PCE, for the
adaptive sequential construction of localized PCEs based on active learning
and decomposition of the input random space. It combines adaptive sequential
sampling based on the recently proposed $\Theta$ criterion to maintain the
balance between exploration of the input random space and exploitation of the
current characteristics of the PCE together with the adaptive sequential
decomposition of the input random space creating sub-domains approximated by
local surrogate models. The methodology offers a general technique that can be
easily adapted or modified for specific functions extending its applicability.
The performance of the proposed methodology was validated on several numerical
examples of increasing complexity investigating different aspects of the
algorithm and leading to superior results in comparison to a single global PCE
and the recently proposed SSE.

## Acknowledgments

The first author acknowledge financial support provided by the Czech Science
Foundation under project number 22-00774S. Additionally, the major part of
this research was conducted during the research stay of the first author at
Johns Hopkins University supported by the project International Mobility of
Researchers of Brno University of Technology, Czechia under project No.
EF18_053/0016962.

## References

  * [1] L. Novák, M. Vořechovský, V. Sadílek, M. D. Shields, Variance-based adaptive sequential sampling for polynomial chaos expansion, Computer Methods in Applied Mechanics and Engineering 386 (2021) 114105.  doi:10.1016/j.cma.2021.114105.
  * [2] S. Marelli, P.-R. Wagner, C. Lataniotis, B. Sudret, STOCHASTIC SPECTRAL EMBEDDING, International Journal for Uncertainty Quantification 11 (2) (2021) 25–47.  doi:10.1615/int.j.uncertaintyquantification.2020034395.
  * [3] N. Wiener, The homogeneous chaos, American Journal of Mathematics 60 (4) (1938) 897–936.  doi:10.2307/2371268.
  * [4] G. Blatman, B. Sudret, Adaptive sparse polynomial chaos expansion based on least angle regression, Journal of Computational Physics 230 (6) (2011) 2345–2367.  doi:10.1016/j.jcp.2010.12.021.
  * [5] R. G. Ghanem, P. D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer New York, 1991.  doi:10.1007/978-1-4612-3094-6.
  * [6] N.-Z. Chen, C. Guedes Soares, Spectral stochastic finite element analysis for laminated composite plates, Computer methods in Applied Mechanics and Engineering 197 (51) (2008) 4830–4839.  doi:10.1016/j.cma.2008.07.003.
  * [7] L. Novak, D. Novak, Surrogate modelling in the stochastic analysis of concrete girders failing in shear, in: Proc. of the Fib Symposium 2019: Concrete - Innovations in Materials, Design and Structures, 2019, pp. 1741–1747.
  * [8] B. Sudret, Global sensitivity analysis using polynomial chaos expansions, Reliability Engineering & System Safety 93 (7) (2008) 964–979.  doi:10.1016/j.ress.2007.04.002.
  * [9] T. Crestaux, O. L. Maître, J.-M. Martinez, Polynomial chaos expansion for sensitivity analysis, Reliability Engineering & System Safety 94 (7) (2009) 1161–1172.  doi:10.1016/j.ress.2008.10.008.
  * [10] A. Cohen, G. Migliorati, Optimal weighted least-squares methods, The SMAI journal of computational mathematics 3 (2017) 181–203.  doi:10.5802/smai-jcm.24.
  * [11] A. C. Narayan, J. Jakeman, T. Zhou, A Christoffel function weighted least squares algorithm for collocation approximations, Math. Comput. 86 (2017) 1913–1947.  doi:10.1090/mcom/3192.
  * [12] N. Lüthen, S. Marelli, B. Sudret, Sparse polynomial chaos expansions: Literature survey and benchmark, SIAM/ASA Journal on Uncertainty Quantification 9 (2) (2021) 593–649.  doi:10.1137/20M1315774.
  * [13] B. Echard, N. Gayton, M. Lemaire, AK-MCS: An active learning reliability method combining kriging and monte carlo simulation, Structural Safety 33 (2) (2011) 145–154.  doi:10.1016/j.strusafe.2011.01.002.
  * [14] L. Shi, B. Sun, D. S. Ibrahim, An active learning reliability method with multiple kernel functions based on radial basis function, Structural and Multidisciplinary Optimization 60 (1) (2019) 211–229.  doi:10.1007/s00158-019-02210-0.
  * [15] X. Yang, X. Cheng, Active learning method combining kriging model and multimodal-optimization-based importance sampling for the estimation of small failure probability, International Journal for Numerical Methods in Engineering 121 (21) (2020) 4843–4864.  doi:10.1002/nme.6495.
  * [16] S. Marelli, B. Sudret, An active-learning algorithm that combines sparse polynomial chaos expansions and bootstrap for structural reliability analysis, Structural Safety 75 (2018) 67–74.  doi:10.1016/j.strusafe.2018.06.003.
  * [17] Y. Zhou, Z. Lu, W. Yun, Active sparse polynomial chaos expansion for system reliability analysis, Reliability Engineering & System Safety 202 (2020) 107025\.  doi:10.1016/j.ress.2020.107025.
  * [18] K. Cheng, Z. Lu, Active learning polynomial chaos expansion for reliability analysis by maximizing expected indicator function prediction error, International Journal for Numerical Methods in Engineering 121 (14) (2020) 3159–3177.  doi:10.1002/nme.6351.
  * [19] N. Fajraoui, S. Marelli, B. Sudret, Sequential design of experiment for sparse polynomial chaos expansions, SIAM/ASA Journal on Uncertainty Quantification 5 (1) (2017) 1061–1085.  doi:10.1137/16m1103488.
  * [20] M. Thapa, S. B. Mulani, R. W. Walters, Adaptive weighted least-squares polynomial chaos expansion with basis adaptivity and sequential adaptive sampling, Computer methods in Applied Mechanics and Engineering 360 (2020) 112759.  doi:10.1016/j.cma.2019.112759.
  * [21] M. D. Shields, Adaptive Monte Carlo analysis for strongly nonlinear stochastic systems, Reliability Engineering & System Safety 175 (2018) 207–224.  doi:10.1016/j.ress.2018.03.018.
  * [22] Y. Zhou, Z. Lu, K. Cheng, C. Ling, An efficient and robust adaptive sampling method for polynomial chaos expansion in sparse bayesian learning framework, Computer Methods in Applied Mechanics and Engineering 352 (2019) 654–674.  doi:10.1016/j.cma.2019.04.046.
  * [23] J. Zhang, W. Gong, X. Yue, M. Shi, L. Chen, Efficient reliability analysis using prediction-oriented active sparse polynomial chaos expansion, Reliability Engineering & System Safety 228 (2022) 108749.  doi:10.1016/j.ress.2022.108749.
  * [24] M. K. Deb, I. M. Babuška, J. Oden, Solution of stochastic partial differential equations using galerkin finite element techniques, Computer Methods in Applied Mechanics and Engineering 190 (48) (2001) 6359–6372.  doi:10.1016/S0045-7825(01)00237-7.
  * [25] J. A. Witteveen, G. Iaccarino, Simplex stochastic collocation with random sampling and extrapolation for nonhypercube probability spaces, SIAM Journal on Scientific Computing 34 (2) (2012) A814–A838.
  * [26] A. Bhaduri, Y. He, M. D. Shields, L. Graham-Brady, R. M. Kirby, Stochastic collocation approach with adaptive mesh refinement for parametric uncertainty analysis, Journal of Computational Physics 371 (2018) 732–750.  doi:10.1016/j.jcp.2018.06.003.
  * [27] P.-R. Wagner, S. Marelli, I. Papaioannou, D. Straub, B. Sudret, Rare event estimation using stochastic spectral embedding, Structural Safety 96 (2022) 102179\.  doi:10.1016/j.strusafe.2021.102179.
  * [28] X. Wan, G. E. Karniadakis, An adaptive multi-element generalized polynomial chaos method for stochastic differential equations, Journal of Computational Physics 209 (2) (2005) 617–642.  doi:10.1016/j.jcp.2005.03.023.
  * [29] D. Xiu, G. E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM Journal on Scientific Computing 24 (2) (2002) 619–644.  doi:10.1137/s1064827501387826.
  * [30] L. Novák, On distribution-based global sensitivity analysis by polynomial chaos expansion, Computers & Structures 267 (2022) 106808.  doi:10.1016/j.compstruc.2022.106808.
  * [31] W. Gautschi, On generating orthogonal polynomials, SIAM Journal on Scientific and Statistical Computing 3 (3) (1982) 289–317.  doi:10.1137/0903018.
  * [32] B. Efron, T. Hastie, I. Johnstone, R. Tibshirani, Least angle regression, The Annals of Statistics 32 (2) (2004) 407–451.  doi:10.2307/3448465.
  * [33] J. A. Tropp, A. C. Gilbert, Signal recovery from random measurements via orthogonal matching pursuit, IEEE Transactions on Information Theory 53 (12) (2007) 4655–4666.  doi:10.1109/tit.2007.909108.
  * [34] S. Ji, Y. Xue, L. Carin, Bayesian compressive sensing, IEEE Transactions on Signal Processing 56 (6) (2008) 2346–2356.  doi:10.1109/TSP.2007.914345.
  * [35] G. Blatman, B. Sudret, An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis, Probabilistic Engineering Mechanics 25 (2) (2010) 183–197.  doi:10.1016/j.probengmech.2009.10.003.
  * [36] A. Olivier, D. Giovanis, B. Aakash, M. Chauhan, L. Vandanapu, M. D. Shields, UQpy: A general purpose python package and development environment for uncertainty quantification, Journal of Computational Science 47 (2020) 101204\.
  * [37] S. Marelli, B. Sudret, UQLab: A framework for uncertainty quantification in Matlab, in: Vulnerability, Uncertainty, and Risk, 2014, pp. 2554–2563.  doi:10.1061/9780784413609.257.
  * [38] M. D. McKay, W. J. Conover, R. J. Beckman, A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics 21 (1979) 239–245.  doi:10.1080/00401706.1979.10489755.
  * [39] W. Conover, On a better method for selecting input variables, unpublished Los Alamos National Laboratories manuscript, reproduced as Appendix A of “Latin Hypercube Sampling and the Propagation of Uncertainty in Analyses of Complex Systems” by J.C. Helton and F.J. Davis, Sandia National Laboratories report SAND2001-0417, printed November 2002. (1975).
URL https://prod-ng.sandia.gov/techlib-noauth/access-
control.cgi/2001/010417.pdf

  * [40] J. Hampton, A. Doostan, Compressive sampling of polynomial chaos expansions: Convergence analysis and sampling strategies, Journal of Computational Physics 280 (2015) 363–386.  doi:10.1016/j.jcp.2014.09.019.
  * [41] P. Diaz, A. Doostan, J. Hampton, Sparse polynomial chaos expansions via compressed sensing and D-optimal design, Computer methods in Applied Mechanics and Engineering 336 (2018) 640–666.  doi:10.1016/j.cma.2018.03.020.
  * [42] J. F. Koksma, Een algemeene stelling uit de theorie der gelijkmatige verdeeling modulo 1, Mathematica B 11 (1942/1943) 7–11.
  * [43] M. Johnson, L. Moore, D. Ylvisaker, Minimax and maximin distance designs, Journal of Statistical Planning and Inference 2 (26) (1990) 131–148.  doi:10.1016/0378-3758(90)90122-B.
  * [44] L. Pronzato, Minimax and maximin space-filling designs: some properties and methods for construction, Journal de la Société Française de Statistique 158 (1) (2017) 7–36.
  * [45] J. Eliáš, M. Vořechovský, V. Sadílek, Periodic version of the minimax distance criterion for Monte Carlo integration, Advances in Engineering Software 149 (2020) 102900.  doi:10.1016/j.advengsoft.2020.102900.
  * [46] M. Rosenblatt, Remarks on a multivariate transformation, The Annals of Mathematical Statistics 23 (3) (1952) 470–472.  doi:10.1214/aoms/1177729394.
  * [47] A. Nataf, Détermination des distributions de probabilité dont les marges sont données, Comptes Rendus de l’Académie des Sciences 225 (1962) 42–43.
  * [48] F. Wang, H. Li, System reliability under prescribed marginals and correlations: Are we correct about the effect of correlations?, Reliability Engineering & System Safety 173 (2018) 94–104.  doi:10.1016/j.ress.2017.12.018.
  * [49] J. M. Davis, P. Hagelstein, Gibbs phenomena for some classical orthogonal polynomials, Journal of Mathematical Analysis and Applications 505 (1) (2022) 125574\.
  * [50] M. Vořechovský, Reliability analysis of discrete-state performance functions via adaptive sequential sampling with detection of failure surfaces, Computer Methods in Applied Mechanics and Engineering 401 (2022) 115606\.  doi:10.1016/j.cma.2022.115606.
  * [51] P. Frantík, FyDik - a software for interactive simulations of dissipative nonlinear dynamical systems based on physical discretization, http://fydik.kitnarf.cz/ (2000–2022).
  * [52] P. Frantík, Simulation of the stability loss of the von Mises truss in an unsymmetrical stress state, Engineering Mechanics 14 (3) (2007) 155–161.
  * [53] M. Vořechovský, J. Eliáš, Modification of the maximin and $\phi_{p}$ (phi) criteria to achieve statistically uniform distribution of sampling points, Technometrics 62 (3) (2020) 371–386.  doi:10.1080/00401706.2019.1639550.
  * [54] M. Vořechovský, J. Mašek, J. Eliáš, Distance-based optimal sampling in a hypercube: Analogies to N-body systems, Advances in Engineering Software 137 (2019) 102709.  doi:10.1016/j.advengsoft.2019.102709.
  * [55] M. Vořechovský, J. Mašek, Distance-based optimal sampling in a hypercube: Energy potentials for high-dimensional and low-saturation designs, Advances in Engineering Software 149 (2020) 102880.  doi:10.1016/j.advengsoft.2020.102880.