Polynomial Chaos Surrogate Construction for Random Fields with Parametric
Uncertainty

Joy N. Mueller, Khachik Sargsyan, Craig Daniels, and Habib N. Najm

Sandia National Laboratories, Livermore, CA

###### Abstract

Engineering and applied science rely on computational experiments to
rigorously study physical systems. The mathematical models used to probe these
systems are highly complex, and sampling-intensive studies often require
prohibitively many simulations for acceptable accuracy. Surrogate models
provide a means of circumventing the high computational expense of sampling
such complex models. In particular, polynomial chaos expansions (PCEs) have
been successfully used for uncertainty quantification studies of deterministic
models where the dominant source of uncertainty is parametric. We discuss an
extension to conventional PCE surrogate modeling to enable surrogate
construction for stochastic computational models that have intrinsic noise in
addition to parametric uncertainty. We develop a PCE surrogate on a joint
space of intrinsic and parametric uncertainty, enabled by Rosenblatt
transformations, and then extend the construction to random field data via the
Karhunen-Loève expansion. We then take advantage of closed-form solutions for
computing PCE Sobol indices to perform a global sensitivity analysis of the
model which quantifies the intrinsic noise contribution to the overall model
output variance. Additionally, the resulting joint PCE is generative in the
sense that it allows generating random realizations at any input parameter
setting that are statistically approximately equivalent to realizations from
the underlying stochastic model. The method is demonstrated on a chemical
catalysis example model.

Keywords: Stochastic surrogate; polynomial chaos expansion; Karhunen-Loève
expansion; global sensitivity analysis; Sobol indices

##  1 Introduction

Engineering and applied science rely on computational simulations to
rigorously study physical systems. These mathematical models are highly
complex and their simulation parameters are often not known exactly. Usually,
rigorous uncertainty quantification (UQ) studies are sample-intensive, i.e.
they require many evaluations of the model at various conditions and parameter
settings. The space of simulation conditions or parameters (‘inputs’) is
sampled extensively and the model evaluated at each input to obtain a model
response. The set of inputs and associated model responses is then probed, for
example through analysis of variance methods, to explore how the input
uncertainty influences the output.

However, computational models of complex physical phenomena are often
expensive and typically not amenable to intensive sampling schemes. In order
to alleviate this expense, one often constructs a _surrogate model_ of the
input-output map of the underlying physical model. The surrogate is
constructed using a number of training simulations spanning a range of model
inputs. The surrogate model then provides a reasonably good approximation to
the original physical model but the cost of model evaluation is comparatively
negligible. Thus, the surrogate can be queried in place of the original model
in order to perform sample-intensive studies necessary for UQ.

Popular surrogate choices that exist for deterministic models with uncertain
inputs include kriging [33, 36], radial basis functions [12], and support
vector regression [47, 44]. Polynomial chaos expansions (PCEs) have been
particularly successful for surrogate modeling and forward/inverse UQ for a
range of engineering applications [11, 20, 30, 27]. They serve as a convenient
technique for uncertainty propagation, moment estimation, and variance
decomposition due to the orthogonality property of polynomial bases [45].

We note that the outlined forward UQ workflow and the conventional surrogate
options above are for _deterministic_ models with uncertain inputs, where only
one model evaluation per input is necessary to fully characterize the model
response. In contrast, in practical applications many examples arise in which
the model itself is inherently random. That is, uncertainty in the output is
only partially explained by the uncertain input, and there is remaining
uncertainty due to intrinsic noise that is not explicitly modeled. Multiple
evaluations of such _stochastic_ models yield an uncertain response for any
fixed input, and therefore the model response itself is a random variable.
Forward UQ for such stochastic systems requires repeated runs of the model for
every computational condition in order to characterize the model response in
distribution. Such intrinsically stochastic models make the computational
burden of obtaining sufficient numbers of evaluations even more prohibitive.

While local/wavelet-based PCE constructions can be used to represent non-
smooth functions, in most practical applications the effectiveness of (global)
PCEs relies on the smoothness of the input-output map comprising the model. As
such, they are not well suited for noisy/non-smooth input-output maps. As a
remedy, PCEs constructed for noisy systems have typically been built for
integrated output quantities of interest (QoIs), e.g. expectations of model
output [25, 4] or quantiles of interest [19, 32, 46].

In contrast, surrogates for stochastic models can be devised using a number of
alternate approaches. For example, random processes can be optimally
represented using the Karhunen-Loève expansions (KLE) [18, 11, 24, 41], a
generative construction that captures the first and second order two-point
statistics of a random process. Alternatively, distribution-based
constructions are available. For example, the Rosenblatt transformation can be
used to provide a mapping between random vectors and associated standard
random variables via conditional cumulative distribution functions (CDFs)
[35]. Further, optimal transport theory and transport maps have been used
extensively to construct mappings between probability density functions (PDFs)
[26]. Analogously, methods which learn a conditional target density from data
such as normalizing flows [34] have found use in machine learning.

It is observed in [3] that recent surrogate constructions designed to predict
the full response PDF, rather than only the output, often have difficulties
producing uncertainties and identifying critical inputs. This limits their
usefulness in a UQ context.

In this work, we generally strive to delay the averaging and consequent loss
of information until after the surrogate is constructed. That is, we would
like to develop a surrogate approximation to the stochastic model that can
generate replica samples with requisite statistics/distributions. We further
seek a surrogate model construction that provides a natural connection between
the model inputs and these output realizations, while also allowing
straightforward decomposition of uncertainties.

Broadly, there are two approaches for achieving this. The first utilizes a
‘multi-stage’ design that captures the influence of the intrinsic noise on the
QoI by extracting a set of summary statistics for the distribution of noisy
outputs. Methods described above are then used to directly map model inputs to
these summary statistics [3, 28, 49]. Typically, a set of sample replica data
is necessary to provide information about the noise distribution.

The second approach more directly fits the model response to a distribution
using either parametric or non-parametric strategies. Non-parametric methods
such a kernel density estimation (KDE) may be used to represent arbitrary
types of response distributions [13, 51] without imposing a priori assumptions
on the form of the PDF. However, non-parametric methods suffer from the curse
of dimensionality, potentially requiring many sample replicas when
representing stochastic models. Parametric methods are more restrictive,
fitting the data to known distribution families, but usually require fewer
samples. Notably, [23] and [49, 51] attempt to balance this by fitting the
data to the highly expressive parametrized class of generalized-lambda
distributions.

We note that several of these approaches (see [28, 49, 23, 51], for example)
utilize PCEs in conjunction with other methods to construct a surrogate model,
resulting in stochastic surrogates which are not PCEs themselves. Similar to
[29], in this work we remain in the PCE framework by employing the Rosenblatt
transformation [35] to incorporate the generative, stochastic dimensions while
keeping the surrogate approximation valid in the parametric dimensions.

A key element of UQ is assessing the sensitivity of the model output to the
uncertainty of the input [38]. While local approaches for sensitivity analysis
exist, the validity of these methods over a range of values of uncertain
variables is limited. Global sensitivity analysis (GSA) methods are required
to accurately represent sensitivities when input uncertainties are large, the
model is non-linear, and/or interactions exist between model inputs [39].
Particularly, variance-based GSA methods exist, most commonly relying on Sobol
sensitivity indices, which decompose the output variance and identify
influential parameters by apportioning output variance to each uncertain
input.

In most cases, Sobol indices cannot be found analytically and are estimated
via sampling methods. However, global sensitivity indices computed from PCEs
are analytically available and allow indices to be computed without sampling
error [41, 40, 45]. By remaining in the PCE framework, we have immediate
access to analytical expressions for global sensitivity indices and consequent
variance decomposition into parametric uncertainty and intrinsic noise
contributions. That is, while developed in the deterministic input-output map
framework, Sobol indices can be extended to stochastic systems in natural
ways. As pointed out in [50], the treatment of the intrinsic noise alters the
definition of the Sobol index. Typically, estimation of sensitivity indices is
approached from one of three common viewpoints: (i) the intrinsic noise is
represented with latent (hidden) variables that are included with the set of
uncertain model inputs [16]; (ii) noise is eschewed by defining integrated
QoIs and constructing the surrogate for these; or (iii) the indices are
defined as functions of the latent (hidden) variables and are random
quantities [14, 17]. Sensitivity indices calculated for our PCE surrogate,
constructed on a joint stochastic-parametric space, can be reformulated under
each of these definitions and are useful for gauging the contribution of
intrinsic noise to output variability.

The following text is organized as follows. In Section 2, we describe the
basic building blocks of the methodology, PCEs and Rosenblatt maps, and the
joint PCE construction for vector QoIs. We then introduce the Karhunen-Loève
expansion (KLE) for random field quantities in Section 3, and briefly show how
the joint PCE extends to field QoIs in a natural way. Next, Section 4
describes the application problem of interest and results of our method
applied to it. We close with concluding remarks in Section 5.

##  2 Methodology

Consider a stochastic model $\bm{y}=f(\bm{\lambda},\omega)$ where the forward
function $f$ depends on a set of input parameters
$\bm{\lambda}\in\mathbb{R}^{\tilde{d}}$ and outputs a random vector
$\bm{y}\in\mathbb{R}^{d}$. With $\omega$ we denote an element of the sample
space explicitly indicating the stochasticity of the model. The underlying
assumption is that $f$ is expensive to evaluate directly, and therefore it is
useful to identify a simpler model whose input-output behavior closely mimics
the underlying forward function. The simple model acts as a surrogate
representation for $f$ and can be used in its place to perform sample-
intensive studies of the system of interest such as uncertainty propagation,
sensitivity analysis, or parameter estimation.

###  2.1 Polynomial Chaos Expansions

Polynomial chaos expansions are functional representations of random variables
and provide a convenient framework for UQ in computational models [11, 21, 7,
48, 30, 9]. Typically, uncertain inputs are written as a polynomial expansion

$\displaystyle\bm{\lambda}\simeq\sum_{p=0}^{P-1}\bm{a}_{p}\Psi_{p}(\bm{\xi})$
(1)

in terms of multivariate polynomials $\Psi_{p}(\bm{\xi})$, where
$\bm{\xi}=(\xi_{1},\dots,\xi_{\tilde{d}})$ is a vector of independent
identically distributed _i.i.d._ standard random variables, and where the
$\Psi_{p}(\bm{\xi})$ polynomials are orthogonal with respect to the density of
$\bm{\xi}$. The multivariate polynomials in turn are written as products of
univariate polynomials
$\Psi_{p}(\bm{\xi})=\prod_{\ell=1}^{\tilde{d}}\psi_{p_{\ell}}(\xi_{\ell})$,
where each scalar index $p=p(p_{1},\dots,p_{\tilde{d}})$ accounts for the
corresponding _multi-index_ or multivariate degree
$(p_{1},\dots,p_{\tilde{d}})$.

We employ Legendre-Uniform PCEs (Legendre polynomials as functions of standard
uniform random variables), given that we have uncertain input parameters with
compact support and are targeting uniform point accuracy of the surrogate.
Moreover, we assume a first-order, independent-component PCE for each
parameter $\lambda_{i}=\frac{b_{i}+a_{i}}{2}+\frac{b_{i}-a_{i}}{2}\xi_{i}$
given a physically meaningful interval of uncertainty $[a_{i},b_{i}]$, for
$i=1,\dots,\tilde{d}$. This simplifying assumption is made without loss of
generality: we largely rely on a correspondence
$\bm{\lambda}=\bm{\lambda}(\bm{\xi})$ between the physical parameter
$\bm{\lambda}$ and the parametric _germ_
$\bm{\xi}\in\text{Unif}[-1,1]^{\tilde{d}}$.

Conventionally, for a _deterministic_ model $g(\bm{\lambda})$, one writes an
associated PCE representation

$\displaystyle
g(\bm{\lambda})\simeq\sum_{p=0}^{P-1}\bm{c}_{p}\Psi_{p}(\bm{\xi})$ (2)

with the same germ $\bm{\xi}$, and finds coefficients $\bm{c}_{p}$ by
regression using a set of $N$ model evaluations $g(\bm{\lambda^{(n)}})$ at
select training inputs $\bm{\lambda^{(n)}}=\bm{\lambda(\bm{\xi^{(n)}})}$, for
$n=1,\dots,N$. Together with the map $\bm{\lambda}(\bm{\xi})$, the PCE in Eq.
(2) acts as a surrogate representation for the deterministic model
$\bm{y}=g(\bm{\lambda})$ [40].

Representing random variables with PCEs offers a number of advantages. Due to
the orthogonality of the basis functions, one can show that the mean of the
random variable is the zeroth-order coefficient, and its variance is a
weighted sum-of-squares of the rest of the coefficients:

$\displaystyle\mathbb{E}[\bm{y}]=c_{0},\qquad\mathbb{V}[\bm{y}]=\sum_{p=1}^{P-1}c_{p}^{2}||\Psi_{p}||^{2},$
(3)

where the norm $\|\Psi_{p}\|$ is defined as weighted with the density of
$\bm{\xi}$, i.e.

$\displaystyle\|\Psi_{p}\|^{2}:=\int\Psi_{p}(\bm{\xi})^{2}p_{\bm{\xi}}(\bm{\xi})\mathrm{d}\bm{\xi}.$
(4)

Furthermore, PCE-based surrogate models enable variance-based decomposition
and subsequent extraction of Sobol indices. The Sobol indices decompose the
output variance into fractional variance contributions of each parameter or
group of parameters [37]. Here we examine the main-effect (first-order)
sensitivities, which measure the influence of each parameter in isolation on
the output variance. These indices are defined as variance ratios and can be
computed exactly from the PCE coefficients:

$\displaystyle S_{i}$
$\displaystyle:=\frac{\mathbb{V}_{\lambda_{i}}[\mathbb{E}_{\lambda_{-i}}[g(\bm{\lambda}\,|\,\lambda_{i})]]}{\mathbb{V}[\bm{y}]}=\frac{1}{\mathbb{V}[\bm{y}]}\sum_{p\in\mathcal{P}_{S_{i}}}c_{p}^{2}||\Psi_{p}||^{2},$
(5) $\displaystyle\text{where }\mathcal{P}_{S_{i}}:=\\{p:p_{i}(p)>0,\
p_{j}(p)=0\text{ for }j\neq i\\}$

where $\mathbb{V}_{\lambda_{i}}$ denotes the variance with respect to the
$i$th parameter and $\mathbb{E}_{\lambda_{-i}}$ is the expectation with
respect to all remaining parameters. Similar to the first two moments in Eq.
(3), the right-hand-side of Eq. (5) is derived using the orthogonality of PCE
bases [5, 40, 45].

###  2.2 The Rosenblatt Map and Stochastic PCEs

Underpinning the PCE surrogate construction for the deterministic forward
model described above is the explicit mapping between uncertain parameters
$\bm{\lambda}$ – and consequently, model outputs $g(\bm{\lambda})$ – and the
PCE germ $\bm{\xi}$. In this work, we refer to these PCE surrogates as
‘parametric PCEs.’

Now, suppose that our model has intrinsic noise, i.e. for each _fixed_ value
of $\bm{\lambda}$ its output is a random quantity $f(\bm{\lambda})$. For ease
of presentation, we write the forward function as $f(\bm{\lambda},\omega)$ to
explicitly highlight a sample space element $\omega$. Surrogate construction
is not as straightforward as described in Section 2.1, since there is no
explicit control on $\omega$. Having said that, one can construct a PCE
representation of a random variable
$\bm{y}=f(\bm{\lambda},\omega)\in\mathbb{R}^{d}$ given $M$ samples
$\\{\bm{y}^{(m)}=f(\bm{\lambda},\omega^{(m)})\\}_{m=1}^{M}$ relying on the
Rosenblatt transformation [35] and kernel density estimation (KDE) [43]. The
Rosenblatt transformation is the multivariate generalization of the cumulative
distribution function (CDF) theorem,

$\displaystyle\zeta_{1}$ $\displaystyle=\mathcal{F}_{1}(y_{1})$
$\displaystyle\zeta_{2}$ $\displaystyle=\mathcal{F}_{2|1}(y_{2}|y_{1})$
$\displaystyle\zeta_{3}$
$\displaystyle=\mathcal{F}_{3|2,1}(y_{3}|y_{2},y_{1})$ (6)
$\displaystyle\hskip 4.97922pt\vdots$ $\displaystyle\zeta_{d}$
$\displaystyle=\mathcal{F}_{d|d-1,\dots,2,1}(y_{d}|y_{d-1},\dots,y_{2},y_{1}),$

mapping any $d$-dimensional random vector $\bm{y}$ to the $d$-dimensional
random vector $\bm{\zeta}$ composed of _i.i.d._ uniform random variables
$(\zeta_{1},\dots,\zeta_{d})$. In fact, the inverse of the Rosenblatt
transformation $R_{\bm{\lambda}}^{-1}(\bm{\zeta})=\bm{y}$ essentially provides
a mapping from $\omega$ to the germ $\bm{\zeta}$ (see [41] for details). Here
we make the modeling choice of representing the possible
$\bm{\zeta}({\bm{\lambda}})$ with a single underlying germ $\bm{\zeta}$ for
all ${\bm{\lambda}}$ that has stochastic dimension $d$. Essentially, this
pushes all dependence on ${\bm{\lambda}}$ into the expansion coefficient which
we denote ${\bm{z}}_{s}({\bm{\lambda}})$. In other words,
$R_{\bm{\lambda}}^{-1}$ acts as a forward map allowing construction of the PCE
similar to Section 2.1. Consequently, the PCE is written, for any
$\bm{\lambda}$, as

$\displaystyle\bm{y}={R}_{\bm{\lambda}}^{-1}(\bm{\zeta})\simeq\sum_{s=0}^{S-1}\bm{z}_{s}({\bm{\lambda}})\Psi_{s}(\bm{\zeta}).$
(7)

This expansion is built purely on samples of the intrinsic noise for a fixed
parameter value ${\bm{\lambda}}$ and, consequently, we refer to it as a
‘stochastic PCE.’

###  2.3 The joint PCE Surrogate

With the parametric and stochastic PCE constructions in mind, we now describe
the algorithm (see also [29]) to arrive at a PCE surrogate on a joint space
$(\bm{\lambda},\omega)$ of parametric and stochastic dimensions.

Suppose we have an ensemble of parameter configurations
$\bm{\lambda}^{(n)}\,(n=1,\dots,N)$ and, for each fixed $\bm{\lambda}^{(n)}$,
we obtain $M$ samples of the stochastic model
$\\{\bm{y}^{(n,m)}=f(\bm{\lambda}^{(n)},\omega^{(n,m)})\\}_{m=1}^{M}$ where
$n$ is fixed. For each parametric sample, we build a multi-output stochastic
PCE according to Eq. (7) which takes the form

$\displaystyle
f(\bm{\lambda}^{(n)},\omega)\simeq\sum_{s=0}^{S-1}\bm{z}_{s}({\bm{\lambda}}^{(n)})\Psi_{s}(\bm{\zeta}),\quad
n=1,\dots,N.$ (8)

The next step in our joint PCE construction is to build a parametric PCE for
the coefficients $\bm{z}_{s}({\bm{\lambda}})$, $s=0,\dots,S-1$, where the
samples
$\\{{\bm{\lambda}}^{(n)},\bm{z}_{s}({\bm{\lambda}}^{(n)})\\}_{n=1}^{N}$ will
serve as training data for the parametric PCE regression. The PCE then takes
the form

$\displaystyle\bm{z}_{s}(\bm{\lambda})\approx\bm{z}_{s}^{\text{\tiny
PCE}}(\bm{\lambda})\simeq\sum_{p=0}^{P-1}\bm{a}_{sp}\Psi_{p}(\bm{\xi})$ (9)

Replacing the coefficients in Eq. (8) with this coefficient surrogate Eq. (9),
we obtain a joint stochastic-parametric PCE representation

$\displaystyle
f(\bm{\lambda},\omega)\simeq\sum_{s=0}^{S-1}\bigg{(}\sum_{p=0}^{P-1}\bm{a}_{sp}\Psi_{p}(\bm{\xi})\bigg{)}\Psi_{s}(\bm{\zeta}).$
(10)

In other words, we arrive at a joint PCE by building regression-based
stochastic PCEs from inverse-Rosenblatt maps for each parameter $\bm{\lambda}$
followed by another regression with respect to the parameter $\bm{\lambda}$
itself. Note that one can rewrite the nested expansion Eq. (10) as a single
joint expansion with a higher-dimensional multivariate basis

$\displaystyle
f(\bm{\lambda},\omega)\simeq\sum_{j=0}^{J-1}\bm{c}_{j}\Psi_{j}(\bm{\xi},\bm{\zeta})$
(11)

where the PCE germs $\bm{\xi}$ and $\bm{\zeta}$ correspond to their ‘physical’
counterparts $\bm{\lambda}$ and $\omega$, respectively.

The mixed PCE in Eq. (11) offers a number of advantages. It can be viewed as a
generative model while remaining statistically consistent with the underlying
model for each parameter value. Besides, it allows additional flexibility of
employing mixed bases. That is, while bounded parameter ranges naturally lead
to LU bases in the $\bm{\xi}$ dimensions, one can employ Gauss-Hermite bases
in the $\bm{\zeta}$-dimensions corresponding to intrinsic stochasticity as the
latter is often expected to have Gaussian behavior for many physical problems
of interest due to aggregation of multiple independent contributions to
intrinsic noise.

Besides, the joint PCE Eq. (11) is flexible enough to allow various
interpretations of GSA for stochastic models. Namely, as categorized in [50],
the variance-based GSA for stochastic models can be approached from one of
three viewpoints: (i) the intrinsic noise is summarized with a (set of) latent
variable(s) that are included with model inputs and the Sobol indices are
functions of these variables [16]; (ii) the intrinsic noise is first
eliminated by averaging processes which represent the model response with one
or more summary statistics (QoIs) for the system and the Sobol indices are
functions of the QoIs; or (iii) the intrinsic noise is associated with
(deterministic) realizations of the model obtained by holding the random seed
fixed, and obtaining Sobol indices as functions of the random seed [14].

Our treatment of the intrinsic noise and resulting Sobol indices most closely
aligns with viewpoint (i), although it is flexible enough to be framed by
viewpoints (ii) and (iii) also. That is, framing the intrinsic noise as a
random variable $\omega$ which we augment to the inputs
$({\bm{\lambda}},\omega)$ is analogous to viewpoint (i). We can extract
variance-based global sensitivity indices using Eq. (5) and attribute the
variance in the output to either uncertainties in the parameters or intrinsic
noise from the stochastic model. Specifically, the Sobol index from the
intrinsic model noise is essentially a _group sensitivity_ that is found by
summing the effects of all $\bm{\zeta}$.

Furthermore, the QoI-based framework (ii) applies also if we treat the
Rosenblatt coefficients from the stochastic PCE construction as summary QoIs.
For example, in the first-order case, the Rosenblatt-PCE coefficients are
identically the mean and the standard deviation, which are summary statistics
for the system. The parametric PCE is constructed for these QoIs, effectively
averaging out the intrinsic noise.

Finally, if we rearrange the sums in Eq. (10)

$\displaystyle
f(\bm{\lambda},\omega)\simeq\sum_{p=0}^{P-1}\bigg{(}\sum_{s=0}^{S-1}\bm{a}_{sp}\Psi_{s}(\bm{\zeta})\bigg{)}\Psi_{p}(\bm{\xi})$

we obtain parametric Sobol indices that are dependent on $\bm{\zeta}$, which
is exactly the form of the indices obtained in viewpoint (iii).

##  3 Representing Random Fields

So far, the stochastic models we have presented are defined over a joint space
$({\bm{\lambda}},\omega)$ of parametric and stochastic dimensions, and the
_deterministic coordinates_ ${\bm{x}}$ of the system – such as the spatial
domain and time evolution – are fixed. This gives a single ‘snapshot’ of the
system at a given location and an instant in the model’s evolution. To study
the system progression across a continuous range of deterministic coordinates,
we now consider stochastic processes
$\textbf{Y}=f({\bm{\lambda}},\omega;{\bm{x}})$ where the forward function $f$
depends explicitly on the system deterministic coordinates
${\bm{x}}\in\mathbb{R}^{L_{\bm{x}}}$ and outputs a random field QoI
$\textbf{Y}\in\mathbb{R}^{d\times L_{\bm{x}}}$. With ${\bm{\lambda}}$ and
$\omega$ we denote the input parameters and model stochasticity, as before.
The dimension $L_{\bm{x}}$ encapsulates any combination of spatial and/or
temporal dimensions such as $N_{\bm{x}}\times N_{\bm{y}}$ or $N_{\bm{x}}\times
N_{\bm{t}}$ or $(N_{\bm{x}}\times N_{\bm{y}})\times N_{\bm{t}}$.

We rely on Karhunen-Loève expansions (KLE) to arrive at a compact
representation of the random field [18, 11, 10, 15, 41]. The algorithm
described in Section 2 is then easily applied to the expansion leading to a
joint PCE surrogate that can be queried as before for sample-intensive UQ
studies of the system.

###  3.1 Compact representation of fields with KLE

The Karhunen-Loève expansion is a well-established spectral expansion
technique that expresses a random field as weighted sums of zero-mean
orthonormal random variables,

$\displaystyle f({\bm{\lambda}},\omega;{\bm{x}})$
$\displaystyle=f_{0}({\bm{x}})+\sum_{l=1}^{L_{\bm{x}}}\eta_{l}({\bm{\lambda}},\omega)\sqrt{\mu_{l}}\,\phi_{l}({\bm{x}}),\qquad
L_{\bm{x}}<\infty.$ (12)

Here, the expectation
$f_{0}({\bm{x}}):=\mathbb{E}_{\lambda,\omega}[f({\bm{\lambda}},\omega;{\bm{x}})]$
centers the process. The $(\mu_{l},\phi_{l}({\bm{x}}))$ are (deterministic)
eigenvalue-eigenfunction pairs of the sample covariance

$\displaystyle\mathbb{C}({\bm{x}},{\bm{x}}^{\prime})=\mathbb{E}_{\lambda,\omega}[(f({\bm{\lambda}},\omega;{\bm{x}})-f_{0}({\bm{x}}))(f({\bm{\lambda}},\omega;{\bm{x}}^{\prime})-f_{0}({\bm{x}}^{\prime}))]$

and the eigenvalues are ordered $\mu_{1}\geq\mu_{2}\geq\dots\geq 0$. The
random coefficients $\eta_{l}({\bm{\lambda}},\omega)$ are uncorrelated with
zero mean and unit variance [18, 11, 10, 15, 41, 23].

While the expansion is infinite by definition, in a discretized $\bm{x}$
context it becomes exact once the number of included eigenterms $L_{\bm{x}}$
equals the number of deterministic coordinates, or the length of the
discretized $\bm{x}$ vector. In practice, the sum in Eq. (12) is truncated to
obtain a lower-dimensional representation of the random field

$\displaystyle f({\bm{\lambda}},\omega;{\bm{x}})$ $\displaystyle\approx
f_{0}({\bm{x}})+\sum_{l=1}^{L}\eta_{l}({\bm{\lambda}},\omega)\sqrt{\mu_{l}}\,\phi_{l}({\bm{x}}),\quad
L\ll L_{\bm{x}}.$ (13)

Identification of a truncation point $L$ exploits simple expressions for the
model variance. First, the pointwise model variance is computed from the
expansion at each ${\bm{x}}$ by taking advantage of independence of $\eta_{l}$
and orthonormality of $\phi_{l}({\bm{x}})$.

$\displaystyle\mathbb{V}_{\bm{\lambda},\omega}[f({\bm{x}})]=\sum_{l=1}^{L_{\bm{x}}}\mu_{l}\phi_{l}({\bm{x}})^{2}$
(14)

Integrating the pointwise variance over the set of deterministic coordinates
${\bm{x}}\in\Omega$ gives the total output variance as the sum of the
eigenvalues.

$\displaystyle\mathbb{V}_{\bm{x}}[f]=\int_{\Omega}\sum_{l=1}^{L_{\bm{x}}}\mu_{l}\phi_{l}({\bm{x}})^{2}d{\bm{x}}=\sum_{l=1}^{L_{\bm{x}}}\mu_{l}$
(15)

We utilize the natural ordering of the eigenvalues to obtain the fraction of
the total output variance explained by the first $L$ terms of the expansion.

$\displaystyle\frac{\sum_{l=1}^{L}\mu_{l}}{\sum_{l=1}^{L_{\bm{x}}}\mu_{l}}=1-\varepsilon,\quad
0<\varepsilon<1$ (16)

From this, the truncation point is chosen so that $\varepsilon$ is negligibly
small, for example $\varepsilon<0.01$ or at least $99\%$ of the explained
variance is preserved by the KLE.

Note that the number of terms needed in the KLE depends intimately on the rate
of decay of the eigenvalues. If the field is smooth, then the spectrum will
rapidly decay and only the first few expansion terms are needed to represent
the random field with reasonable accuracy. On the other hand, spectral decay
is slow and a large number of expansion terms are needed to represent rough
fields (e.g. Gaussian white noise).

###  3.2 Joint PCE surrogates for random fields

We now extend the joint PCE surrogate construction to random field QoIs
$\textbf{Y}=f({\bm{\lambda}},\omega;{\bm{x}})\in\mathbb{R}^{d\times
L_{\bm{x}}}$. Suppose we have a collection of parameter configurations
${\bm{\lambda}}^{(n)},n=1,\dots,N$, and for each fixed ${\bm{\lambda}}^{(n)}$
we obtain $M$ samples of the random field
$f({\bm{\lambda}}^{(n)},\omega^{(n,m)};{\bm{x}})$. From this, we wish to
construct a joint PCE surrogate in the form of Eq. (11) using the methods
described in Section 2.

We first find a low dimensional representation via the KLE for the random
field QoI(s) using Eq. (13). The expansion is _global_ in the sense that the
set of _combined_ stochastic-parametric samples are taken together in this
expansion. The expansion coefficients $\eta_{l}({\bm{\lambda}},\omega)$ are
random variables and will be represented with joint PCE surrogates:

$\displaystyle\eta_{l}({\bm{\lambda}},\omega)\simeq\sum_{j=0}^{J-1}\bm{c}_{lj}\Psi_{j}(\bm{\xi},\bm{\zeta}).$
(17)

Substituting this joint PCE for $\eta_{l}({\bm{\lambda}},\omega)$ into the KLE
yields an expression for the random field QoIs.

$\displaystyle f({\bm{\lambda}},\omega;{\bm{x}})$
$\displaystyle=f_{0}({\bm{x}})+\sum_{l=1}^{L}\bigg{(}\underbrace{\sum_{j=0}^{J-1}\bm{c}_{lj}\Psi_{j}({\bm{\xi}},{\bm{\zeta}})}_{\eta_{l}({\bm{\lambda}},\omega)}\bigg{)}\sqrt{\mu_{l}}\,\phi_{l}({\bm{x}})$
(18)

Simple rearrangement of terms allows us to identify the surrogate model Eq.
(18) as a PCE with coefficients that depend on the deterministic coordinates
$\bm{x}$. We refer to the result as a _KLPC surrogate_ , (i.e. a surrogate
constructed via KLE and PCEs)

$\displaystyle f({\bm{\lambda}},\omega;{\bm{x}})$
$\displaystyle=\sum_{j=0}^{J-1}\bigg{(}\underbrace{f_{0}({\bm{x}})\delta_{j0}+\sum_{l=1}^{L}\bm{c}_{lj}\sqrt{\mu_{l}}\,\phi_{l}({\bm{x}})}_{{\bm{c}}^{\text{\tiny
KLPC}}_{j}({\bm{x}})}\bigg{)}\Psi_{j}({\bm{\xi}},{\bm{\zeta}}).$ (19)

The Kronecker delta function $\delta_{j0}$ includes the process mean
$f_{0}({\bm{x}})$ in only the constant $(j=0)$ term of the KLPC.

Since the KLPC is indeed a PCE, we take advantage of the analytical
expressions for the moments and Sobol indices in terms of the PCE coefficients
to perform GSA without incurring sampling error. Starting from Eq. (3), the
first and second moments of the KLPC given in Eq. (19) are

$\displaystyle\mathbb{E}[\bm{y}({\bm{x}})]$ $\displaystyle=c^{\text{\tiny
KLPC}}_{0}({\bm{x}})=f_{0}({\bm{x}})+\sum_{l=1}^{L}\bm{c}_{l0}\sqrt{\mu_{l}}\,\phi_{l}({\bm{x}})$
(20) $\displaystyle\mathbb{V}[\bm{y}({\bm{x}})]$
$\displaystyle=\sum_{j=1}^{J-1}\big{(}c^{\text{\tiny
KLPC}}_{j}({\bm{x}})\big{)}^{2}||\Psi_{j}||^{2}=\sum_{j=1}^{J-1}\bigg{(}\sum_{l=1}^{L}\bm{c}_{lj}\sqrt{\mu_{l}}\,\phi_{l}({\bm{x}})\bigg{)}^{2}||\Psi_{j}||^{2}.$
(21)

Unlike the moment equations given in Eq. (3), the moments for the KLPC now
depend on a particular deterministic coordinate of the system ${\bm{x}}$ via
the KLPC coefficients $\bm{c}_{j}^{\text{\tiny KLPC}}({\bm{x}})$.

The main effect sensitivity indices are similarly found from Eq. (5) by
substituting in the KLPC coefficients.

$\displaystyle S_{i}$
$\displaystyle=\frac{1}{\mathbb{V}[\bm{y}({\bm{x}})]]}\sum_{j\in\mathcal{J}_{S_{i}}}\big{(}c^{\text{\tiny
KLPC}}_{j}({\bm{x}})\big{)}^{2}||\Psi_{j}||^{2}=\frac{1}{\mathbb{V}[\bm{y}({\bm{x}})]}\sum_{j\in\mathcal{J}_{S_{i}}}\bigg{(}\sum_{l=1}^{L}\bm{c}_{lj}\sqrt{\mu_{l}}\,\phi_{l}({\bm{x}})\bigg{)}^{2}||\Psi_{j}||^{2},$
(22) $\displaystyle\text{where
}\mathcal{J}_{S_{i}}:=\\{j:p_{i}(p)>0,p_{k}(p)=0\text{ for }k\neq i\\}.$

In the expressions for both the output variance and Sobol indices, the
squared-norm $||\Psi_{j}||^{2}$ is defined analogously to Eq. (4)

##  4 Results

To illustrate the methodology, we consider a heterogeneous gas-surface
chemical catalysis system, which we model using a mesoscopic lattice kinetic
Monte Carlo (KMC) construction [1, 31]. We presume an elementary-step chemical
kinetic mechanism employing a number of reactions, with uncertain rate
constants. At the same time, KMC employs intrinsic noise that models the
randomness of chemical reaction at the microscale.

Forward Processes | Reverse Processe
---|---
$[1]$ Adsorption: | CO $\xrightarrow[]{k_{1}}$ CO(`cus`) | Desorption: | CO(`cus`) $\xrightarrow[]{k_{-1}}$ CO
$[2]$ Adsorption: | CO $\xrightarrow[]{k_{2}}$ CO(`br`) | Desorption: | CO(`br`) $\xrightarrow[]{k_{-2}}$ CO
$[3]$ Adsorption: | O2 $\xrightarrow[]{k_{3}}$ O(`cus`) + O(`cus`) | Desorption: | O(`cus`) + O(`cus`) $\xrightarrow[]{k_{-3}}$ O2
$[4]$ Adsorption: | O2 $\xrightarrow[]{k_{4}}$ O(`br`) + O(`br`) | Desorption: | O(`br`) + O(`br`) $\xrightarrow[]{k_{-4}}$ O2
$[5]$ Adsorption: | O2 $\xrightarrow[]{k_{5}}$ O(`br`) + O(`cus`) | Desorption: | O(`br`) + O(`cus`) $\xrightarrow[]{k_{-5}}$ O2
$[6]$ Diffusion: | CO(`cus`) $\xrightarrow[]{k_{6}}$ CO(`cus`) |  |
$[7]$ Diffusion: | CO(`br`) $\xrightarrow[]{k_{7}}$ CO(`br`) |  |
$[8]$ Diffusion: | CO(`cus`) $\xrightarrow[]{k_{8}}$ CO(`br`) | Diffusion: | CO(`br`) $\xrightarrow[]{k_{-8}}$ CO(`cus`)
$[9]$ Diffusion: | O(`cus`) $\xrightarrow[]{k_{9}}$ O(`cus`) |  |
$[10]$ Diffusion: | O(`br`) $\xrightarrow[]{k_{10}}$ O(`br`) |  |
$[11]$ Diffusion: | O(`cus`) $\xrightarrow[]{k_{11}}$ O(`br`) | Diffusion: | O(`br`) $\xrightarrow[]{k_{-11}}$ O(`cus`)
$[12]$ Formation: | CO(`cus`) + O(`cus`) $\xrightarrow[]{k_{12}}$ CO2 |  |
$[13]$ Formation: | CO(`br`) + O(`br`) $\xrightarrow[]{k_{13}}$ CO2 |  |
$[14]$ Formation: | CO(`br`) + O(`cus`) $\xrightarrow[]{k_{14}}$ CO2 |  |
$[15]$ Formation: | CO(`cus`) + O(`br`) $\xrightarrow[]{k_{15}}$ CO2 |  |
Table 1: Summary of processes involved in the oxidation of CO to CO2 via a
RuO2 catalyst. Forward reactions are given in the left column, and any
associated reverse reactions listed in the right column. Adsorbed species are
indicated by the type of binding sites: br or cus. The rates of the forward
(reverse) reactions are denoted above the reaction arrows by $k_{i}$
($k_{-i}$). While diffusion is not typically written as a reaction in this
way, we include the diffusion events for completeness of the mechanism. Here
the reaction arrows associated with diffusion events indicate change of
position on the lattice and not change to chemical structures.

As an example problem, we consider CO oxidation on a RuO2(110) surface,
represented at atomic scales as a lattice with alternating coordinatively
unsaturated (`cus`) sites and bridge (`br`) binding sites where unimolecular
adsorption (desorption) of CO and dissociative adsorption (desorption) of O2
may occur. The adsorbed species can make diffusive hops to neighboring vacant
sites and nearest-neighbor O and CO can react, producing gaseous CO2. The
system is modeled using 22 permissible processes, where a ‘process’ on the
lattice is specified by event type, chemical species, and lattice site(s)
involved – for example, diffusion of CO from a `cus` site to a `br` site (see
[1, 8] for a more thorough discussion of the mechanism).

Associated with each process is a reaction rate constant which is treated as
an uncertain parameter. The CO2 production reactions, as well as the species
diffusion processes between identical lattice sites, are modeled as
irreversible ‘forward’ processes with associated forward rate constants
$k_{f}$, while diffusion between non-identical lattice sites and
adsorption/desorption are treated as reversible processes. The reversible
processes have forward and reverse rate constants $k_{f}$ and $k_{r}$ that are
thermodynamically balanced. The set of included reactions are summarized in
Table 1.

The uncertain input parameter vector $\bm{\lambda}$ consists of the logarithms
of the $15$ forward reaction rates, specified as uniform random variables
around given nominal values with a scale factor vector $\bm{r}$:

$\bm{\lambda}:=\\{\bm{k}\,|\,\ln\bm{k}\sim\text{Unif}(\ln\bm{k}_{\text{\tiny
nom}}-\ln\bm{r},\ln\bm{k}_{\text{\tiny nom}}+\ln\bm{r})\\}.$

The QoI considered here is the production rate of CO2 over time which we
define as $y(t)=n_{p}/(\alpha\Delta t)$ where $n_{p}/\Delta t$ is the number
of CO2 formations per time increment and $\alpha$ is the lattice area. That
is, the production rate is given by the number of CO2 formation events per
time per unit area. Clearly, $\alpha$ will scale with the number of lattice
cells $n_{c}$.

Figure 1 illustrates the time-evolution of the QoI for varying parameter
values. Clearly, the production rate of the chemical system changes
significantly with time. The effect of the parameters ${\bm{\lambda}}$ is
illustrated by the divergence of the sample solution medians to distinct
trajectories over time, while the magnitude of intrinsic noise $\omega$ is
indicated by the colored band surrounding each sample median. Looking ahead to
global sensitivity analysis, the fractional variance contributions of
intrinsic noise and parametric uncertainty are expected to be dominated by the
latter over time, since the color bands seem well separated, except possibly
where path crossing compresses the overall variance of the model.

Fig. 1: Production rate over time for CO2. Each color corresponds to a
different parameter configuration ${\bm{\lambda}}^{(n)}$ where $n$ is fixed.
Here the intrinsic stochasticity is highlighted by plotting individual
trajectories which correspond to different random seeds $\omega$. The original
trajectories are plotted on the left to illustrate the value separation within
and between parametric samples. However, we eschew the large value ranges by
working with logarithms of the QoI, as shown on the right.

###  4.1 Quantitative measures for results

We illustrate construction of the joint PCE for
$\eta_{l}({\bm{\lambda}},\omega)$ and the KLPC surrogate for the production
rate of CO2 over time as our main QoI. The uncertain parameter vector is a set
of 15 forward (log) reaction rate constants ${\bm{\lambda}}\,=\,(\ln
k_{1},...,\ln k_{15})$. Here we generate $N=1000$ parametric configurations
sampled around nominal rates proposed in [8] with an uncertainty scale factor
of $r=1.50$ for each parameter. For each fixed parameter configuration we use
our open-souce Kinetic Catalysis KMC software `KinCat` [6] to generate $M=500$
replica solutions of the stochastic model on a $50\times 50$ lattice. The
initial condition of the surface is generated by selecting 25% of the sites at
random and scattering CO molecules and O atoms on these sites in even
proportions.

In this work, the joint PCE construction for $\eta_{l}({\bm{\lambda}},\omega)$
is undertaken as a set of sequential steps. We first examine the quality of
the stochastic and parametric PCEs separately to demonstrate that each step of
the surrogate construction is sufficiently accurate. It follows that we must
consider multiple measures of error related to the PCE goodness-of-fit, as the
stochastic and parametric PCEs define fundamentally different approximations
of the model outputs and cannot be evaluated in the same way.

Whenever point estimates are appropriate, we compare surrogate predictions to
the training data via the relative root mean-squared error (rRMSE) defined as

$\displaystyle E_{\text{\tiny
rRMSE}}:=\sqrt{\frac{\sum_{n=1}^{N}(v^{\text{\tiny train}}_{n}-v^{\text{\tiny
pred}}_{n})^{2}}{\sum_{n=1}^{N}\big{(}v^{\text{\tiny train}}_{n}\big{)}^{2}}}$
(23)

where $v^{\text{\tiny train}}_{n}$ and $v^{\text{\tiny pred}}_{n}$ represent
the training and predicted values, respectively.

###  4.2 KLE for production rate over time

The first step is to find a suitable compact representation of the time-
evolving system. Random fields are well represented by the KLE, and highly
correlated fields can usually withstand substantial truncation of the
expansion without suffering significant loss of information. As noted in
Section 3, the rate of spectral decay dictates how early truncation can occur.
That is, the rate of variance accumulation in the KLE terms is tied to the
rate of decrease in the eigenvalue magnitudes, since the captured variance is
given by the sum of the eigenvalues.

Fig. 2: Spectral decay is shown by plotting the (log) eigenvalue against
expansion index $l$ (left). Accumulation of explained variance in the KLE is
demonstrated by plotting the explained variance given by Eq. (16) against
proposed truncation points $L$ (right). The crossing point of the gray lines
indicates where the KLE captures at least 99% of the explained model variance.

Figure 2 clarifies the relationship between the spectral decay and
accumulation of explained variance. Here we take advantage of the ordering of
the eigenvalues and plot the (log) eigenvalues against the expansion index
$l$. The early eigenvalues exhibit orders-of-magnitude value changes that
result in gaps in the spectrum and large changes in the accumulated variance.
As the rate of spectral decay slows, the accumulation of variance similarly
levels off. Typically the truncation point is chosen to preserve some
proportion of the output variance.

Fig. 3: PDFs of the model samples at the final simulation time (gray) and KLE
representation at different truncation points. The number of terms in the KLE
increases from left to right.

The effect of the truncation point $L$ on the error in the KLE of the training
samples is illustrated in Figures 3 and 4. PDFs of the _combined_ stochastic-
parametric training samples are plotted in gray while KLE reconstructions of
these data are plotted in pink ($L=1$), blue ($L=3$), and green ($L=39$). Note
that $L=1$ produces the most compact representation of the data possible,
while $L=3$ and $L=39$, according to Figure 2, capture at least 99% and 99.99%
of the total explained variance, respectively.

Fig. 4: PDFs of the model samples (gray) and KLE representation over time at
different truncation points: $L=1$ (pink), $L=3$ (blue), and $L=39$ (green).
Note that the vertical distances between PDFs is a log-scale of the simulation
time to improve visibility of PDFs. The values along the $x$-axis indicate the
mean of each PDF over time.

Figure 4 illustrates the quality of the KL representation of these models
relative to the training data over time, while Figure 3 compares the PDFs of
the QoI and its KL representation at the final KMC simulation time once the
trajectories have (approximately) reached equilibrium. Increasing the number
of KLE terms results in a more accurate KL approximation across all simulation
times, but comes at the expense of higher-dimensional representations of the
data. The key is to find a reasonable balance between the accuracy of the
model and the number of features we must construct surrogates for. From this
comparison, we select the KLE with $L=3$ terms. This captures slightly more
than 99% of the explained variance and produces a reasonably accurate
approximation of the training data, while curbing the number of elements that
must be estimated with joint PCEs.

Fig. 5: Production rate over time for CO2. Each color corresponds to a
different parameter configuration ${\bm{\lambda}}^{(n)}$, $n$ fixed. The
ensemble medians are plotted as thick lines and the $5$th/$95$th quantiles as
a band around the median. The original model samples (left) are compared to
their KLE representation when $L=3$ (right). We see that the KLE solution
trajectories reproduce the correct behavior over time for this level of
truncation.

We compare the selected KLE representation to the training samples in Figure 5
and see that the KLE solution trajectories reproduce the correct behavior over
time. The medians of the solution trajectories are plotted as thick, colored
lines for individual parametric samples ${\bm{\lambda}}^{(n)}$, $n$ fixed, and
the variation in samples indicated by the colored band surrounding each median
which represent the $5$th/$95$th percentiles.

###  4.3 Stochastic PCE

With a reasonable low-dimensional representation of the random field data in
hand, we now seek stochastic PCE representations of the random coefficients
$\eta_{l}({\bm{\lambda}},\omega)$ for each $l=1,2,3$. The PCEs take the form
given in Eq. (7).

We limited the stochastic PCEs to first-order expansions of the germ data, as
higher-order expansions appeared to offer little gain in accuracy for the QoI
while rapidly increasing the number of parametric surrogates needed to
construct the joint PCE. Since the sample replicas capture intrinsic model
noise and are (approximately) normally distributed, we use the Gauss-Hermite
polynomials as the expansion basis.

Due to lack of control over $\omega$ and $\bm{\zeta}$, the training samples
and the stochastic PCE predictions of $\bm{y}(t)$ are only guaranteed to agree
in distribution. Therefore, we assess the overall performance of the
stochastic PCEs by comparing statistics of the training model simulations and
the PCE-predicted random coefficients. For each fixed ${\bm{\lambda}}^{(n)}$,
we take sample replicas to find the sample means and standard deviations for
PCE predictions of $\bm{\eta}^{\text{\tiny
PCE}}({\bm{\lambda}}^{(n)},\omega)$. We then plot each component of the
predicted means (standard deviations) against the true sample means (standard
deviations) in parity plots.

Fig. 6: Parity plots of the mean (top row) and standard deviation (bottom row)
of the random KLE coefficients $\eta_{l}({\bm{\lambda}},\omega)$ predicted by
the stochastic PCE. The rRMSE for the mean and standard deviation of the set
of combined samples are $0.0043$ and $0.0272$, respectively.

Figure 6 demonstrates this comparison, indicating that the stochastic PCEs
predict the sample means with high accuracy for all three KLE coefficient
variables $\eta_{1}({\bm{\lambda}},\omega),\eta_{2}({\bm{\lambda}},\omega)$,
and $\eta_{3}({\bm{\lambda}},\omega)$. A similar comparison of the predicted
standard deviations reveals that the stochastic PCEs predict the variance of
the distributions well, on average, but the spread of the predictions
increases with the amount of noise in the random coefficient (see Figure 6).
We reduce the spread by increasing the number of samples in the training
dataset. Further, the bandwidth used to construct the Rosenblatt maps must be
chosen with care to avoid severely biasing the predicted standard deviations
of samples. Particularly, a bandwidth based on Silverman’s rule-of-thumb [43]
is not optimal for moment approximation and leads to systematic over-
estimation of predicted variances in the stochastic samples. Here the
bandwidth was chosen empirically, and automated optimization of these values
will be addressed in future work.

###  4.4 Parametric PCE

Construction of stochastic PCEs for the $\eta_{l}({\bm{\lambda}},\omega)$,

$\displaystyle\eta_{l}({\bm{\lambda}},\omega)=\sum_{s=0}^{3}\bm{z}_{ls}({\bm{\lambda}})\Psi_{s}(\bm{\zeta}),\quad
l=1,2,3$ (24)

yields model outputs $\bm{z}_{ls}({\bm{\lambda}})$ that directly correspond to
the parametric input samples ${\bm{\lambda}}$. Here the subscript $l$
associates explicitly a particular $\eta_{l}$ with the coefficients of its
representation, while subscript $s$ indicates the particular element of the
sum. We now construct parametric PCEs to represent this input-output mapping.

Since the elements of $\bm{\lambda}$ are assumed to be distributed uniformly
on $[\ln\bm{k}_{\text{\tiny nom}}-\ln\bm{r},\ln\bm{k}_{\text{\tiny
nom}}+\ln\bm{r}]$, we use Legendre polynomials as the expansion basis and
obtain the simple rescaling relationship as a first-order Legendre-Uniform
PCE,

$\displaystyle\lambda_{i}=\ln k_{\text{\tiny nom,}\,i}+\xi_{i}\ln r_{i},\text{
for }i=1,\dots,\tilde{d}.$ (25)

As described in Section 2, we then construct a parametric PCE of the form Eq.
(9),

$\displaystyle\bm{z}_{ls}=\sum_{p=0}^{P_{ls}-1}\bm{a}_{lsp}\Psi_{p}({\bm{\xi}}).$
(26)

This gives the stochastic PCE coefficients $\bm{z}_{ls}$ as functions of the
rescaled input, or germ, $\bm{\xi}$.

The parametric PCEs are chosen to have at most second-order expansions of the
coefficient data $\bm{z}_{ls}$, and the optimal orders for each $s$ are
selected via model evidence which is analytically available for this
regression step [40]. We reduce the number of expansion terms in the
parametric PCE by estimating the coefficients with Bayesian Compressive
Sensing (BCS), which adaptively finds a sparse, optimal basis set for the PCE
[2, 22, 42]. We examined the accuracy of the parametric PCEs using randomized
75%:25% training/testing splits of the data with the rRMSE metric on
coefficient approximation, followed by a final training of the parametric
surrogate using all samples.

Unlike the stochastic PCEs, the parametric PCE preserves pointwise mappings
between the PCE germ $\bm{\xi}$, parametric inputs $\bm{\lambda}$, and
$\bm{z}_{ls}(\bm{\lambda})$ coefficient values, hence we do not need to resort
to moment comparisons and can examine pointwise parity plots of the true and
PCE-predicted outputs.

Fig. 7: Parity plots of parametric PCE predicted coefficients
$\bm{z}_{ls}^{\text{\tiny PCE}}(\bm{\lambda})$ constructed in Eq. (9) vs true
coefficients $\bm{z}_{ls}(\bm{\lambda})$ found in Eq. (8) for all ensemble
members. The coefficients $\bm{z}_{l}(\bm{\lambda})$ are grouped by random KLE
element $\eta_{l}({\bm{\lambda}},\omega)$.

Figure 7 illustrates the approximation quality of each parametric PCE
coefficient by plotting the predicted values $\bm{z}_{ls}^{\text{\tiny
PCE}}({\bm{\lambda}})$ against the true $\bm{z}_{ls}({\bm{\lambda}})$ for each
$s=0,1,2,3$. We reduce the total number of plots needed to show this result by
collecting all $\bm{z}_{ls}$ coefficients of a particular $\eta_{l}$ together
in one plot.

The close fit of the models to the training data for each coefficient is
supported by the computed rRMSE for each collection of coefficients:
rRMSE$(\bm{z}_{1})=0.019$, rRMSE$(\bm{z}_{2})=0.028$,
rRMSE$(\bm{z}_{3})=0.048$.

###  4.5 Joint PCE, KLPC, and Global Sensitivity Analysis

Finally, we assemble the joint PCE predictions sequentially. We substitute the
parametric PCE results for $\bm{z}_{ls}$ into the stochastic PCE to obtain the
joint PCE predictions of $\eta_{l}({\bm{\lambda}},\omega)$:

$\displaystyle\eta_{l}({\bm{\lambda}},\omega)=\sum_{s=0}^{3}\bigg{(}\underbrace{\sum_{p=0}^{P_{ls}-1}\bm{a}_{lsp}\Psi_{p}({\bm{\xi}})}_{\bm{z}_{ls}}\bigg{)}\Psi_{s}(\bm{\zeta})$
(27)

These are further folded into Eq. (18), the KLE, to obtain physical QoI
samples of the form of Eq. (19) from the surrogate.

As described in Section 3, the KLPC surrogate is itself a PCE and its quality
is dependent upon both the quality of the samples of
$\eta_{l}({\bm{\lambda}},\omega)$ generated by the joint PCE and the number of
KLE modes included in the initial low-dimensional representation of the random
fields. In Figures 8-9 we compare these surrogate predictions to the training
data on the set of combined stochastic-parametric samples.

Fig. 8: Violin plots demonstrating the accuracy of the joint PCE predictions
in spectral space.

Examining the quality of both the joint PCE predictions of
$\eta_{l}({\bm{\lambda}},\omega)$ (see Figure 8) and the KLPC predictions of
$f({\bm{\lambda}},\omega;{\bm{t}})$ (see Figure 9) show that the surrogate is
accurate on both the spectral space of KLE elements and physical space of
QoIs.

Fig. 9: Ridge plots compare the combined KLPC predictions over time to the
training samples. The values along the $x$-axis indicate the mean of each PDF
over time.

As a bi-product of the KLPC construction (19), one can extract global
sensitivity indices and attribute the variance in the output QoI to either the
uncertainties in the parameters or intrinsic noise from the stochastic model.
As illustrated in Figure 10, we find the adsorption rates to be the largest
contributors to the output variance for all times in the simulation, although
the importance of the formation rates increases as the lattice fills and the
number of formation events stabilizes. The Sobol index from the intrinsic
model noise is a group sensitivity that is found by summing the effects of all
elements of $\bm{\zeta}$.

The parametric uncertainties overall seem to dominate the intrinsic noise as
one would expect intuitively based on Figure 1. Effectively, the high
parametric uncertainty washes out the influence of intrinsic noise in the
system. Meanwhile, the highest noise contribution to the output uncertainty
occurs where the number of CO2 formation events rapidly increases. Here, path
crossing compresses the overall variance of the model output due to the
parameters and increases the sensitivity of the QoI to the influence of
intrinsic noise.

Finally, for this GSA result, we point out the effect of cross-interactions
introduced by the second order expansion terms. It is known that main
sensitivities will sum to one if cross terms are absent from the polynomial
expansions. When cross terms are present, they siphon off a fraction of the
main sensitivities and the sum of the main indices is less than one. The
results in Figure 10 indicate that the contribution of cross terms is small
initially, growing slightly larger in time.

Fig. 10: Sobol indices over time for uncertain rate constants $k_{i}$ of
processes described in Table 1. Only adsorption rates (blues), certain
formation rates (greens), and noise (red) contribute appreciably to the
sensitivity analysis. Sensitivity indices of these most important parameters
are labeled by associated process.

One strength of this approach for stochastic surrogate construction is the
following: the effect of intrinsic noise on the model output can be
systematically quantified and thereby controlled. That is, we can use the
joint PCE (KLPC) construction as a model selection tool that guides, for
example, the ensemble size chosen to reduce intrinsic noise when estimating
integrated QoIs or the size of the lattice used in a KMC simulation to
generate data samples. This is demonstrated in Figure 11, where the Sobol
indices over time are plotted for KLPC surrogates constructed from data
generated from KMC simulations which differ only in the lattice size $n_{c}$.
As the lattice size increases, the presence of intrinsic noise (plotted with
red) is diminished until it is negligible. From this we conclude that the
largest lattice with dimension $n_{c}=100$ is optimal, in that it is the
smallest lattice which allows nearly noise-free samples to be generated for
the catalytic system under investigation.

Fig. 11: Sobol indices over time for uncertain adsorption rates (blues),
formation rates (greens), and noise (red) of processes described in Table 1.
The lattice sizes increases from left-to-right and have dimensions $n_{c}=25$
(left), $n_{c}=50$ (center), and $n_{c}$=100 (right). Sensitivity indices of
the most important parameters are labeled by associated process.

In Figure 11, we ensure consistency amongst the training datasets generated
with KMC by using the same set of parametric inputs ${\bm{\lambda}}$ when
constructing trajectories. We further initiate each simulation from the same
set of initial random seed values. Thus, the only significant difference in
the results comes through the influence of the lattice size, $n_{c}$. To
reduce the computational effort of generating training data for this
comparison, we reduced the parametric uncertainty and used a smaller number of
parametric inputs samples than in the previous results.

###  4.6 Using the KLPC to obtain sensitivity indices for log-normal rates

Constructing the surrogate model from independent, log-uniform rate samples
ensures that the joint PCE surrogate is uniformly accurate over the hypercube
$[a_{i},b_{i}],i=1,\dots,15$. Here the bounds are given from the nominal rates
and uncertainty scale factors via $\ln k_{\text{\tiny nom},i}\pm\ln r_{i}$.
However, chemical reaction rates are best described with log-normal
distributions. To obtain more physically meaningful GSA results, we make use
of the following procedure, presuming independent lognormal rate
uncertainties.

  * 1.

With the surrogate model Eq. (18), we generate $N^{\prime}$ samples from the
standard normal germ $\bm{\nu}\sim\mathcal{N}(0,1)$. These can be mapped to
physical rate values via $\ln k_{i}=\mu_{i}+\nu_{i}\sigma_{i}$ – so that $\ln
k_{i}\sim\mathcal{N}(\mu_{i},\sigma_{i}^{2})$ – where $\mu_{i}$ and
$\sigma_{i}$ are the mean and standard deviation of the $i$th process rate,
respectively. Denote the collection of log-normally distributed rates
${\bm{\theta}}=[\ln k_{1},\dots,\ln k_{15}]$.

    * (a)

The mean is naturally given by the nominal rate values $\mu_{i}=\ln
k_{\text{\tiny nom},i}$.

    * (b)

The sample standard deviation is obtained from confidence intervals. That is,
given a sample set of size $N^{\prime}$ with mean $\ln k_{\text{\tiny nom},i}$
and standard deviation $\sigma_{i}$, there is some probability $p(z)$ that a
value $\ln k_{i}$ lies within $z$ standard deviations of the mean:

$\displaystyle\ln k_{i}=\ln k_{\text{\tiny
nom},i}\pm\frac{z\sigma_{i}}{\sqrt{N^{\prime}}}.$

If we wish to ensure that $\ln k_{i}\in[\ln k_{\text{\tiny nom},i}-\ln
r_{i},\ln k_{\text{\tiny nom},i}+\ln r_{i}]$ with probability $p(z)$, then set
$\ln r_{i}=\frac{z\sigma_{i}}{\sqrt{N^{\prime}}}$. Solving for $\sigma_{i}$
gives the sample standard deviation we should choose:

$\displaystyle\sigma_{i}=\sqrt{N^{\prime}}\,\frac{\ln r_{i}}{z}.$ (28)

Note that any sample $\ln k_{i}$ may fall outside of the desired interval with
probability $1-p(z)$. Standard normal tables can provide the one-to-one
correspondence between $p$ and $z$, and one can chose $z$ according to a given
threshold value of $p$. With $p$ being close to one, we still obtain a
reasonable estimate for the Sobol indices for GSA with only small number of
normal samples falling outside the uniform range. While the GSA is accurate
enough, the model predictions of the QoI at these inputs is untrustworthy as
they are outside the range of validity of the pre-constructed surrogate.

  * 2.

With the surrogate model Eq. (18), germ samples $\bm{\nu}$, and associated
physical rate values ${\bm{\theta}}$ in hand, we evaluate the surrogate model
at values of ${\bm{\theta}}$ to obtain associated samples of the QoI.

$\displaystyle g({\bm{\theta}},\omega;{\bm{t}})=f^{\text{\tiny
PCE}}({\bm{\xi}}({\bm{\theta}}),{\bm{\zeta}};{\bm{t}})$ (29)

As before, we transform independent elements of the input rates $\theta_{i}$
to the compact interval of support $[-1,1]$ via rescaling $\theta_{i}=\ln
k_{\text{\tiny nom},i}+\xi_{i}\ln r_{i}$. These transformed values
${\bm{\xi}}({\bm{\theta}})$ are evaluated in the surrogate, Eq. (29).

  * 3.

Now, given the new dataset of log-normal rates and associated QoIs,
$\\{{\bm{\theta}},g({\bm{\theta}},\omega;{\bm{t}})\\}$, we can construct a
_second_ joint PCE surrogate such that _both_ the stochastic and parametric
PCEs utilize a _standard normal_ germ.

$\displaystyle g({\bm{\theta}},\omega;{\bm{t}})\approx g^{\text{\tiny
PCE}}(\bm{\nu},{\bm{\zeta}};{\bm{t}})=\sum_{j=0}^{J-1}\bigg{(}\underbrace{f_{0}({\bm{t}})\delta_{0j}+\sum_{l=1}^{L}\tilde{\bm{b}}_{lj}\sqrt{\mu_{l}}\phi_{i}({\bm{t}})}_{\bm{c}^{\text{\tiny
KLPC}}({\bm{t}})}\bigg{)}\Psi_{j}(\bm{\nu},{\bm{\zeta}})$ (30)

The new surrogate Eq. (30) is used to propagate uncertainty from the physical
rates ${\bm{\theta}}$ and Sobol indices are obtained from the PCE coefficients
$\bm{c}^{\text{\tiny KLPC}}({\bm{t}})$.

$\displaystyle S_{i}$
$\displaystyle=\frac{1}{\mathbb{V}[g({\bm{t}})]}\sum_{j\in\mathcal{J}_{S_{i}}}c_{j}(t)^{2}||\Psi_{j}||^{2}=\frac{1}{\mathbb{V}[g({\bm{t}})]}\sum_{j\in\mathcal{J}_{S_{i}}}\small{\bigg{(}\sum_{l\leq
L}\tilde{\bm{b}}_{lj}\sqrt{\mu_{l}}\,\phi_{l}({\bm{t}})\bigg{)}^{2}}||\Psi_{j}||^{2}$
(31)

As before, $||\Psi_{j}||^{2}$ is defined analogously to Eq. (4).

In summary, given nominal rates and uncertainty scale factors, we can generate
log-normal reaction rate samples ${\bm{\theta}}$, and then evaluate these
rates in the original Legendre-Uniform KLPC surrogate to quickly obtain
samples of the QoI. As noted previously, the generative nature of the KLPC
allows us to generate an arbitrary number of sample replicas for each fixed
${\bm{\theta}}^{(n)}$.

We use these surrogate-generated samples of the input/output data
$\\{{\bm{\theta}},f^{\text{\tiny
PCE}}({\bm{\xi}}({\bm{\theta}}),\omega;{\bm{t}})\\}$ to construct a Gauss-
Hermite KLPC Eq. (30) which exclusively uses standard normal germs throughout
the construction. The coefficients $\tilde{\bm{b}}$ of this normal-based
surrogate are then used to compute Sobol indices over time. We plot these
indices in Figure 12 and compare to the Sobol indices obtained from the
original Legendre-Uniform surrogate and to the Sobol indices obtained from
directly constructing a Gauss-Hermite PCE-based surrogate on KMC samples
generated with the same log-normal rates.

Fig. 12: Comparison of Sobol indices over time for KLPC surrogates constructed
from three different training datasets. Left: Training data are generated with
KMC and the input parameters ${\bm{\lambda}}$ follow log-uniform
distributions. Right: Training data are generated with KMC and the input
parameters $\bm{\theta}$ follow log-normal distributions. Center: Training
data are generated for the same log-normal input parameters $\bm{\theta}$ by
sampling the uniformly-accurate KLPC. Sensitivity indices of the most
important parameters are labeled by associated process as given in Table 1.

The GSA result obtained from the Gauss-Hermite KLPC that is trained on data
sampled from the original Legendre-Uniform surrogate does not see an
appreciable amplification of noise. Regardless of the distribution of the
input parameters, the identity of and significance of the sensitive parameters
are preserved over time.

##  5 Conclusion

We developed a joint PCE model for representing both parametric uncertainty
and intrinsic noise of a stochastic model in a compact surrogate form. The
construction combines the conventional parametric PCE with a stochastic PCE
component constructed with the help of the inverse Rosenblatt transformation.
The joint PCE surrogate is generative as the underlying stochastic germ is
represented in a way that can be resampled to generate new realizations of the
underlying stochastic model. We extended the joint PCE construction to spatio-
temporal random field realizations and showed that the resulting KLPC
surrogate still fits within the PCE framework while preserving the generative
qualities of the joint PCE.

The joint PCE (KLPC) construction enables analytical extraction of global
sensitivity indices via variance decomposition. This allows separating the two
sources of uncertainty: intrinsic noise and parametric uncertainty.
Decomposing the output variance in this way allows us to assess the impact
that intrinsic noise has on the model predictions and can be used to inform
the size and quality of the data set used to construct the surrogate. The
parametric PCE component of the algorithm, while relatively routine, also
includes model selection and basis selection strategies to address the
factorial growth of the basis set with the model order and input dimension.

We demonstrated the results on an example heterogeneous catalysis model.
Improving elements of the surrogate construction workflow, such as removing
the parameter independence assumption, investigating alternative bandwidth
selection strategies in the Rosenblatt map construction, and implementing a
rigorous cross-validation study to complement BCS are subjects of ongoing
work.

## Acknowledgments

We would like to thank our colleagues Judit Zádor and Kyungjoo Kim for their
contributions to this work.
Funding: This work was done within the Exascale Catalytic Chemistry (ECC)
Project, which is supported by the U.S. Department of Energy, Office of
Science, Basic Energy Sciences, Chemical Sciences, Geosciences and Biosciences
Division, as part of the Computational Chemistry Sciences Program.
This article has been authored by employees of National Technology &
Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the
U.S. Department of Energy (DOE). The employees co-own right, title and
interest in and to the article and are responsible for its contents. The
United States Government retains and the publisher, by accepting the article
for publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this article or allow others to do so, for United States
Government purposes. The DOE will provide public access to these results of
federally sponsored research in accordance with the DOE Public Access Plan
https://www.energy.gov/downloads/doe-public-access-plan.

## References

  * [1] M. Andersen, C. Panosetti, and K. Reuter, A practical guide to surface kinetic monte carlo simulations, Front. Chem., 7 (2019), pp. 1–24, doi: 10.3389/fchem.2019.00202.
  * [2] S. Babacan, R. Molina, and A. Katsaggelos, Bayesian compressive sensing using Laplace priors, IEEE T Image Process, 19 (2010), pp. 53–63.
  * [3] E. Baker, P. Barbillon, A. Fadikar, R. B. Gramacy, R. Herbei, D. Higdon, J. Huang, L. R. Johnson, P. Ma, A. Mondal, B. Pires, J. Sacks, and V. Sokolov, Analyzing stochastic computer models: A review with opportunities, Statist. Sci., 37 (2022), pp. 64–89.
  * [4] M. Binois, R. B. Gramacy, and M. Ludkovski, Practical heteroscedastic gaussian process modeling for large simulation experiments, J. Comput. Graph. Statist., 27 (2018), pp. 808–821.
  * [5] T. Crestaux, O. Le Maître, and J. Martinez, Polynomial chaos expansion for sensitivity analysis, Reliab. Eng. Syst. Saf., 94 (2009), pp. 1161–1172.
  * [6] C. Daniels, H. N. Najm, K. Kim, N. V. Roberts, and J. Zador, Kinetic catalysis (KinCat) software for kinetic Monte Carlo simulations of surface kinetics on heterogeneous computing architectures.  https://github.com/sandialabs/KinCat, 2023.
  * [7] B. Debusschere, H. Najm, P. Pébay, O. Knio, R. Ghanem, and O. Le Maître, Numerical challenges in the use of polynomial chaos representations for stochastic processes, SIAM J. Sci. Comput., 26 (2004), pp. 698–719.
  * [8] S. Döpking and S. Matera, Error propagation in first-principles kinetic monte carlo simulation, Chem. Phys. Lett., 674 (2017), pp. 28–32.
  * [9] R. Ghanem, D. Higdon, and H. Owhadi, eds., Handbook of Uncertainty Quantification, Springer, 2017.
  * [10] R. Ghanem and J. RedHorse, Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach., Phys. D, 133 (1999), pp. 137 – 144.
  * [11] R. Ghanem and P. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer Verlag, New York, 1991.
  * [12] H. Gutmann, A radial basis function method for global optimization, J. Glob. Optim., 19 (2001), pp. 201–227.
  * [13] P. Hall, J. Racine, and Q. Li, Cross-validation and the estimation of conditional probability densities, J. Amer. Statist. Assoc., 99 (2004), pp. 1015–1026.
  * [14] J. L. Hart, A. Alexanderian, and P. A. Gremaud, Efficient computation of Sobol’ indices for stochastic models, SIAM J. Sci. Comput., 39 (2017), pp. A1514–A1530.
  * [15] S. P. Huang, S. T. Quek, and K. K. Phoon, Convergence study of the truncated Karhunen-Lòeve expansion for simulation of stochastic processes, Internat. J. Numer. Methods Engrg., 52 (2001), pp. 1029 – 1043.
  * [16] B. Iooss and M. Ribatet, Global sensitivity analysis of computer models with functional inputs, Reliab. Eng. Syst. Saf., 94 (2009), pp. 1194–1204.
  * [17] M. N. Jimenez, O. P. Le Maître, and O. M. Knio, Nonintrusive polynomial chaos expansions for sensitivity analysis in stochastic differential equations, SIAM-ASA J Uncertain, 5 (2017), pp. 378–402.
  * [18] K. Karhunen, Zur spektraltheorie stochastischer prozesse, Ann. Acad. Sci. Fennicae, 37 (1946).
  * [19] R. Koenker and G. B. Jr., Regression quantiles, Econometrica, 46 (1978), pp. 33–50.
  * [20] O. Le Maître and O. Knio, Spectral Methods for Uncertainty Quantification, Springer, New York, NY, 2010.
  * [21] O. Le Maître, O. Knio, H. Najm, and R. Ghanem, A stochastic projection method for fluid flow I. Basic formulation, J. Comput. Phys., 173 (2001), pp. 481–511.
  * [22] N. Lüthen, S. Marelli, and B. Sudret, Sparse polynomial chaos expansions: Literature survey and benchmark, SIAM-ASA J Uncertain, 9 (2021), pp. 593–649.
  * [23] N. Lüthen, S. Marelli, and B. Sudret, A spectral surrogate model for stochastic simulators computed from trajectory samples, Comput. Methods. Appl. Mech. Engrg., 406 (2023), p. 115875.
  * [24] O. L. Maître and O. Knio, Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics (Scientific Computation), Springer, 1st edition. ed., April 2010.
  * [25] A. Marrel, B. Iooss, S. D. Veiga, and M. Ribatet, Global sensitivity analysis of stochastic computer models with joint metamodels, Stat. Comput., 22 (2012), pp. 833–847.
  * [26] Y. Marzouk, T. Moselhy, M. Parno, and A. Spantini, Sampling via measure transport: An introduction, in Handbook of Uncertainty Quantification, . H. O. R. Ghanem, D. Higdon, ed., Springer, 2017.
  * [27] Y. Marzouk, H. Najm, and L. Rahn, Stochastic spectral methods for efficient Bayesian solution of inverse problems, J. Comput. Phys., 224 (2007), pp. 560–586.
  * [28] Moutoussamy, Vincent, Nanty, Simon, and Pauwels, Benoît, Emulators for stochastic simulation codes, ESAIM: Proc., 48 (2015), pp. 116–155.
  * [29] J. N. Mueller, K. Sargsyan, and H. N. Najm, Polynomial chaos surrogate construction for stochastic models with parametric uncertainty, in 14th International Conference on Applications of Statistics and Probability in Civil Engineering (ICASP14), Dublin, Ireland, 2023.
  * [30] H. Najm, Uncertainty Quantification and Polynomial Chaos Techniques in Computational Fluid Dynamics, Annu. Rev. Fluid Mech., 41 (2009), pp. 35–52.
  * [31] M. Pineda and M. Stamatakis, Kinetic monte carlo simulations for heterogeneous catalysis: fundamentals, current status, and challenges, J. Chem. Phys., 156 (2022).
  * [32] M. Plumlee and R. Tuo, Building accelerate emulators for stochastic simulations via quantile kriging, Technometrics, 56 (2014), pp. 466–473.
  * [33] C. E. Rasmussen and J. Quinonero-Candela, Healing the relevance vector machine through augmentation, in Proceedings of the 22nd International Conference on Machine learning, ACM, 2005, pp. 689–696.
  * [34] D. J. Rezende and S. Mohamed, Variational inference with normalizing flows, 2015.
  * [35] M. Rosenblatt, Remarks on a multivariate transformation, Ann. Math. Stat., 23 (1952), pp. 470 – 472.
  * [36] J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn, Design and analysis of computer experiments (C/R: p423–435), Statist. Sci., 4 (1989), pp. 409–423.
  * [37] A. Saltelli, M. Ratto, S. Tarantola, F. Campolongo, E. Commission, and J. R. C. Ispra, Sensitivity analysis practices: Strategies for model-based inference, Reliab. Eng. Syst. Saf., 91 (2006), pp. 1109–1125.
  * [38] A. Saltelli, S. Tarantola, F. Campolongo, and M. Ratto, Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models, John Wiley & Sons, 2004.
  * [39] A. Saltellia, K. Aleksankina, W. Becker, P. Fennell, F. Ferretti, N. Holst, S. Li, and Q. Wu, Why so many published sensitivity analyses are false: A systematic review of sensitivity analysis practices, Environ. Model. Softw., 114 (2019), pp. 29–39.
  * [40] K. Sargsyan, Surrogate models for uncertainty propagation and sensitivity analysis, in Handbook of Uncertainty Quantification, . H. O. R. Ghanem, D. Higdon, ed., Springer, 2017, pp. 673–698.
  * [41] K. Sargsyan, B. Debusschere, H. Najm, and O. Le Maître, Spectral representation and reduced order modeling of the dynamics of stochastic reaction networks via adaptive data partitioning, SIAM J Sci. Comput., 31 (2010), pp. 4395–4421.
  * [42] K. Sargsyan, C. Safta, H. N. Najm, B. J. Debusschere, D. Ricciuto, and P. Thornton, Dimensionality reduction for complex models via Bayesian compressive sensing, Int. J Uncert. Quant., 4 (2014), pp. 63–93.
  * [43] B. Silverman, Density Estimation for Statistics and Data Analysis, Chapman and Hall, London, 1986.
  * [44] A. Smola, , and B. Schölkopf, A tutorial on support vector regression, Stat. Comp., 14 (2004), pp. 199–222.
  * [45] B. Sudret, Global sensitivity analysis using polynomial chaos expansions, Reliab. Eng. Syst. Saf., 93 (2008), pp. 964–979.
  * [46] L. Torossian, V. Picheny, R. Faivre, and A. Garivier, A review on quantile regression for stochastic computer experiments, Reliab. Eng. Syst. Saf., 201 (2020), p. 106858.
  * [47] V. Vapnik, The nature of statistical learning theory, (1998).
  * [48] D. Xiu and G. Karniadakis, Modeling uncertainty in flow simulations via generalized polynomial chaos, J. Comput. Phys., 187 (2003), pp. 137–167.
  * [49] X. Zhu and B. Sudret, Replication-based emulation of the response distribution of stochastic simulators using generalized lambda distributions, Int. J Uncert. Quant., 10 (2020), pp. 249–275.
  * [50] X. Zhu and B. Sudret, Global sensitivity analysis for stochastic simulators based on generalized lambda surrogate models, Reliab. Eng. Syst. Saf., 214 (2021), p. 107815.
  * [51] X. Zhu and B. Sudret, Stochastic polynomial chaos expansions to emulate stochastic simulators, Int. J Uncert. Quant., (2023).