# Sparse Polynomial Chaos Expansion for Universal Stochastic Kriging

J. C. García-Merino  C. Calvo-Jurado  E. García-Macías  Department of
Mathematics, School of Technology, 10003, Cáceres, Spain  Department of
Structural Mechanics and Hydraulic Engineering, University of Granada, Campus
de Fuentenueva s/n, 18071 Granada, Spain

###### Abstract

Surrogate modelling techniques have opened up new possibilities to overcome
the limitations of computationally intensive numerical models in various areas
of engineering and science. However, while fundamental in many engineering
applications and decision-making, the incorporation of uncertainty
quantification into meta-models remains a challenging open area of research.
To address this issue, this paper presents a novel stochastic simulation
approach combining sparse polynomial chaos expansion (PCE) and Stochastic
Kriging (SK). Specifically, the proposed approach adopts adaptive sparse PCE
as the trend model in SK, achieving both global and local prediction
capabilities and maximizing the role of the stochastic term to conduct
uncertainty quantification. To maximize the generalization and computational
efficiency of the meta-model, the Least Angle Regression (LAR) algorithm is
adopted to automatically select the optimal polynomial basis in the PCE. The
computational effectiveness and accuracy of the proposed approach are
appraised through a comprehensive set of case studies and different quality
metrics. The presented numerical results and discussion demonstrate the
superior performance of the proposed approach compared to the classical
ordinary SK model, offering high flexibility for the characterization of both
extrinsic and intrinsic uncertainty for a wide variety of problems.

###### keywords:

Kriging , Least Angle Regression , Polynomial chaos expansion , Stochastic
simulation , Surrogate modelling , Uncertainty quantification

††journal: Journal of Computational and Applied Mathematics

##  1 Introduction

Computer models, commonly referred to as simulators, have become ubiquitous
assets in most areas of engineering and science. In cases where systems or
processes are assumed to be predictable and devoid of inherent randomness, a
deterministic simulator can be adopted. Nonetheless, deterministic approaches
appear inadequate for facing real-world problems characterized by intrinsic or
non-parameterisable randomness, for which the adoption of stochastic
simulators becomes necessary for proper modelling. A stochastic simulator
incorporates randomness or uncertainty in its predictions. Thus, while a
deterministic simulator always produces the same prediction for a fixed input,
in stochastic simulators, each simulation run may produce different results as
a result of the inherent randomness in the model. Such models have gained
significant popularity with manifold applications in diverse fields such as
epidemiology [1, 2], materials science [3, 4], or reliability analysis [5],
enabling researchers to make robust predictions to assist decision-making for
complex systems. However, for both deterministic and stochastic
computationally demanding simulators, formidable challenges arise when
implemented in iterative procedures such as optimization [6], sensitivity
analysis [7], parameter inference [8], or uncertainty propagation [9]. Despite
the introduction of numerous surrogate models or meta-models in the literature
over the past few years to replace intensive numerical models and alleviate
the computational burden, their success mostly limits to deterministic
simulators, while the consideration of stochasticity remains an open research
topic.

There is a vast literature on the development of surrogate models for
deterministic simulators [10, 11]. In particular, non-intrusive models such as
polynomial chaos expansion (PCE) [12], radial basis functions [13] or Kriging
[14] have received significant attention from the scientific community due to
their effectiveness in addressing resource-intensive simulations. The
construction of these meta-models is usually conducted on the basis of a
training dataset, also referred to as the experimental design (ED) [14]. Such
a training dataset contains the evaluations of the forward model, representing
the quantity of interest (QoI), for a set of combinations of the input
parameters covering the parameter space. In the context of stochastic
simulators, the QoI usually presents a certain level of uncertainty, stemming
from the stochastic nature of the model parameters and/or the presence of
epistemic uncertainties in the model predictions [15]. A first approach to
meta-model such processes involves estimating the probability density function
(PDF) of the simulator’s predictions through functional decomposition [16].
Alternatively, Stochastic Kriging (SK) [17] has gained significant recognition
as a general meta-modelling tool for representing stochastic simulation
response surfaces. Defined as an extension of Kriging, SK accounts for
heteroscedastic noise in the response by estimating the sample variance at
each point of the ED using multiple runs of the simulator for each specific
input. The effectiveness of SK has been demonstrated on a broad range of
applications, including risk and reliability analysis [18, 19], industrial
welding processes [20], and game theoretic simulations [21], to mention a few.

One of the primary challenges encountered in SK is the requirement of a
significant number of replications at each ED point to accurately estimate the
sample variance [15], which often leads to serious computational demands. One
approach to alleviate such a limitation regards the optimal definition of the
ED. In the case of stochastic simulators, sampling efforts are divided into
two parts: (i) sampling new ED locations, and (ii) acquiring additional
replications at existing points. In this light, numerous sequential sampling
techniques have been presented (refer to [22, 23] for an extensive state-of-
the-art review). On the other hand, various extensions of the classical SK
have been proposed in the literature. It is worth noting the work by Binois et
al. [24], who proposed modelling the input-dependent noise as a separate
Gaussian process. Additionally, they employed computationally efficient
techniques such as the Woodbury identity to address the challenges encountered
in SK for large EDs. In a similar vein, Hao et al. [19] introduced a Nested SK
model, which decouples in an iterative way the response noise parameters from
the deterministic Kriging parameters used for calibrating the meta-model. In a
similar context, Cheng et al. [25] proposed a generalized polynomial chaos-
stochastic Kriging (gPCE SK) approach, which combines the benefits of PCE and
SK. The gPCE SK method proposed by those authors utilizes a local PCE
surrogate model at each design point, which is calibrated using a small number
of simulation runs. By relying on these local surrogate models, numerous
replications can be performed at each training point, significantly reducing
the computational burden involved in the calibration of the SK model. However,
it is important to note that this model assumes a distinction between
controllable design variables and non-controllable noise variables, the latter
causing the variations between replicates at each fixed design point.
Furthermore, it relies on knowledge of the PDFs of the uncontrollable
variables, which may be a strong assumption in many real-life applications.

In light of the discussion above, it is apparent that the development and
application of surrogate modelling approaches for stochastic simulators have
experienced considerable advances in recent years. Nevertheless, there are
some aspects of SK that remain not fully addressed in the literature. Notably,
while most applications of SK in the literature are limited to the use of
constant trend models (Ordinary Kriging), the use of more advanced models for
Universal SK remains largely unexplored in the literature, to the best of the
author’s knowledge. To address this gap, this paper presents an extension of
the deterministic Polynomial-Chaos Kriging (PCK) method originally presented
in [26] for the meta-modelling of stochastic simulators with predictions
contaminated by noise. The PCK meta-model is a powerful technique that
combines PCE with Kriging to develop highly efficient surrogate models for
complex systems. Notably, a prior work by the authors [27] demonstrated the
ability of PCK meta-modelling to successfully capture the behaviour of highly
non-linear functions. Drawing from this positive outcome, the current study
aims to apply the principles of PCK to the realm of stochastic simulators.
Specifically, the proposed meta-model involves a Universal SK model with an
adaptive sparse PCE model as the trend term. The orthonormal polynomial basis
of the PCE is optimally selected by a model selection technique for sparse
linear models, the least-angle regression (LAR) algorithm set out by Efron
[28]. Afterwards, the adjusted LAR-PC model is inserted as the trend term in a
Universal SK model, and the hyper-parameters of the resulting meta-model are
calibrated through a genetic algorithm (GA) global optimization approach. The
effectiveness of the proposed approach is validated through a set of well-
known benchmark case studies in the surrogate modelling literature. The
presented results and discussion demonstrate the potential of the proposed
approach and the importance of selecting a suitable trend model in the context
of SK.

The remainder of this paper is organized as follows. Section 2 outlines the
theoretical formulation of Kriging and SK respectively. Section 3 provides a
concise review of PCE and LAR and illustrates the proposed LAR-PCE SK
approach. Section 4 presents the numerical results and discussion and,
finally, Section 5 closes the paper with concluding remarks and potential
future developments.

##  2 Surrogate modelling: Deterministic and Stochastic Kriging

Let $(\Omega,{\cal F},\mu)$ be a probability space, where $\Omega$ denotes the
event space equipped with a $\sigma$-algebra $\mathcal{F}$ of subsets of
$\Omega$ and the probability measure $\mu$. Given
$\Omega\subset\mathbb{R}^{M}$ as a connected set, let
${\mathcal{M}}:\Omega\subset\mathbb{R}^{M}\to\mathbb{R}$ be a computational
model mapping an $M$-dimensional vector of input (design) variables
$\bm{x}=\left[x_{1},\ldots,x_{M}\right]^{\textrm{T}}$ into the output variable
$y\in\mathbb{R}$ (response surface or QoI). The aim of a surrogate model is to
replace the original model $\mathcal{M}$ by an approximated one
$\mathcal{M}^{*}$, whose evaluation needs lower computational cost. When
defined in a non-intrusive framework, the hyper-parameters of ${\cal M}^{*}$
must be calibrated on a basis of a certain training dataset. Such a training
dataset consists in the realization of the forward model $\cal M$ over a set
ED=$\left\\{\bm{x}^{(1)},\ldots,\bm{x}^{(k)}\right\\}$ of $k$ samples of the
input variable $\bm{x}$ covering the design space, often referred to as the
experimental design (ED). Note that, when the model is deterministic, the
purpose of the surrogate model is simply to approximate the response surface
$y=\mathcal{M}({\bm{x}})$ $\forall\bm{x}\in\Omega$. Conversely, when the model
$\mathcal{M}$ exhibits a stochastic behaviour and its response is affected by
a zero-mean random noise, the goal is to approximate the mean surface
$\mathbb{E}[\mathcal{M}(\Omega)]$ and, possibly, provide a certain
characterization of the uncertainty [17]. In this context, Kriging or
_Gaussian process_ meta-modelling has proved to represent a powerful meta-
modelling approach commonly used in engineering or geostatistics. In a
deterministic computer experiment, the Kriging model assumes that the model
response ${\mathcal{M}}(\bm{x})$ can be conceived as the sum of a
deterministic and a stochastic model as:

$\mathcal{M}\left(\bm{x}\right)\approx\mathcal{M}^{K}\left(\bm{x}\right)=\mathcal{T}\left(\bm{x}\right)+\mathcal{Z}\left(\bm{x}\right),$
(1)

where $\mathcal{Z}\left(\bm{x}\right)$ is a zero-mean stochastic process and
the trend
$\mathcal{T}\left(\bm{x}\right)=\bm{f}\left(\bm{x}\right)^{\textrm{T}}\bm{\beta}=\displaystyle\sum_{r=1}^{s}\beta_{r}f_{r}\left(\bm{x}\right)$
is a linear regression model with $s$ regression coefficients $\beta_{r}$ and
user-selected regression functions $f_{r}\left(\bm{x}\right)$ (e.g
$f_{r}\left(\bm{x}\right)=\bm{x}^{r}$). The formulation in Eq. (1) corresponds
to the most general version of Kriging, often referred to as Universal Kriging
(UK). Simpler expressions known as Simple Kriging (SK) and Ordinary Kriging
(OK) are deduced when the trend term is assumed to be constant throughout the
design space. In SK, it is assumed that the mean of the variable being
predicted is known and constant across the entire area of study, whereas (OK)
estimates the mean based on the available sample data. An interpretation of
Eq. (1) is that deviations from the regression model may resemble a sample
path of a properly chosen stochastic process [29]. The stochastic process
$\mathcal{Z}\left(\bm{x}\right)$ is determined by the process variance
$\sigma$ and an auto-correlation function
$\textbf{R}\left(\bm{x},\,\bm{x}^{\prime}\right)=\mathbf{R}\left(\left|{\bm{x}-\bm{x}^{\prime}}\right|;\bm{\theta}\right)$,
which depends on the spatial distance $\left|\bm{x}^{(i)}-\bm{x}^{(j)}\right|$
and some hyper-parameters $\bm{\theta}$ to be determined. The auto-correlation
function R must satisfy ${\bf{R}}\left(\bm{0};\bm{\theta}\right)=1$ and
${\bf{R}}\left(\bm{x}^{(i)}-\bm{x}^{(j)};\bm{\theta}\right)\to 0$ when
$\left|\bm{x}^{(i)}-\bm{x}^{(j)}\right|\to+\infty$. Thus, ${\cal Z}$ is
characterized by the covariance matrix [30, 31]:

$\bm{\Sigma}_{\mathcal{Z}}:=\left\\{{\rm
Cov}\left({\mathcal{Z}}\left(\bm{x}^{(i)}\right),{\mathcal{Z}}\left(\bm{x}^{(j)}\right)\right)\right\\}_{ij}={\sigma}^{2}\,\bf{R}\left(\left|\bm{x}^{(i)}-\bm{x}^{(j)}\right|;\bm{\theta}\right).$
(2)

Among the possible choices of the function R, a common alternative and the one
used in the present work is the Gaussian correlation function
$\mathbf{R}\left(\left|{\bm{x}-\bm{x}^{\prime}}\right|;\bm{\theta}\right)=\rm{exp}\left(-\bm{\theta}\left|\bm{x}-\bm{x}^{\prime}\right|^{2}\right)$
[29].

The Kriging model is thus fully determined by the hyper-parameters
$\left\\{\bm{\beta},\bm{\theta},\sigma\right\\}$, which are usually estimated
by maximizing the likelihood function over the $k$ samples of the ED [32]:

$\mathcal{L}\left(\bm{\theta},\sigma,{\bm{\beta}};\bf{y}\right)=\frac{1}{\sqrt{\left(2\pi\sigma^{2}\right)^{k}\left|\textbf{R}\right|}}\,\rm{exp}\left[-\frac{1}{2\sigma^{2}}\left(\textbf{y}-\textbf{F}\bm{\beta}\right)^{\textrm{T}}\textbf{R}^{-1}\left(\textbf{y}-\textbf{F}\bm{\beta}\right)\right],$
(3)

where
${\textbf{y}}=\left[{y}^{(1)}={\mathcal{M}}\left(\bm{x}^{(1)}\right),\ldots,{y}^{(k)}={\mathcal{M}}\left(\bm{x}^{(k)}\right)\right]^{\rm
T}$ and $\textbf{F}_{ij}=f_{j}(\bm{x}^{(i)})$ is referred to as the
information matrix. Nevertheless, it is possible to extract estimates in
closed form of the hyper-parameters $\bm{\hat{\beta}}$ and $\bm{\hat{\sigma}}$
as functions of the estimates of $\hat{\bm{\theta}}$ by the best linear
unbiased estimator (BLUE) as:

$\bm{\hat{\beta}}(\hat{\bm{\theta}})=\left(\textbf{F}^{\textrm{T}}\textbf{R}^{-1}\textbf{F}\right)^{-1}\textbf{F}^{\textrm{T}}\textbf{R}^{-1}\textbf{y},\quad\hat{\sigma}^{2}\left(\hat{\bm{\theta}}\right)=\frac{1}{k}\left(\textbf{y}-\textbf{F}{\bm{\beta}}\right)^{\textrm{T}}\textbf{R}^{-1}\left(\textbf{y}-\textbf{F}{\bm{\beta}}\right).$
(4)

Substituting Eq. (4) into (3) and taking logarithm, a simplified optimisation
problem that depends only on ${\bm{\theta}}$ is derived as:

$\hat{\bm{\theta}}=\underset{\bm{\theta}}{\textrm{arg\,min}}\left[\frac{1}{k}\left(\textbf{y}-\textbf{F}{\bm{\beta}}\right)^{\textrm{T}}\textbf{R}^{-1}\left(\textbf{y}-\textbf{F}{\bm{\beta}}\right)\left|\textbf{R}\right|^{1/k}\right].$
(5)

Therefore, the Kriging mean and mean-squared error (MSE) at any point $\bm{x}$
in the parameter space can be computed as [31]:

$\begin{array}[]{c}\hat{\cal
M}^{K}(\bm{x})=\bm{f}(\bm{x})^{\textrm{T}}\hat{\bm{\beta}}+\bm{r}(\bm{x})^{\textrm{T}}\textbf{R}^{-1}\left({\bm{y}}-\textbf{F}\hat{\bm{\beta}}\right),\\\
\widehat{\rm{MSE}}(\bm{x})=\sigma^{2}\left[1-\bm{r}(\bm{x})^{\textrm{T}}\textbf{R}^{-1}\bm{r}(\bm{x})+\bm{u}(\bm{x})^{\textrm{T}}\left(\textbf{F}^{\textrm{T}}\textbf{R}^{-1}\textbf{F}\right)^{-1}\bm{u}(\bm{x})\right],\end{array}$
(6)

where
$\bm{u}\left(\bm{x}\right)=\textbf{F}^{\textrm{T}}\textbf{R}^{-1}\bm{r}(\bm{x})-\bm{f}(\bm{x})$
and
$r_{i}(\bm{x})=\textbf{R}\left(\left|\bm{x}-\bm{x}^{(i)}\right|;\hat{\bm{\theta}}\right)$,
$i=1,\ldots,k$.

To accommodate the existence of stochasticity in the response surface model
${\cal M}(x)$, Ankenman et al. [17] proposed Stochastic Kriging as an
extension of the classical Kriging scheme specified in Eq. (1) to define the
response surface $\mathcal{M}(\bm{x})$ in stochastic terms. In this context,
several simulation runs in the same design point are performed, so the
experimental design must be redefined as
$\textrm{ED}=\left\\{\left(\bm{x}^{(i)},n_{i}\right)\right\\}_{i=1}^{k}$,
where $n_{i}\in\mathbb{N}$ is the number of replications at the design point
$\bm{x}^{(i)}$. In particular, the output of the $l$-th simulation replication
at $\bm{x}^{(i)}$, $i=1,\ldots,k$ is given by:

$\mathcal{M}_{i}\left(\bm{x}^{(i)}\right)={\bm{f}}(\bm{x})^{\textrm{T}}{\bm{\beta}}+\mathcal{Z}\left(\bm{x}^{(i)}\right)+\varepsilon_{l}\left(\bm{x}^{(i)}\right),\quad
l=1,\ldots,n_{i},$ (7)

where the terms $\mathcal{Z}\left(\bm{x}^{(i)}\right)$ and
$\varepsilon_{l}\left(\bm{x}^{(i)}\right)$ represent the model extrinsic
uncertainty and the intrinsic uncertainty for the $l$-th replication,
respectively. In other words, the changes in the model response are due to the
simulator evaluation at different points in the parameter space, or due to the
stochastic nature of the model itself, respectively. Analogously to Eq. (2),
the intrinsic uncertainty can be characterized by the $k\times k$ covariance
matrix:

$\bm{\Sigma}_{\varepsilon}=\left\\{{\rm{Cov}}\left(\bar{\varepsilon}\left({\bm{x}}^{(i)}\right),\bar{\varepsilon}\left({\bm{x}}^{(j)}\right)\right)\right\\}_{ij},\quad\overline{\varepsilon}\left(\bm{x}^{(i)}\right)=\dfrac{1}{n_{i}}\displaystyle\sum_{l=1}^{n_{i}}{\varepsilon_{l}}\left(\bm{x}^{(i)}\right),\,i=1,\ldots,k.$
(8)

In this light, similarly to Kriging, the calibration of the surrogate model
involves maximizing a likelihood function, whose expression is the natural
consequence of incorporating the intrinsic uncertainty in the Kriging
methodology:

$\mathcal{L}\left(\bm{\beta},\sigma^{2},\bm{\theta};\overline{\cal
M}\right)=\frac{1}{\sqrt{\left|\bm{\Sigma}_{\cal
Z}+\bm{\Sigma}_{\varepsilon}\right|\,(2\pi)^{k}}}\,{\rm
exp}\left(-\frac{1}{2}\left(\overline{\cal
M}-\textbf{F}\bm{\beta}\right)^{\textrm{T}}\left(\bm{\Sigma}_{\cal
Z}+\bm{\Sigma}_{\varepsilon}\right)^{-1}\left(\overline{\cal
M}-\textbf{F}\bm{\beta}\right)\right),$ (9)

where

$\overline{\mathcal{M}}=\left[\overline{\mathcal{M}}\left(\bm{x}^{(1)}\right),\ldots,\overline{\mathcal{M}}\left(\bm{x}^{(k)}\right)\right]^{\textrm{T}},\quad\overline{\mathcal{M}}\left(\bm{x}^{(i)}\right)=\frac{1}{n_{i}}\displaystyle\sum_{l=1}^{n_{i}}\mathcal{M}_{l}(\bm{x}^{(i)}),$
(10)

and $\bm{\Sigma}_{\mathcal{Z}}(\bm{x},\cdot)$ represents the $k\times 1$
cross-correlation vector
$\left\\{{\rm{Cov}}\left(\mathcal{Z}\left(\bm{x}\right),\mathcal{Z}\left(\bm{x}^{(j)}\right)\right)\right\\}_{j=1}^{k}$
between the prediction point $\bm{x}$ and each one of the ED points. Upon the
computation of the SK parameters, the SK predictor along with its associated
mean squared error (MSE) can be derived as follows [33]:

$\begin{array}[]{c}\hat{\mathcal{M}}^{SK}\left(\bm{x}\right)={\bm{f}}\left(\bm{x}\right)^{\textrm{T}}\hat{\bm{\beta}}+\bm{\Sigma}_{\mathcal{Z}}\left(\bm{x},\cdot\right)\left[\bm{\Sigma}_{\mathcal{Z}}+\bm{\Sigma}_{\varepsilon}\right]^{-1}\left(\overline{\mathcal{M}}-\textbf{F}\hat{\bm{\beta}}\right),\\\
\widehat{\rm
MSE}\left(\bm{x}\right)=\bm{\Sigma}_{\mathcal{Z}}\left(\bm{x},\bm{x}\right)-\bm{\Sigma}_{\mathcal{Z}}\left(\bm{x},\cdot\right)^{\textrm{T}}\left[\bm{\Sigma}_{\mathcal{Z}}+\bm{\Sigma}_{\varepsilon}\right]^{-1}\bm{\Sigma}_{\mathcal{Z}}\left(\bm{x},\cdot\right)+{\bm{\zeta}}^{\textrm{T}}\left(\textbf{F}^{\textrm{T}}\left[\bm{\Sigma}_{\mathcal{Z}}+\bm{\Sigma}_{\varepsilon}\right]^{-1}\textbf{F}\right)^{-1}{\bm{\zeta}},\end{array}$
(11)

with
$\bm{\zeta}={\bm{f}}\left(\bm{x}\right)-\textbf{F}^{\textrm{T}}\left[\bm{\Sigma}_{\mathcal{Z}}+\bm{\Sigma}_{\varepsilon}\right]^{-1}{\bm{\Sigma}_{\cal
Z}\left(\bm{x},\cdot\right)}$. It is worth mentioning that the last term in
Eq. (11) represents the variability inflation of the SK predictor resulting
from the estimation of $\hat{\bm{\beta}}$. Additionally, it is important to
highlight that, in deterministic scenarios where the intrinsic uncertainty is
negligible, $\bm{\Sigma}_{\varepsilon}$ vanishes and Eq. (11) naturally
converges to Eq. (6). In other words, SK reduces to the standard Kriging
metamodel when $\bm{\Sigma}_{\varepsilon}$ tends to zero.

##  3 LAR-PCE based Stochastic Kriging

The higher complexity of UK compared to OK might suggest that it leads to
improved estimates. However, as reported in the literature, a trend
misspecification can lead UK not to improve the accuracy in its estimates [34]
or even worsen its performance [35]. This limitation can be attributed to the
fact that UK involves the estimation of additional parameters ($\bm{\beta}$),
which can potentially induce over-fitting [36]. Therefore, to harness the
advantages of UK while mitigating its potential drawbacks, it is crucial to
choose a minimal trend basis that effectively captures the complexity of the
response surface. Motivated by these considerations, the objective of this
section is to identify an optimal choice of the trend model $\bm{f}$ in Eq.
(7) by combining PCE and LAR. Hereafter, the fundamentals of PCE and the LAR
algorithm are briefly presented, followed by their combination into Universal
SK.

###  3.1 Polynomial-chaos expansion

The PCE model represents the model output response $y\in\mathbb{R}$ as its
expansion through a basis of orthonormal multivariate polynomials
$\Psi_{\bm{\alpha}}$ as:

$y=\mathcal{M}\left(\bm{x}\right)=\displaystyle\sum_{\bm{\alpha}\,\in\,\mathbb{N}^{M}}a_{\bm{\alpha}}\Psi_{\bm{\alpha}}(\bm{x}),$
(12)

where $a_{\bm{\alpha}}$ are the coefficients of the expansion, and the multi-
dimensional index notation
$\bm{\alpha}=\left[\alpha_{1},\ldots,\alpha_{M}\right]^{\textrm{T}}$,
$\alpha_{i}\in\mathbb{N}$, has been adopted. Depending on the nature of the
input variables, classical families of orthonormal polynomials
$\Psi_{\bm{\alpha}}$ include the Hermite or Legendre polynomials [37].
Although the expression (12) is exact on deterministic simulators for an
infinite number of terms, only a finite number of terms can be computed in
practice. Consequently, different truncation schemes have been proposed in the
literature. One of the simplest approaches consists in selecting the set of
all polynomials whose total degree
$\left|\bm{\alpha}\right|=\displaystyle\sum_{i=1}^{M}\alpha_{i}$ is up to $p$,
that is
$\mathcal{A}^{M,p}=\\{\bm{\alpha}\in\mathbb{N}^{M}:\,0\leq\left|\bm{\alpha}\right|\leq
p\\}$ [38]. This leads to a total number of ${\cal P}$ terms in the expansion:

${\cal P}:=\rm{card}\left(\mathcal{A}^{M,p}\right)=\left(\begin{matrix}M+p\\\
p\\\ \end{matrix}\right).$ (13)

Note however that, for high-dimensional or non-linear problems, this
truncation scheme leads to computing a large number of coefficients. As a
solution, a more restrictive hyperbolic truncation scheme was proposed by
Blatman and Sudret [12]. Based on the common observation in practice that
coefficients associated with high-order interaction terms are usually close to
zero (often referred to as the _sparsity-of-effect principle_), this
truncation scheme selects all the multi-indices belonging to the set
$\mathcal{A}^{M,p,q}=\\{\bm{\alpha}\in\mathbb{N}:\,0\leq\left\|\bm{\alpha}\right\|_{q}\leq
p\\}$ with $q$-norm
$\left\|\bm{\alpha}\right\|_{q}=\left(\displaystyle\sum_{i=1}^{M}\left|\alpha_{i}\right|^{q}\right)^{\frac{1}{q}}$.
Note that, as the $q$-norm decreases, the number of mixed-order terms in the
expansion also decreases until, ultimately leaving only the univariate
polynomials remain.

Once the set of candidate polynomials and an ED of $k$ samples are defined,
the expansion coefficients
${\bm{a}}=\left\\{{\bm{a}_{\alpha}},\;{\bm{\alpha}}\in\mathcal{A}^{M,p}\subset\mathbb{N}^{M}\right\\}$
are obtained by minimizing the expectation of the least squares residuals:

$\bm{a}={\textrm{arg}\min_{\bm{a}\,\in\,\mathbb{R}^{\mathcal{P}}}}\,\frac{1}{k}\displaystyle\sum_{i=1}^{k}\left[{\mathcal{M}}\left(\bm{x}^{(i)}\right)-\displaystyle\sum_{{\bm{\alpha}}\,\in\,{\mathcal{A}}^{M,p}}{\bm{a}}_{\bm{\alpha}}{\bm{\Psi}}_{\bm{\alpha}}\left(\bm{x}^{(i)}\right)\right]^{2}.$
(14)

Collecting the realizations of the QoI in the ED in vector form as
${\textbf{y}}=\left[{y}^{(1)}={\mathcal{M}}(\bm{x}^{(1)}),\ldots,{y}^{(k)}={\mathcal{M}}(\bm{x}^{(k)})\right]^{\rm
T}$, the solution of the optimization problem (14) is given by:

$\begin{array}[]{c}\hat{\bm{a}}=\left({\bm{\Psi}}^{\textrm{T}}{\bm{\Psi}}\right)^{-1}{\bm{\Psi}}^{\textrm{T}}{\bm{y}},\quad{\bm{\Psi}}=\left({\Psi}_{ij}\right)=\left[\Psi_{\bm{\alpha}_{j}}\left(\bm{x}^{(i)}\right)\right]_{\,i=1,\ldots,\,k}^{\,j=1,\ldots,\,{\mathcal{P}}},\end{array}$
(15)

where the information matrix $\bm{\Psi}$ contains the evaluation of the basis
polynomials on the ED. Once Eq. (14) has been solved, the estimates by the
resulting PCE surrogate model read:

${\mathcal{M}}^{PC}(\bm{x})=\displaystyle\sum_{\bm{\alpha}\,\in\,\mathcal{A}^{M,p}}\hat{\bm{a}}_{\bm{\alpha}}{\bm{\Psi}}_{\bm{\alpha}}(\bm{x}).$
(16)

The quality of the resulting expansion is typically assessed through cross-
validation. Among the different metrics available in the literature, the
leave-one-out (LOO) error $\epsilon_{LOO}$ is particularly interesting since
it can be estimated in closed form for linear models. In the case of the PCE
model as given in Eq. (16), $\epsilon_{LOO}$ can be expressed as follows:

$\epsilon_{LOO}=\frac{\displaystyle\sum_{i=1}^{k}\left(\frac{\mathcal{M}\left(\bm{x}^{(i)}\right)-\mathcal{M}^{PC}(\bm{x}^{(i)})}{1-h_{i}}\right)^{2}}{\displaystyle\sum_{i=1}^{k}\left(\mathcal{M}\left(\bm{x}^{(i)}\right)-\hat{\mu}_{\bm{y}}\right)^{2}},$
(17)

where $h_{i}$ is the $i$-th component of the vector
$\mathbf{h}={\rm{diag}}\left\\{\bm{\Psi}\left({\bm{\Psi}}^{\textrm{T}}{\bm{\Psi}}\right)^{-1}{\bm{\Psi}}^{\textrm{T}}\right\\}$
and
$\hat{\mu}_{\bm{y}}=\frac{1}{k}\displaystyle\sum_{i=1}^{k}\mathcal{M}\left(\bm{x}^{(i)}\right)$.

###  3.2 Least angle regression (LAR)

The LAR algorithm is a model selection approach [28] that provides a very
efficient procedure to implement the Least Absolute Shrinkage and Selection
Operator (LASSO) method introduced by Tibshirani [39]. LASSO is a procedure
for shrinking some regression coefficients to zero in a regression model by
adding a restriction or penalty term to the least squares problem as:

$\min_{{\bm{\beta}}\,\in\,{\rm I\kern-1.5ptR}^{\cal
P}}\left(\left\|{\bm{\Psi}}\left(\bm{x}\right){\bm{\beta}}-{\cal
M}\left(\bm{x}\right)\right\|_{2}+\eta\|{\bm{\beta}}\|_{1}\right),\quad\eta\in{\rm
I\kern-1.5ptR},\quad\eta\geq 0.$ (18)

The norm $\left\|\cdot\right\|_{1}$ adds a penalty equal to the absolute value
of the magnitude of the coefficients. In other words, it limits the size of
the coefficients leading to sparse models, i.e. models with few coefficients.
The computation of the solution of the optimisation problem in Eq. (18) is a
convex quadratic programming problem with linear inequality constraints, which
can be solved very efficiently by means of the LAR algorithm [28]. It is based
on a variable selection procedure called Incremental Forward Stepwise
Regression [40] that builds the model sequentially, adding one basis
polynomial function $\Psi_{\bm{\alpha}_{i}}$ at a time to the set ${\cal
A}_{i}$ of active basis functions, i.e. the set of predictors most correlated
with the current residual $r_{i}=y-{\bm{\Psi}}_{{\cal
A}_{i}}\,{\bm{\beta}}_{{\cal A}_{i}}$ at the $i$-th iteration. The LAR
algorithm can be summarily described by the following steps:

  1. 1.

Set the vector residuals ${r}_{0}={y}$ and initialize the coefficients vector
${\bm{\beta}}_{\bm{\alpha}_{i}}=0$, $i=1,\ldots,{\cal P}-1$. Set all the
predictors $\Psi_{\bm{\alpha}_{i}}$ standardized to have mean zero and unit
norm.

  2. 2.

Find the predictor $\Psi_{\bm{\alpha}_{i}}$ most correlated with the current
residual $r_{i}=y-{\bm{\Psi}}_{{\cal A}_{i}}\,{\bm{\beta}}_{{\cal A}_{i}}$
(with ${\bm{\beta}}_{{\cal A}_{i}}$ the vector of coefficients for the set of
the active basis functions ${{\mathcal{A}}_{i}}$ at the $i$-th iteration).

  3. 3.

Move ${\bm{\beta}}_{\bm{\alpha}_{i}}$ from 0 towards its least-squares
coefficient on $\Psi_{\bm{\alpha}_{i}}$, until another predictor
$\Psi_{\bm{\alpha}_{j}}$ has as much correlation with the current residual
$r_{i}$ as $\Psi_{\bm{\alpha}_{i}}$.

  4. 4.

Move $\\{{\bm{\beta}}_{\bm{\alpha}_{i}},\,{\bm{\beta}}_{\bm{\alpha}_{j}}\\}$
in the direction defined by their joint least squares coefficient of the
current residual on $\\{\Psi_{\bm{\alpha}_{i}},\,\Psi_{\bm{\alpha}_{j}}\\}$,
until another predictor $\Psi_{\bm{\alpha}_{l}}$ has the same correlation with
the current residual.

  5. 5.

Continue the process until $m$= min$\,({\cal P},k-1)$ predictors have been
added.

Steps 3 and 4 update the set of active coefficients toward their least-square
value in the form ${\bm{\beta}}_{{\cal A}_{i+1}}={\bm{\beta}}_{{\cal
A}_{i}}+\zeta_{i}\bm{u}_{{\cal A}_{i}}$. The factor $\zeta_{i}$ is the
smallest positive scalar such that some new index $j$ joins the active set
${\cal A}_{i+1}={\cal A}_{i}\cup\\{j\\}$, and the vector $\bm{u}_{{\cal
A}_{i}}$ is known as the descent direction (refer to [28] for further
details). The name least angle arises from a geometric interpretation of this
process, whereby $\bm{u}_{{\cal A}_{i}}$ provides the smallest (and equal)
angle with each of the predictors in ${\cal A}_{i}$. Using the LAR scheme, a
collection of $m$ PCE meta-models of decreasing sparsity is extracted. Then,
the optimal meta-model can be selected by a cross-validation scheme, to which
the LOO error from Eq. (17) has been adopted in this work. Typically, the
optimal model maximizing prediction accuracy and minimizing overfitting
presents a cardinality, i.e. the number of active basis funcions, lower that
the cardinality $\cal P$ of the whole candidate basis. It is known that LAR
only requires ${\mathcal{O}}({\cal P}^{3}+k{\cal P}^{2})$ computations [28].
Therefore, when ${\cal P}\ll k$, it follows that ${\cal P}^{3}<k{\cal P}^{2}$
and ${\mathcal{O}}(k{\cal P}^{2})\sim{\mathcal{O}}(k)$. Thus, the
computational complexity added by the PCE model when inserted as the trend
term in the SK model is marginal with respect to the optimization problem
related to the determination of the hyper-parameters in Eq. (5), which is
${\mathcal{O}}(k^{3})$.

###  3.3 LAR-PCE Stochastic Kriging

Figure 1: Flowchart of the proposed PCE SK surrogate modelling approach for
stochastic simulators.

This section presents the proposed PCE SK surrogate modelling approach,
synthetically described in Fig. 1. One of the most remarkable features of the
proposed model regards the automatic optimal definition of the trend model
$\bm{f}$ for Universal SK. This is achieved by the PCE meta-model introduced
in Eq. (16) using the hyperbolic truncation and the LAR algorithm overviewed
in Sections 3.1 and 3.2, respectively. Note that the resulting approach
provides a trend model with high flexibility in virtue of the LAR algorithm,
in such a way that the terms in the PCE expansion are automatically tuned
depending on the specific problem at hand. For the development of the proposed
meta-model, once the optimal set $\bm{\Psi}$ is identified, a universal SK
surrogate model is fitted (see Eq. (9)) through a GA global optimization
approach to avoid stagnation in local extrema. Then, the prediction of the
response surface and the associated MSE are extracted from Eqs. (25) and (26)
respectively.

Under the common assumption in the literature that
$\left\\{\varepsilon_{l}\left(\bm{x}^{(i)}\right)\right\\}_{l=1}^{n_{i}}$ are
i.i.d. and independent of the random field $\mathcal{Z}$ [17, 25, 41], the
proposed optimal LAR-PCE based SK predictor reads:

$\hat{\mathcal{M}}^{PCESK}\left(\bm{x}\right)={\bm{\Psi}}\left(\bm{x}\right){\bm{\beta}}+\bm{\Sigma}_{\mathcal{Z}}\left(\bm{x},\cdot\right)^{\textrm{T}}\left[\bm{\Sigma}_{\mathcal{Z}}+\bm{\Sigma}_{\varepsilon}\right]^{-1}\left(\overline{{\mathcal{M}}}-\bm{\Psi}{\bm{\beta}}\right),$
(19)

and the corresponding optimal MSE reads

${\rm
MSE}\left(\hat{\mathcal{M}}^{PCESK}\right)=\bm{\Sigma}_{\mathcal{Z}}\left(\bm{x},\bm{x}\right)-\bm{\Sigma}_{\mathcal{Z}}\left(\bm{x},\cdot\right)^{\textrm{T}}\left[\bm{\Sigma}_{\mathcal{Z}}+\bm{\Sigma}_{\varepsilon}\right]^{-1}\bm{\Sigma}_{\mathcal{Z}}\left(\bm{x},\cdot\right).$
(20)

It is straightforward to check that the estimator proposed above is the best
linear estimator of the form $\hat{\cal
M}^{PCESK}\left({\bm{x}}\right)=\lambda_{0}\left({\bm{x}}\right)+\bm{\lambda}^{\textrm{T}}\left({\bm{x}}\right)\overline{\cal
M}$, where $\lambda_{0}$ and $\bm{\lambda}$ are weights depending on $\bm{x}$
and chosen in the sense of minimizing the MSE. Moreover, the estimator is
unbiased. This can be deduced by following the general reasoning given in [42]
as follows:

$\begin{array}[]{c}{\rm{MSE}}={\mathbb{E}}\left[\left({\cal
M}^{K}\left(\bm{x}\right)-\hat{\cal
M}^{PCESK}\left(\bm{x}\right)\right)^{2}\right]=\mathbb{E}\left[\left({\bm{\Psi}}\left(\bm{x}\right){\bm{\beta}}+{\mathcal{Z}}\left(\bm{x}\right)-\lambda_{0}\left(\bm{x}\right)-\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right)\overline{\mathcal{M}}\pm\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right){\bm{\Psi}}{\bm{\beta}}\right)^{2}\right]\\\
\\\
=\left({\bm{\Psi}}\left(\bm{x}\right){\bm{\beta}}-\lambda_{0}\left(\bm{x}\right)-\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right){\bm{\Psi}}{\bm{\beta}}\right)^{2}+{\rm{Var}}\left({\mathcal{Z}}\left(\bm{x}\right)-\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right)\overline{\cal
M}+\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right){\bm{\Psi}}{\bm{\beta}}\right)+\mathbb{E}\left[\left({\mathcal{Z}}\left(\bm{x}\right)-\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right)\overline{\cal
M}+\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right){\bm{\Psi}}{\bm{\beta}}\right)\right]^{2}\\\
\\\
=\left({\bm{\Psi}}\left(\bm{x}\right){\bm{\beta}}-\lambda_{0}\left(\bm{x}\right)-\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right){\bm{\Psi}}{\bm{\beta}}\right)^{2}+\Sigma_{\cal
Z}\left(\bm{x},\bm{x}\right)+\displaystyle\sum_{i,j=1}^{k}\lambda_{i}\left(\bm{x}\right)\lambda_{j}\left(\bm{x}\right){\rm{Cov}}\left({\cal
Z}\left(\bm{x}^{(i)}\right),{\cal
Z}\left(\bm{x}^{(j)}\right)\right)+\displaystyle\sum_{i,j=1}^{k}\lambda_{i}\left(\bm{x}\right)\lambda_{j}\left(\bm{x}\right){\rm{Cov}}\left(\overline{\varepsilon}\left(\bm{x}^{(i)}\right),\overline{\varepsilon}\left(\bm{x}^{(j)}\right)\right)\\\
\\\
-2\displaystyle\sum_{i=1}^{k}\lambda_{i}\left(\bm{x}\right)\lambda_{j}\left(\bm{x}\right){\rm{Cov}}\left({\cal
Z}\left(\bm{x}^{(i)}\right),{\cal
Z}\left(\bm{x}\right))\right)-2\displaystyle\sum_{i=1}^{k}\lambda_{i}\left(\bm{x}\right)\lambda_{j}\left(\bm{x}\right){\rm{Cov}}\left({\cal
Z}\left(\bm{x}^{(i)}\right),\overline{\varepsilon}\left(\bm{x}^{(i)}\right)\right)=\\\
\\\
=\left({\bm{\Psi}}\left(\bm{x}\right){\bm{\beta}}-\lambda_{0}\left(\bm{x}\right)-\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right){\bm{\Psi}}{\bm{\beta}}\right)^{2}+\Sigma_{\cal
Z}\left(\bm{x},\bm{x}\right)+\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right)\left[\Sigma_{\cal
Z}+\Sigma_{\varepsilon}\right]\bm{\lambda}\left(\bm{x}\right)-2\bm{\lambda}^{\textrm{T}}\left(\bm{x}\right)\Sigma_{\cal
Z}({\bm{x}},\cdot),\end{array}$

where $\lambda_{i}\left({\bm{x}}\right)$ denotes the $i$-th component of
$\bm{\lambda}\left({\bm{x}}\right)$. Minimising the MSE with respect to
$\lambda_{0}$ and $\bm{\lambda}$, i.e. taking
$\partial\rm{MSE}/\partial\lambda_{0}=\partial\rm{MSE}/\partial\bm{\lambda}=0$,
extract the solutions
$\lambda_{0}\left({\bm{x}}\right)=\Psi\left({\bm{x}}\right){\bm{\beta}}-\bm{\lambda}^{\textrm{T}}\left({\bm{x}}\right){\bm{\Psi}}{\bm{\beta}}$
and $\bm{\lambda}\left({\bm{x}}\right)=\left[\Sigma_{\cal
Z}+\Sigma_{\varepsilon}\right]^{-1}\Sigma_{\cal Z}\left({\bm{x}},\cdot\right)$
can be readily extracted. With these weights, the optimal linear estimator and
MSE are directly given by Eqs. (19) and (20), respectively. Furthermore,
following that ${\cal Z}$ is a zero-mean stochastic process, i.e.
$\mathbb{E}\left[{\cal Z}\left(\bm{x}\right)\right]=0$, and
$\mathbb{E}\left[\overline{\mathcal{M}}-{\bm{\Psi}}{\bm{\beta}}\right]=0$, it
follows that:

${\mathbb{E}}\left[\hat{\cal M}^{PCESK}\left({\bm{x}}\right)-{\cal
M}^{K}\left(\bm{x}\right)\right]={\bm{\Psi}}\left({\bm{x}}\right){\bm{\beta}}+\Sigma_{\cal
Z}\left({\bm{x}},\cdot\right)^{\textrm{T}}\left[\Sigma_{\cal
Z}+\Sigma_{\varepsilon}\right]^{-1}\mathbb{E}\left[\overline{\mathcal{M}}-{\bm{\Psi}}{\bm{\beta}}\right]-{\bm{\Psi}}\left({\bm{x}}\right){\bm{\beta}}-\mathbb{E}\left[{\cal
Z}\left({\bm{x}}\right)\right]=0,$ (21)

and, therefore, it is extracted that $\hat{\cal M}^{PCESK}\left(\bm{x}\right)$
is an unbiased estimation for ${\cal M}^{K}\left(\bm{x}\right)$. The
estimation of the hyper-parameters of the proposed LAR-PCE SK surrogate model
requires maximizing the likelihood function described in Eq. (9), where the
information matrix $\mathbf{F}$ is simply replaced by the optimal trend matrix
$\bm{\Psi}$. However, some considerations must be made. Firstly, note that the
covariance matrix of the noise in the response surface (intrinsic bias between
the surrogate model and the true simulator’s response),
$\Sigma_{\varepsilon}$, in Eq. (9) is unknown in general. Nonetheless,
$\varepsilon$ is commonly assumed in practice to be normally distributed, so
one can write:

$\left\\{\varepsilon_{l}\left(\bm{x}^{(i)}\right)\right\\}_{l=1}^{n_{i}}\sim{\cal
N}\left(0,V_{i}\right),\quad\mbox{with}\,\,V_{i}\,\,\mbox{ estimated as }\quad
V_{i}=\frac{1}{{n_{i}-1}}\displaystyle\sum_{l=1}^{n_{i}}\left({\cal
M}_{l}\left({\bm{x}}^{(i)}\right)-\overline{\cal
M}_{l}\left({\bm{x}}^{(i)}\right)\right)^{2},$ (22)

leading to a positive definite diagonal covariance matrix
$\Sigma_{\varepsilon}$ (see Theorem 1 in [17]). Secondly, it should be noted
that, when calibrating a Kriging meta-model for deterministic simulators, Eq.
(3) can be simplified as a function of $\bm{\theta}$ by leveraging on the
BLUE, while in SK this fact is no longer valid. Instead, ${\bm{\beta}}$ must
be first approximated as $\hat{\bm{a}}$ in Eq. (15). This entails limited
computational overhead because the trend term has been optimally adjusted
beforehand by LAR-PCE. Then, the extrinsic variance ${\sigma}^{2}$ needs to be
simultaneously optimized with the hyper-parameters ${\bm{\theta}}$ as:

$\hat{\bm{\theta}},\,\hat{{\sigma}}^{2}=\underset{\bm{\theta},\;{\sigma}^{2}}{\textrm{arg\,max}}\,{\rm
log}\left({\cal
L}\,\left({\bm{\beta}},{\sigma}^{2},\bm{\theta}\right)\right).$ (23)

Finally, and based on Eq. (23), the optimal estimate for ${\bm{\beta}}$ as a
function of $\hat{\bm{\theta}}$ and $\hat{\sigma}^{2}$ can be written
similarly to Eq. (4) as:

$\hat{\bm{\beta}}\left(\hat{\bm{\theta}},\hat{{\sigma}}^{2},{\varepsilon}\right)=\left[\bm{\Psi}^{\textrm{T}}\left(\hat{\bm{\Sigma}}_{\cal
Z}+\hat{\bm{\Sigma}}_{\varepsilon}\right)^{-1}\bm{\Psi}\right]^{-1}\bm{\Psi}^{\textrm{T}}\left(\hat{\bm{\Sigma}}_{\cal
Z}+{\hat{\bm{\Sigma}}}_{\varepsilon}\right)^{-1}\overline{\cal M}.$ (24)

Once the optimal parameters have been derived, the optimal predictor can be
re-written as:

$\hat{\mathcal{M}}^{PCESK}\left(\bm{x}\right)={\bm{\Psi}}\left(\bm{x}\right)\hat{\bm{\beta}}+\hat{\bm{\Sigma}}_{\mathcal{Z}}\left(\bm{x},\cdot\right)^{\textrm{T}}\left[\hat{\bm{\Sigma}}_{\mathcal{Z}}+{\hat{\bm{\Sigma}}}_{\varepsilon}\right]^{-1}\left(\overline{{\mathcal{M}}}-\bm{\Psi}\hat{\bm{\beta}}\right).$
(25)

with the MSE estimator

$\widehat{\rm
MSE}\left(\bm{x}\right)=\hat{\bm{\Sigma}}_{\mathcal{Z}}\left(\bm{x},\bm{x}\right)-\hat{\bm{\Sigma}}_{\mathcal{Z}}\left(\bm{x},\cdot\right)^{\textrm{T}}\left[\hat{\bm{\Sigma}}_{\mathcal{Z}}+\hat{\bm{\Sigma}}_{\varepsilon}\right]^{-1}\hat{\bm{\Sigma}}_{\mathcal{Z}}\left(\bm{x},\cdot\right)+{\bm{\gamma}}^{\textrm{T}}\left(\bm{\Psi}^{\textrm{T}}\left[\hat{\bm{\Sigma}}_{\mathcal{Z}}+\hat{\bm{\Sigma}}_{\varepsilon}\right]^{-1}\bm{\Psi}\right)^{-1}{\bm{\gamma}},$
(26)

and
$\bm{\gamma}={\bm{f}}\left(\bm{x}\right)-\bm{\Psi}^{\textrm{T}}\left[\hat{\bm{\Sigma}}_{\mathcal{Z}}+\hat{\bm{\Sigma}}_{\varepsilon}\right]^{-1}{\hat{{\bm{\Sigma}}}_{\cal
Z}\left(\bm{x},\cdot\right)}$. The expression (26) follows from the
development in [42] for standard Kriging by replacing the covariance matrix of
the random field by
$\hat{\bm{\Sigma}}_{\mathcal{Z}}+\hat{\bm{\Sigma}}_{\varepsilon}$.

Let us lastly focus on the common assumption of the stochastic process ${\cal
Z}$ following a Gaussian distribution. In this case, when considering an ED of
$k$ samples, the hypotheses in Eq. (2) imply that $\left({\cal
M}\left(\bm{x}\right),\,\overline{\cal
M}\left({\bm{x}}^{(1)}\right),\ldots,\,\overline{\cal
M}\left({\bm{x}}^{(k)}\right)\right)\sim{\cal
N}_{k+1}\left(\bm{\mu},\bm{\Sigma}\right)$ follow a multivariate Gaussian
distribution with mean vector and covariance matrix given by:

$\bm{\mu}=\left({\cal
M}\left(\bm{x}\right),{\bm{\Psi}}(\bm{x})\,{\bm{\beta}}\right)\,\in\,{\rm
I\kern-1.5ptR}^{k+1}\quad\mbox{and}\quad{\bm{\Sigma}}=\left(\begin{array}[]{cc}{\sigma}^{2}&{\bm{\Sigma}}_{\cal
Z}(\bm{x},\cdot)^{\textrm{T}}\\\ {\bm{\Sigma}}_{\cal
Z}(\bm{x},\cdot)&{\bm{\Sigma}}_{\cal
Z}+{{\bm{\Sigma}}}_{\varepsilon}\end{array}\right)\,\in\,{\rm
I\kern-1.5ptR}^{(k+1)\,\times\,(k+1)},$ (27)

and Eq. (19) coincides with the conditional mean of ${\cal
M}\left(\bm{x}\right)$ given $\overline{\cal M}$, i.e.
$\mathbb{E}\left[\left.{\cal M}\left(\bm{x}\right)\right|\overline{\cal
M}\right]$ (see [30]).

To conclude, the quality of the proposed metamodel must be assessed. In
practical contexts where the computational burden is high, cross-validation
techniques are popular approaches for evaluation the accuracy of stochastic
simulators (refer to refernce [43] for an extensive discussion). Nevertheless,
in this work, to conduct a robust assessment of the proposed PCE SK metamodel,
its accuracy has been evaluated on the basis of a dense validation set ${\rm
VS}:=\\{\bm{x}^{(1)},\ldots,\bm{x}^{(K)}\\}$, with $K\gg k$, covering the
parameter space, and using diverse quality metrics reported hereafter.

##  4 Numerical results and discussion

This section presents an extensive numerical assessment of the proposed LAR-
PCE SK meta-model. The evaluation includes a detailed comparison against
ordinary SK through various well-known benchmark examples from the literature.
These presented numerical studies cover both single (Section 4.1) and multi-
dimensional functions (Sections 4.2 and 4.3) contaminated by noise. For each
case study, given a fixed computational budget, several sampling strategies to
obtain the ED are discussed. These approaches range from (i) using a low
number of sampling points with many replications to represent the uncertainty
to (ii) strategies with numerous sampling points with less budget devoted to
the characterization of the uncertainty. Additionally, intermediate strategies
that fall between these two extremes are also explored. The purpose of
defining different sampling scenarios is to demonstrate that the proposed LAR-
PCE SK method consistently outperforms ordinary SK, regardless of whether a
sparse and deep ED (with few points but high accuracy) or a denser but
shallower ED is chosen. The computational effort allocated at each point is
proportional to its corresponding noise variance in every instance. In all the
analyses, the model performance is assessed by means of two error metrics,
namely the Empirical Root Mean Squared Error (ERMSE) and the Normalized
Maximum Absolute Error (NMAE) over a large VS given by:

${\rm
ERMSE}=\sqrt{\frac{1}{K}\displaystyle\sum_{i=1}^{K}\left(\widehat{\mathcal{M}}\left(\bm{x}^{(i)}\right)-\mathcal{M}\left(\bm{x}^{(i)}\right)\right)^{2}},\quad{\rm
NMAE}=\underset{1\leq i\leq
K}{\max}\frac{\left|\widehat{\mathcal{M}}\left(\bm{x}^{(i)}\right)-\mathcal{M}\left(\bm{x}^{(i)}\right)\right|}{K\sigma_{VS}},\
x^{(i)}\in VS,$ (28)

where:

$\sigma_{VS}=\sqrt{\frac{1}{K}\displaystyle\sum_{i=1}^{K}\left(\overline{\mathcal{M}}\left(\bm{x}^{(i)}\right)-\mathcal{M}\left(\bm{x}^{(i)}\right)\right)^{2}}.$
(29)

Note that the ERMSE metric evaluates the global quality of the meta-model,
while the NMAE metric is an indicator of the local error.

###  4.1 Case study I: M/M/1 Queue

The M/M/1 queue model is a simple yet powerful model used to analyse the
behaviour of a wide range of real-world queuing systems [44], from call
centers to manufacturing lines. The model represents a system where there is a
single server and a line of entities to be served (e.g., packets, customers,
etc.). In this system, entities arrive randomly at a certain rate, called
arrival rate, and wait in a queue until the server becomes available. The
waiting time is determined by the arrival rate and by the service rate, which
refer to the time it takes for the server to process an entity and complete
its service, respectively.

The present case study consists of the simulation of an M/M/1 queue where the
service rate is fixed at one and the arrival rate ranges from $0.3$ to $0.9$.
The response surface $y$ represents the expected waiting time as a function of
the arrival rate $x$. The accuracy of the simulation at each design point is
determined by its time run-length $T$, that is, the period of time over which
the system’s behaviour is observed and analysed. This system is one of the
most studied examples in the SK literature [17, 33, 45]. It is particularly
notable for its distinct and progressively increasing intrinsic variance
across the parameter space. This characteristic naturally provides a
formidable yet simple example to understand the potential benefits of the
proposed approach. Despite the stochastic nature of the model, analytical
expressions for the response surface and its intrinsic variance are known [17,
33]:

$y(x)=\frac{1}{1-x},\ V(x)=\frac{2x\left(1+x\right)}{T\left(1-x\right)^{4}}.$
(30)

In the following, we introduce two variations of the experiment. In the first
one in Section 4.1.1, the intrinsic variance is assumed as known following Eq.
(30), and one single simulation is run for each design point. In the second
variation in Section 4.1.2, the intrinsic variance of the simulator is
estimated by performing multiple replications at each of the ED points. For
both cases, the results obtained by the classical Ordinary SK are compared
with those obtained by the proposed LAR-PCE SK approach.

####  4.1.1 Known intrinsic noise

For this first numerical investigation, we will consider the intrinsic
variance as
$\Sigma_{\varepsilon}={\rm{diag}}\left(V\left(\bm{x}^{(i)}\right)\right)_{i=1}^{k}$.
Additionally, $\overline{\mathcal{M}}$ is directly taken as
$\left(\mathcal{M}\left(x^{(i)}\right)\right)_{i=1}^{k}$ since one single
replication at each design point is considered in this first numerical study.
It is important to note that, although this situation may not reflect the
complexity of most real-world applications, it serves as a useful analysis to
shed light on the characteristics of the proposed LAR-PCE SK approach. Three
different scenarios with varying number of design points and run-lengths but
with constant computational resources have been considered. The first one
involves 10 equispaced design points with a run-length of 6000; the second one
consists of 30 design points with a run-length of 2000; and the third scenario
encompasses 50 design points with a run-length of 1200. In the case of LAR-PCE
SK, following the recommendations provided in reference [46], the maximum
degree of the polynomials is considered as $5$, $10$ and $16$ for the
scenarios one to three, respectively. It is worth noting that, given the one-
dimensional nature of the problem, the LAR algorithm simply selects the
optimal degree for the uni-variate regression model.

Figure 2: Case Study I: M/M/1 example. Comparison between the true simulator
(a,d) against the predictions by the proposed LAR-PCE SK (c,f) and Ordinary SK
(b,e) under the assumption of known variance (one replication per ED sample)
in Scenarios 1 (a,b,c) (10 ED samples and run-length=6000) and 2 (d,e,f) (30
ED samples and run-length=2000). Scatter points denote the ED samples, and
light blue shaded areas represent the 95% confidence interval.

Figures 2 (a,d) depict the true response surface along with the analytically
computed $95\%$ confidence intervals represented as $\left(y(x)\pm
z_{\alpha}\sqrt{V(x)}\right)$ for Scenarios 1 and 2, respectively. Similarly,
Figs. 2 (b,e) and (c,f) illustrate the estimates produced by Ordinary SK and
LAR-PCE SK, respectively. Estimated $95\%$ confidence intervals in Figs. 2
(b,e) and (c,f) have been computed as $\left(\hat{\mathcal{M}}^{SK}\pm
z_{\alpha}\smash{\sqrt{\widehat{MSE}(x)}}\right)$ and
$\left(\hat{\mathcal{M}}^{PCESK}\pm
z_{\alpha}\smash{\sqrt{\widehat{MSE}(x)}}\right)$, respectively. It can be
seen that for low values of $x$, where the intrinsic uncertainty is small,
both meta-models exhibit a similar behaviour. In contrast, for high values of
$x$ where both the QoI and the intrinsic variance experience a rapid growth,
Ordinary SK fails to capture the variation in the response, whereas LAR-PCE SK
exhibits a significantly better performance. This is to be expected since, as
reported in [35], Ordinary SK virtually ignores those simulation outputs whose
estimated intrinsic variance are large compared to the estimated extrinsic
variance. Thus, these results demonstrate that introducing a proper trend
model effectively mitigates this issue.

Figure 3: Performance analysis of Ordinary SK and LAR-PCE SK in Case Study I in terms of ERMS and NMAE error metrics using 100 known intrinsic noise M/M/1 experiments and sampling Scenarios 1 to 3. | SK |  | LAR-PCE SK
---|---|---|---
| Mean ERMSE | Mean NMAE |  | Mean ERMSE | Mean NMAE
Scenario 1 | 0.787 | 0.0076 |  | 0.516 | 0.0049
Scenario 2 | 0.810 | 0.0079 |  | 0.519 | 0.0052
Scenario 3 | 0.832 | 0.0080 |  | 0.650 | 0.0074
Table 1: Summary of the accuracy in the estimates of the M/M/1 model with
known intrinsic noise variance by Ordinary SK and LAR-PCE SK adopting sampling
Scenarios 1 to 3.

To provide a further comparison in performance of LAR-PCE SK with respect to
Ordinary SK, the results of the fitting of one hundred independent experiments
are reported in Fig. 3 and Table 1. Significant enhancements have been
achieved in terms of ERMSE, with improvements ranging from 20% to 35% on
average. Similarly, the NMAE showed average improvements ranging from 8% to
35%. It is noteworthy that, when considering the median values, the
improvements are even more significant, ranging from 28% to 42% for ERMSE and
from 28% to 38% for NMAE. Finally, it can be noted in Fig. 2 that, in terms of
variance estimation, both approaches similarly fail to capture the true
variability. This issue has been recognized and documented in the SK
literature [19, 36]. Although there are methods available to address this
problem (refer e.g. to [36]), they are beyond the scope of our study, in which
the primary focus is on the evaluation of the importance of properly selecting
the trend model in SK.

####  4.1.2 Unknown intrinsic noise

This section introduces a modified version of the previous investigation,
where the intrinsic variance is assumed to be unknown. Consequently, following
Eq. (22), the estimation of the $\Sigma_{\varepsilon}$ matrix is carried out
through a series of replications at each ED point. In this particular
analysis, a total computational budget $C$ of 500 replications to be
distributed along the training sites has been considered. Consistently with
reference [45], the number of replications assigned to each design point
$x_{i}$ is determined in a single step as proportional to the model variance
following
$n_{i}=\frac{\sqrt{V\left(\bm{x}_{i}\right)}}{\sum_{i=1}^{k}\sqrt{V\left(\bm{x}_{i}\right)}}C$.
It should be noted that the use of the actual variance to determine the number
of replications at each point in our analysis has been chosen for the sake of
simplicity. However, alternative methods which employ a multi-step approach
for this purpose such as those described in reference [33] may be considered
for bypassing this issue.

Figure 4: Performance analysis of Ordinary SK and LAR-PCE SK in Case Study I in terms of ERMS and NMAE error metrics using 100 unknown intrinsic noise M/M/1 experiments and sampling Scenarios 1 to 3. | SK |  | LAR-PCE SK
---|---|---|---
| Mean ERMSE | Mean NMAE |  | Mean ERMSE | Mean NMAE
Scenario 1 | 0.073 | 5.16E-04 |  | 0.067 | 4.67E-04
Scenario 2 | 0.107 | 1.26E-03 |  | 0.084 | 7.53E-04
Scenario 3 | 0.142 | 1.80E-03 |  | 0.127 | 1.34E-03
Table 2: Summary of the accuracy in the estimates of the M/M/1 model with
unknown intrinsic noise variance by Ordinary SK and LAR-PCE SK adopting
sampling Scenarios 1 to 3.

The same three scenarios as discussed in Section 4.1.1 have been considered
over a macro-replication of one hundred experiments. The results are
summarized in Table 2 and Fig. 4. As in the previous case, a relevant
improvement in performance can be observed in each situation, ranging from 8%
to 21% on average with respect to ERMSE, and from 9% to 40% with respect to
NMAE. Considering the median values, similar conclusions to the previous
analyses can be extracted, as can be seen in Fig. 4.

###  4.2 Case study II: two dimensional example.

This second case study investigates the egg-box shaped surface suggested in
reference [45]. The simulation output at a design point
${\bm{x}}=\left[x_{1},x_{2}\right]^{\textrm{T}}\in\left[-1,1\right]^{2}$ on
the $j$-th replication is given by:

$\mathcal{M}_{j}\left(x_{1},x_{2}\right)=\sin\left(9x_{1}^{2}\right)+\sin\left(9x_{2}^{2}\right)+\varepsilon_{j}\left(x_{1},x_{2}\right),$
(31)

with $\varepsilon_{j}\sim{\cal
N}\left(0,2+\cos\left(\pi+\left(x_{1}+x_{2}\right)/2\right)\right)$.

Figure 5: Response surface of the Case Study II (a), and predictions by SK
(b), LAR-PCE SK (c) and full PCE SK (d) surrogate models with $k=64$ design
points.

Likewise in the previous case study, three different scenarios have been
considered by varying the number of design points but keeping a fixed
computational budget $C$ of 1280 replications to be distributed among them.
Specifically, the three different settings are characterized by an ED size of
$k$=32, 64, and 128 points, respectively. The locations of the design points
were determined by adopting the Latin Hypercube Sampling (LHS) algorithm, a
commonly used technique for generating diverse and evenly distributed design
points in multi-dimensional space. The number of replications $n_{i}$ at each
design point has been adjusted proportionally to its standard deviation as in
Section 4.1.2. For each scenario, the performance of Ordinary SK is
benchmarked against the proposed LAR-PCE SK meta-model. Additionally, to
stress the importance of choosing an appropriate sparse basis of polynomials
in the trend model, a full PCE SK model without LAR has been also included in
the comparison. In this case, the entire polynomial basis is included up to a
certain degree given by the hyperbolic truncation defined in Section 3.1.
Regarding the polynomial basis settings, a $q$-norm value of 0.8 has been
employed in all cases. With respect to the degree of the polynomial basis to
be considered, in compliance with the recommendations by Sudret [46], the
maximum degree of the polynomial basis for the full PCE SK case has been set
at 5, 8, and 10 for Scenarios 1, 2, and 3, correspondingly. On the other hand,
due to the sparse nature of the basis selected by the LAR algorithm, it is
possible to consider polynomials of slightly higher degree for LAR-PCE SK.
Specifically, degrees of 9, 11, and 12 have been considered in Scenarios 1 to
3, respectively.

Figure 6: Performance analysis of Ordinary SK, LAR-PCE SK and full PCE SK in
terms of ERMS and NMAE error metrics using 100 experiments of Case Study II
and sampling Scenarios 1 to 3.

Figure 5 illustrates the outcomes of an arbitrary replication of Scenario 2.
It depicts the true response surface to be approximated (a), alongside the
meta-models constructed using Ordinary SK (b), LAR-PCE SK (c), and full PCE SK
(d) with a total of $k=64$ design points. Figure 6 and Table 3 present the
results of the comparison between the different approaches discussed in the
considered scenarios over a macro-replication of one hundred independent
experiments. In terms of performance, LAR-PCE SK consistently outperforms the
other methods, while full PCE SK showed the poorest performance based on both
ERMSE and NMAE metrics. Comparing Ordinary SK with the proposed LAR-PCE SK,
average reductions around 50% in ERMSE are found in all the considered
scenarios. With respect to the local accuracy, NMAE enhancements vary from 23%
in Scenarios 1 and 2, to 36% in Scenario 3. Taking into account the median,
similar conclusions can be inferred. On the other hand, it is noteworthy that
blindly setting an inadequate polynomial basis as the trend model resulted in
notable errors and instability in the results across all the scenarios,
highlighting the crucial importance of a proper basis selection. Remarkably,
LAR-PCE SK included a modest number of terms in the basis with a median of 9,
11, and 14 terms for Scenarios 1 to 3, respectively. In contrast, the complete
expansion included 17, 32, and 49 polynomials in the basis for each respective
scenario.

| SK |  | LAR-PCE SK
---|---|---|---
| Mean ERMSE | Mean NMAE |  | Mean ERMSE | Mean NMAE
Scenario 1 | 0.782 | 9.13E-04 |  | 0.385 | 6.98E-04
Scenario 2 | 0.623 | 8.32E-04 |  | 0.316 | 6.42E-04
Scenario 3 | 0.495 | 7.69E-04 |  | 0.228 | 4.90E-04
Table 3: Summary of the accuracy in Case Study II adopting sampling Scenarios
1 to 3.

###  4.3 Case study III: Ishigami function.

This last case study investigates the Ishigami function [26], a classical
benchmark function widely used in the literature. This function is
characterized by its non-linear and non-monotonic behaviour, making it a
challenging test for meta-modelling techniques. The function takes as input a
three-dimensional vector $\bm{x}=\left[x_{1},x_{2},x_{3}\right]^{\textrm{T}}$
with each component $x_{i}$ ranging from $-\pi$ to $\pi$, and outputs a scalar
value $\mathcal{M}_{j}\left(\bm{x}\right)$ given by the non-linear sum of
three sinusoidal functions. For the purpose of this work, a stochastic version
of the Ishigami function is defined by incorporating a normally distributed
random noise term $\varepsilon_{j}\left(\bm{x}\right)$ with zero mean and
standard deviation proportional to the absolute value of the function. This
results in the following stochastic Ishigami function:

$\begin{array}[]{c}\mathcal{M}_{j}\left(x_{1},x_{2},x_{3}\right)=\sin\left(x_{1}\right)+7\sin^{2}\left(x_{2}\right)+0.1x_{3}^{4}\sin\left(x_{1}\right)+\varepsilon_{j}\left(x_{1},x_{2},x_{3}\right),\
x_{1},x_{2},x_{3}\in\left[-\pi,\pi\right],\\\ \varepsilon_{j}\sim{\cal
N}\left(0,\sqrt{\left|\sin\left(x_{1}\right)+7\sin^{2}\left(x_{2}\right)+0.1x_{3}^{4}\sin\left(x_{1}\right)\right|}\right).\end{array}$
(32)

The ED has been defined with a total computational budget of $C=2560$. Three
distinct scenarios considering different ED sizes have been investigated:
$k=64$ (Scenario 1), $k=128$ (Scenario 2), and $k=256$ (Scenario 3). The
allocation of replications has been defined proportionally to the standard
deviation, as outlined in Section 4.2. On this basis, Ordinary SK, LAR-PCE SK,
and full-PCE SK meta-models have been constructed and compared in terms of
accuracy. The calibrations for Scenarios 1 to 3 have been performed with
different maximum degrees for the polynomial basis. Specifically, the full-PCE
SK meta-models have been calibrated with maximum degrees of 4, 6, and 7,
respectively. On the other hand, the LAR-PCE SK meta-models have been
calibrated with slightly higher maximum degrees of 6, 8, and 9 for Scenarios 1
to 3, respectively. Additionally, a value of 0.8 is specified for the $q$-norm
for all the considered scenarios.

Figure 7 shows an arbitrary replication of Scenario 1 ($k=64$ design points),
displaying the true response surface (a) and the predictions of the meta-
models constructed using SK (b), LAR-PCE SK (c), and full PCE SK (d). The
comparison among the different approaches based on a replication of 100
independent experiments is presented in Fig. 6 and Table 3. It can be observed
that LAR-PCE SK consistently outperforms the other methods in terms of both
global (ERMSE) and local (NMAE) accuracy. In this case study, similarly to
Case Study II, full PCE SK exhibits the poorest performance. Nonetheless, in
this particular case study, the difference between full PCE SK and Ordinary SK
is not as pronounced. Instead, when comparing SK with the proposed LAR-PCE SK,
significant enhancements are noticeable in terms of ERMSE, with average
improvements ranging from 35% in Scenario 1 to 74% in Scenario 3. Similar
enhancements are also observed when considering the NMAE metric or the
improvement in median terms. Lastly, it is worth mentioning the sparsity of
the polynomial bases used in the trend model. While the full-PCE expansions
consist of a larger number of terms (23, 50, and 68 in Scenarios 1 to 3), the
LAR-selected polynomials are significantly more sparse, with only 10, 14, and
17 terms on average, respectively. This allows capturing the trend without
encountering over-fitting issues.

Figure 7: Response surfaces of the Ishigami function (a,e), and predictions by Ordinary SK (b,f), LAR-PCE SK (c,g) and full PCE SK (d,h) meta-models with $k=64$ design points. The top and bottom row panels denote the response surfaces obtained for $x_{3}=\pi$ and $x_{1}=\pi/2$, respectively. Figure 8: Performance analysis of Ordinary SK, LAR-PCE SK and full PCE SK in terms of ERMS and NMAE error metrics using 100 experiments of the Ishigami function and sampling Scenarios 1 to 3. | SK |  | LAR-PCE SK
---|---|---|---
| Mean ERMSE | Mean NMAE |  | Mean ERMSE | Mean NMAE
Scenario 1 | 2.169 | 1.77E-03 |  | 1.418 | 1.20E-03
Scenario 2 | 1.421 | 1.46E-03 |  | 0.772 | 7.32E-04
Scenario 3 | 0.983 | 1.11E-03 |  | 0.257 | 2.92E-04
Table 4: Summary of the accuracy in the estimates of the Ishigami function for
sampling Scenarios 1 ($k=64$), 2 ($k=128$), and 3 ($k=256$).

##  5 Conclusions

This study has proposed the use of Universal Stochastic Kriging for
constructing surrogate models for stochastic simulators. Specifically, the
proposed approach, LAR-PCE SK, adopts adaptive sparse PCE as the trend model
in a Universal SK model. With the aim of minimizing over-fitting, the LAR
algorithm is adopted to automatically identify the optimal polynomial basis in
the PCE, providing the resulting model with local/global prediction abilities
and high flexibility for a broad variety of problems. Lastly, the hyper-
parameters of the resulting Universal SK model are fitted through a global GA
optimization. The effectiveness of the proposed approach has been validated
through three classical benchmark case studies, namely: (i) the M/M/1 queue;
(ii) a noisy eggbox-shaped surface; and (iii) a stochastic version of the
Ishigami function. The numerical results and discussion have proved the
appropriateness of the developed scheme to analyse and develop surrogate
models in a wide range of situations, mitigating some of the deficiencies of
Ordinary SK when dealing with stochastic processes with high variability. The
key findings and contributions of this work include the following:

  * 1.

The proposed LAR-PCE SK approach has demonstrated superior performance
compared to Ordinary SK across various scenarios, including both single and
multiple dimensions, with known and unknown intrinsic variance. Benefiting
from the model selection capabilities of the LAR algorithm, LAR-PCE SK has
shown particularly promising for multidimensional applications.

  * 2.

Since the computational complexity added by the PCE trend model in Universal
SK is marginal, the enhancements in accuracy by the proposed LAR-PCE SK do not
entail significant increases in the computational burden compared to Ordinary
SK.

  * 3.

The analysis of Case Studies II and III has revealed that the consideration of
a full PCE as a trend basis leads to unacceptable results due to over-
overfitting. This stresses the need to consider an optimal sparse basis to
effectively enhance the performance of the ordinary SK.

The presented numerical results and discussion suggest the ability of the
proposed methodology to be adapted to a large variety of problems in science
and engineering. In this light, future developments regard the extension of
the presented formulation to other procedures for sparse basis selection.
Another interesting objective for future work involves addressing the inherent
underestimation of the variance within the SK methodology. A well-selected
trend term could have the potential to improve the ability of the stochastic
component in SK to accurately capture the overall stochastic nature of the
model. Thus, developing a strategy to effectively handle this underestimation
could further enhance the performance and robustness of the surrogate
modelling approach.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgement

This work has been partially supported by the Spanish Ministry of Science and
Innovation [PID2020-116809GB-I00], by the Junta de Extremadura (Spain) through
Research Group Grants [GR18023] by the European Regional Development Fund
(ERDF) and by “Junta de Extremadura” (Ref. IB20040). E. García-Macías was
partially supported by the Spanish Ministry of Science and Innovation through
the research project “BRIDGEXT - Life-extension of ageing bridges: Towards a
long-term sustainable Structural Health Monitoring” (Ref.
PID2020-116644RB-I00).

## References

  * [1] L.J. Allen, A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis, Infect. Dis. Model. 2 (2) (2017) 128–142.
  * [2] A. Mourad, F. Mroue, Z. Taha, Stochastic mathematical models for the spread of COVID-19: a novel epidemiological approach, Math. Med. Biol. 39 (1) (2021) 49–76.
  * [3] R. Thiedmann, O. Stenzel, A. Spettl, P.R. Shearing, S.J. Harris, N.P. Brandon, V. Schmidt, Stochastic simulation model for the 3d morphology of composite materials in li–ion batteries, Comput. Mater. Sci. 50 (12) (2011) 3365–3376.
  * [4] S. Meyer, A. Brückner-Foit, A. Möslang, A stochastic simulation model for microcrack initiation in a martensitic steel, Comput. Mater. Sci. 26 (2003) 102–110.
  * [5] G.K. Sakki, I. Tsoukalas, P. Kossieris, C. Makropoulos, A. Efstratiadis, Stochastic simulation–optimization framework for the design and assessment of renewable energy systems under uncertainty, Renew. Sustain. Energy Rev. 168 (2022) 112886.
  * [6] C. Luo, S.L. Zhang, C. Wang, Z. Jiang, A metamodel–assisted evolutionary algorithm for expensive optimization, J. Comput. Appl. Math. 236 (5) (2011) 759–764.
  * [7] [7] K. Khaledi, E. Mahmoudi, M. Datcheva, D. König, T. Schanz, Sensitivity analysis and parameter identification of a time dependent constitutive model for rock salt, J. Comput. Appl. Math. 293 (2016) 128–138.
  * [8] E. García-Macíıas, I. Venanzi, F. Ubertini, Metamodel-based pattern recognition approach for real–time identification of earthquake–induced damage in historic masonry structures, Autom. Constr. 120 (2020) 103389.
  * [9] J. García–Merino, C. Calvo–Jurado, E. García–Macías, Polynomial chaos expansion for uncertainty propagation analysis in numerical homogenization of 2d/3d periodic composite microstructures, Compos. Struct. 300 (2022) 116130.
  * [10] R. Alizadeh, J.K. Allen, F. Mistree, Managing computational complexity using surrogate models: a critical review, Res. Eng. Des. 31 (2020) 275–298.
  * [11] J. Kudela, R. Matousek, Recent advances and applications of surrogate models for finite element method computations: a review, Soft. Computing 26 (24) (2022) 13709–13733.
  * [12] G. Blatman, B. Sudret, Adaptive sparse polynomial chaos expansion based on least angle regression, J. Comput. Phys. 230 (6) (2011) 2345–2367.
  * [13] L. Nuñez, R.G. Regis, K. Varela, Accelerated random search for constrained global optimization assisted by radial basis function surrogates, J. Comput. Appl. Math. 340 (2018) 276–295.
  * [14] J. P. Kleijnen, Kriging metamodeling in simulation: A review, Eur. J. Oper. Res. 192 (3) (2009) 707–716.
  * [15] E. Baker, P. Barbillon, A. Fadikar, R.B. Gramacy, R. Herbei, D. Higdon, J. Huang, L.R. Johnson, P. Ma, A. Mondal, et al., Analyzing stochastic computer models: A review with opportunities, Stat. Sci. 37 (1) (2022) 64–89.
  * [16] V. Moutoussamy, S. Nanty, B. Pauwels, Emulators for stochastic simulation codes, ESAIM: Proceedings and Surveys 48 (2015) 116–155.
  * [17] X. Chen, B.E. Ankenman, B.L. Nelson, Enhancing stochastic kriging metamodels with gradient estimators, Oper. Res. 61 (2) (2013) 512–528.
  * [18] X. Chen, K.K. Kim, Efficient var and cvar measurement via stochastic kriging, Informs. J. Comput. 28 (4) (2016) 629–644.
  * [19] P. Hao, S. Feng, H. Liu, Y. Wang, B. Wang, B. Wang, A novel nested stochastic kriging model for response noise quantification and reliability analysis, Comput. Methods Appl. Mech. Eng. 384 (2021) 113941.
  * [20] X. Ruan, Q. Zhou, L. Shu, J. Hu, L. Cao, Accurate prediction of the weld bead characteristic in laser keyhole welding based on the stochastic kriging model, Metals 8 (7) (2018) 486.
  * [21] J. Pousi, J. Poropudas, K. Virtanen, Game theoretic simulation metamodeling using stochastic kriging, in: Proceedings of the 2010 Winter Simulation Conference, IEEE, 2010, pp. 1456–1467.
  * [22] X. Chen, Q. Zhou, Sequential design strategies for mean response surface metamodeling via stochastic kriging with adaptive exploration and exploitation, Eur. J. Oper. Res. 262 (2) (2017) 575–585.
  * [23] Z.Wang, M. Ierapetritou, Constrained optimization of black-box stochastic systems using a novel feasibility enhanced kriging-based method, Comput. Chem. Eng. 118 (2018) 210–223.
  * [24] M. Binois, R. B. Gramacy, M. Ludkovski, Practical heteroscedastic gaussian process modeling for large simulation experiments, J. Comput. Graph. Stat. 27 (4) (2018) 808–821.
  * [25] Y. Che, Z. Guo, C. Cheng, Generalized polynomial chaos-informed efficient stochastic kriging, J. Comput. Phys. 445 (2021) 110598.
  * [26] R. Schobi, B. Sudret, J. Wiart, Polynomial–chaos-based kriging, Int. J. Uncertain. Quantif. 5 (2) (2015).
  * [27] J. García-Merino, C. Calvo-Jurado, E. Martínez–Pañeda, E. García–Macías, Multielement polynomial chaos kriging–based metamodelling for bayesian inference of nonsmooth systems, Appl. Math. Model. 116 (2023) 510–531.
  * [28] B. Efron, T. Hastie, I. Johnstone, R. Tibshirani, Least angle regression, Ann. Stat. 32 (2) (2004) 407–499.
  * [29] S.N. Lophaven, H.B. Nielsen, J. Sondergaard, et al., DACE: a Matlab Kriging toolbox, Vol. 2, Citeseer, 2002.
  * [30] T.J. Santner, B.J. Williams, W.I. Notz, B.J. Williams, The design and analysis of computer experiments, Vol. 1, Springer, 2003.
  * [31] C.E. Rasmussen, C.K. Williams, et al., Gaussian processes for machine learning, Vol. 1, Springer, 2006.
  * [32] A. Marrel, B. Iooss, F. Van Dorpe, E. Volkova, An efficient methodology for modeling complex computer codes with gaussian processes, Comput. Stat. Data. Anal. 52 (10) (2008) 4731–4744.
  * [33] X. Chen, Q. Zhou, Sequential experimental designs for stochastic kriging, in: Proceedings of the Winter Simulation Conference 2014, IEEE, 2014, 3821–3832.
  * [34] S.D. Tajbakhsh, E. Del Castillo, J.L. Rosenberger, A bayesian approach to sequential optimization based on computer experiments, Qual. Reliab. Eng. Int. 31 (6) (2015) 1001–1012.
  * [35] J. Staum, Better simulation metamodeling: The why, what, and how of stochastic kriging, in: Proceedings of the 2009 Winter Simulation Conference (WSC), IEEE, 2009, 119–133.
  * [36] J.P. Kleijnen, E. Mehdad, Estimating the variance of the predictor in stochastic kriging, Simul. Model. Pract. Theory 66 (2016) 166–173.
  * [37] D. Xiu, G. E. Karniadakis, The wiener–askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2) (2002) 619–644.
  * [38] N. Fajraoui, S. Marelli, B. Sudret, Sequential design of experiment for sparse polynomial chaos expansions, ASA J. Uncertain. Quantif. 5 (1) (2017) 1061–1085.
  * [39] R. Tibshirani, Regression shrinkage and selection via the lasso, J. R. Stat. Soc. Series B 58 (1) (1996) 267– 288.
  * [40] T. Hastie, R. Tibshirani, J.H. Friedman, J.H. Friedman, The elements of statistical learning: data mining, inference, and prediction, Vol. 2, Springer, 2009.
  * [41] X. Chen, B. E. Ankenman, B.L. Nelson, The effects of common random numbers on stochastic kriging metamodels, ACM Trans. Model. Comput. Simul. 22 (2) (2012) 1–20.
  * [42] M.L. Stein, Interpolation of spatial data: some theory for Kriging, Springer Science and Business Media, 1999.
  * [43] J.P.C. Kleijnen, Design and analysis of simulation experiments, Springer, 2018.
  * [44] S. Zheng, A.F. Seila, Some well–behaved estimators for the m/m/1 queue, Oper. Res. Lett. 26 (5) (2000) 231–235.
  * [45] G. Xie, X. Chen, Uniform error bounds for stochastic kriging, in: 2020 Winter Simulation Conference (WSC), IEEE, 2020, 361–372.
  * [46] [46] B. Sudret, Polynomial chaos expansions and stochastic finite–element methods, in: Risk and Rel. Geoth. Eng., CRC Press, 2018, pp. 295–330.