# On Active Learning for Gaussian Process-based Global Sensitivity Analysis

Mohit Chauhan, Mariel Ojeda-Tuz, Ryan Catarelli, Kurtis Gurley,
Dimitrios Tsapetis, Michael D. Shields

## Abstract

This paper explores the application of active learning strategies to
adaptively learn Sobol indices for global sensitivity analysis. We demonstrate
that active learning for Sobol indices poses unique challenges due to the
definition of the Sobol index as a ratio of variances estimated from Gaussian
process surrogates. Consequently, learning strategies must either focus on
convergence in the numerator or the denominator of this ratio. However, rapid
convergence in either one does not guarantee convergence in the Sobol index.
We propose a novel strategy for active learning that focuses on resolving the
main effects of the Gaussian process (associated with the numerator of the
Sobol index) and compare this with existing strategies based on convergence in
the total variance (the denominator of the Sobol index). The new strategy,
implemented through a new learning function termed the MUSIC (minimize
uncertainty in Sobol index convergence), generally converges in Sobol index
error more rapidly than the existing strategies based on the Expected
Improvement for Global Fit (EIGF) and the Variance Improvement for Global Fit
(VIGF). Both strategies are compared with simple sequential random sampling
and the MUSIC learning function generally converges most rapidly for low-
dimensional problems. However, for high-dimensional problems, the performance
is comparable to random sampling. The new learning strategy is demonstrated
for a practical case of adaptive experimental design for large-scale Boundary
Layer Wind Tunnel experiments.

Keywords: Sobol Index, Active Learning, Global Sensitivity Analysis, Gaussian
Process Regression, Kriging

##  1 Introduction

Active learning is widely used in computational physics-based modeling, and
increasingly in experimental studies [1], to serve a variety of objectives. In
these schemes, a machine learning (ML) model is employed in an iterative
routine to extract information from previous simulations/experiments and
identify, through some learning function the most informative
simulation/experiment to run next. This has been especially beneficial in the
fields of optimization and reliability, although it has also found wider
applications in uncertainty quantification (UQ) more generally.

In optimization, active learning is at the heart of Bayesian Optimization
(BO). In BO, Bayes Rule is used to define a learning function (termed an
acquisition function in BO) that seeks to identify the minimum (or maximum) of
some function. Comprehensive reviews of BO can be found in the following
references [2, 3, 4]. A detailed discussion of BO is beyond the scope of this
work, but the most common approach in these methods is to approximate the
function with a Gaussian Process (GP) and construct a learning function that
exploits the jointly Gaussian nature of the function to adaptively learn where
the sought extremum has the highest probability of lying. Among these
developments, the most widely used learning functions are the Expected
Improvement Function (EIF) [5, 6, 7] the Knowlege Gradient [8, 9, 10] and the
Entropy Search [11, 12] algorithms.

For reliability analysis, the objective is to estimate the probability of
failure of a physical system, which involves solving a high-dimensional
probability integral (expectation) over a binary classifier that indicates
either safety or failure. This integral is often solved statistically using
Monte Carlo methods wherein the integral is estimated as the statistical
expectation of the failure indicator by sampling from the input distribution
and assessing the indicator at each sample point. Active learning is then used
to approximate this indicator function efficiently to minimize the
computational burden associated with assessing system failure at each sample
with an expensive computational model. Numerous learning functions have been
developed for this task, with each attempting to efficiently approximate the
limit surface (surface that separates safe and failure domains) in different
ways, but in most cases using a GP. The first such method, termed the Expected
Feasibility Function (EFF) was proposed by Bichon et. al in 2008 [13]. This
was followed by the widely used U-function developed by Echard et al. [14] in
the Adaptive Kriging with Monte Carlo Simulation (AK-MCS) method. Since then,
various modifications have been made and these concepts have been integrated
with other reliability methods (such as subset simulation [15, 16, 17]) to
improve the efficiency, convergence, and accuracy of these methods [18, 19,
20, 21].

Although optimization and reliability are the most common application domains
for active learning, it has been applied more broadly in UQ, for example to
efficiently construct surrogate models (aka metamodels or emulators) for
expensive physics-based computational models. Lam [22], for example, proposed
the Expected Improvement for Global Fit (EIGF) learning function to improve
the convergence of GP surrogate models for UQ purposes. Active learning has
also been applied using other ML surrogate model forms, such as Polynomial
Chaos Expansions (PCE) [23, 24, 25] and Deep Neural Networks (DNNs) [26, 27,
28]. Furthermore, recent efforts have integrated active learning concepts into
the design of physical experiments for wind tunnel experiments[29, 30, 31, 1],
materials discovery [32, 33, 34], and biological sciences [35, 36, 37, 38].

Regardless of their application or objective, a common challenge in active
learning is to balance exploration and exploitation. That is, any active
learning method (i.e. learning function) must simultaneously use (exploit) the
information obtained from prior simulations/experiments to guide further
analyses and, at the same time, explore new regions of the parameter space
where little or no data has been collected and uncertainty is correspondingly
high. Balancing these objectives is a primary challenge in active learning
[39], especially in high-dimensional parameter spaces where exploration can be
particularly difficult [40].

In this work, we explore the potential for active learning for Global
Sensitivity Analysis (GSA), which to the authors’ knowledge, has not
previously been studied. Active learning for GSA is of particular interest
here because of its potential to aid in other UQ studies. In particular, GSA
can be used as an initial filter in UQ studies to reduce the number of random
variables that need to be considered in the study. For example, in Section 6
we consider a practical boundary layer wind tunnel (BLWT) experimental design
problem that has 10 input random variables and, through active learning,
determine that only 3 random variables contribute appreciably to the variance
of the response. We can then practically ignore uncertainty in the other 7
random variables, thus making UQ more computationally tractable in this
complex physical setting. Hence, active learning for GSA allows us to
adaptively reduce the dimension rapidly on-the-fly and therefore speed up the
overall UQ task.

We show that active learning for GSA presents some unique challenges that do
not arise in other active learning domains. These challenges arise from the
definition of Global Sensitivity Indices (specifically Sobol indices [41]) as
the ratio of two variances, which makes the construction of a learning
function that directly targets the Sobol index difficult to derive. Instead,
we explore two different learning functions – one that targets efficient
learning for the numerator and one that targets the denominator. We show that
each can be beneficial under certain circumstances, but neither approach
universally outperforms random sampling. We explore the reasons for this, and
specifically demonstrate that this results from the challenges of minimizing
the error in a ratio. We discuss when the different approaches appear to be
beneficial using a series of analytical functions of varying dimensions. We
then demonstrate our use of active learning for GSA on the BLWT application of
interest.

##  2 Preliminaries

Here, we briefly review the foundations of Global Sensitivity Analysis using
Sobol Indices and Gaussian Process (GP) regression, upon which all further
analyses will depend.

###  2.1 Sobol Sensitivity Indices

Variance-based methods for GSA aim to decompose the variance of the model
output $y(\bm{X})$ into distinct contributions from each of the individual
random variables $X_{i}$ (main effects) and the interaction of multiple random
variables. In the Sobol method [41, 42, 43], this is done by first
constructing the high-dimensional model representation (HDMR) of the model as:

$y(\bm{X})=f_{0}+\sum_{i=1}^{d}f_{i}(X_{i})+\sum_{i=1}^{d}\sum_{j>i}f_{ij}(X_{i},X_{j})+\dots$
(1)

where

$\displaystyle f_{0}$ $\displaystyle=E[Y]$ (2) $\displaystyle f_{i}$
$\displaystyle=E[Y|X_{i}]-E[Y]$ $\displaystyle f_{ij}$
$\displaystyle=E[Y|X_{i},X_{j}]-f_{i}-f_{j}-E[Y]$

Taking the variance of Eq. (1) and exploiting the orthogonality between terms
in the expansion yields the following variance decomposition:

$\displaystyle\text{Var}(Y)$
$\displaystyle=\sum_{i=1}^{d}\text{Var}_{X_{i}}(E_{X_{\sim
i}}[Y|X_{i}])+\sum_{i=1}^{d}\sum_{j>i}\text{Var}_{X_{ij}}(E_{X_{\sim
ij}}[Y|X_{i},X_{j}])+\dots$ (3) $\displaystyle V(Y)$
$\displaystyle=\sum_{i=1}^{d}V_{i}+\sum_{i=1}^{d}\sum_{j>i}V_{ij}+\dots$

where $V_{i}$ represent the main effect variances and $V_{ij}$ represent the
variance contribution from interactions between variables $X_{i}$ and $X_{j}$.
The Sobol sensitivity indices reflect the relative contribution of each term
in Eq. (3) to the total variance and are given by:

$S_{i}=\dfrac{V_{i}}{V(Y)}=\dfrac{\text{Var}_{X_{i}}(E_{X_{\sim
i}}[Y|X_{i}])}{\text{Var}(Y)}$ (4)
$S_{ij}=\dfrac{V_{ij}}{V(Y)}=\dfrac{\text{Var}_{X_{ij}}(E_{X_{\sim
ij}}[Y|X_{i},X_{j}])-\text{Var}_{X_{i}}(E_{X_{\sim
i}}[Y|X_{i}])-\text{Var}_{X_{j}}(E_{X_{\sim j}}[Y|X_{j}])}{\text{Var}(Y)}$ (5)

where $S_{i}$ are the main effect sensitivity indices and $S_{ij}$ are the
interaction sensitivities and we have that

$\sum_{i=1}^{d}S_{i}+\sum_{i=1}^{d}\sum_{j>i}S_{ij}+\dots=1$ (6)

Here, we will focus predominantly on the main effect sensitivity indices in
Eq. (4), with some supporting discussion of interaction sensitivities.
Importantly, we highlight that the sensitivity index $S_{i}$ is defined as a
ratio of variances (the so-called main effect variance in the numerator and
total variance in the denominator), which both must be estimated in an active
learning framework.

###  2.2 Gaussian Process Regression / Kriging

Gaussian Process (GP) regression is a machine learning (ML) method that seeks
to identify the best fit Gaussian stochastic process to a set of data points
[44]. Kriging [45, 46, 47], a variation of GP regression, is a special case
wherein the GP regressor serves to interpolate between the set of training
points. Kriging is a widely-used ML method to construct approximate surrogate
models for complex computational models. It is particularly attractive for
active learning because it provides a model predictor as well as a measure of
the uncertainty in the prediction of the Gaussian posterior standard
deviation.

Consider a computational model $y(\bm{x})$ having input vector
$\bm{x}=\\{x_{1},x_{2},\dots,x_{d}\\}\in\mathbb{R}^{d}$. Next, consider that
we have $n$ realizations of $\bm{x}$, denoted by
$\bm{X}=\\{\bm{x}^{(1)},\bm{x}^{(2)},\dots,\bm{x}^{(n)}\\}$ and corresponding
model outputs $\bm{Y}=\\{y^{(1)},y^{(2)},\dots,y^{(n)}\\}$ evaluated by
$y^{(i)}=y(\bm{x}^{(i)}$). The pairs, $\bm{X},\bm{Y}$ serve as the training
data for the Gaussian process regression wherein the model $y(\bm{x})$ is
approximated by

$\displaystyle\mathcal{Y}(\bm{x},\omega)=\mathcal{F}(\bm{x})+Z(\bm{x},\omega)$
(7)

where $\mathcal{F}(\cdot)$ is a regression model and $Z(\cdot)$ is a zero-mean
Gaussian random process having sample space indexed by $\omega\in\Omega$. We
consider that the regression model $\mathcal{F}(\cdot)$ is defined through a
linear combination of basis functions
$\bm{f}(\bm{x})=\\{f_{1}(\bm{x}),f_{2}(\bm{x}),\dots,f_{p}(\bm{x})\\}$ having
coefficients $\bm{\beta}=\\{\beta_{1},\beta_{2},\dots,\beta_{p}\\}$ as

$\mathcal{F}(x)=\bm{\beta}^{T}\bm{f}(\bm{x})$ (8)

The Gaussian random process $Z(\cdot)$ is considered to have zero mean and
covariance:

$\displaystyle\mathbb{E}[Z(\bm{x}_{1})Z(\bm{x}_{2})]=\sigma^{2}_{z}\mathcal{R}(\bm{x}_{1},\bm{x}_{2}|\bm{\theta})$
(9)

where $\mathcal{R}(\cdot)$ is the autocorrelation function having
hyperparameters $\bm{\theta}$. In this work, we will consider the Gaussian
correlation model having form

$\mathcal{R}(\bm{x}_{1},\bm{x}_{2}|\bm{\theta})=\prod_{k=1}^{d}\exp(-\theta_{k}|x_{2k}-x_{1k}|^{2})=\prod_{k=1}^{d}r_{k}(x_{1k},x_{2k}|\theta_{k}).$
(10)

The GP is then fit to the normalized training data by identifying the
posterior distribution of the random process given the data. Given the
Gaussian assumption, the posterior distribution is determined by considering
that the joint distribution of the GP predictions $\mathcal{Y}(\bm{x})$ and
the observations $\bm{Y}$ is a multivariate Gaussian and can be expressed as

$\left\\{\begin{array}[]{c}\mathcal{Y}(\bm{x})\\\
\bm{Y}\end{array}\right\\}\sim
N\left(\left\\{\begin{array}[]{c}\bm{f}(\bm{x})^{T}\bm{\beta}\\\
\bm{F}\bm{\beta}\end{array}\right\\},\sigma_{z}^{2}\left\\{\begin{array}[]{cc}1&\bm{r}^{T}(\bm{x})\\\
\bm{r}(\bm{x})&\bm{R}\end{array}\right\\}\right)$ (11)

where $\bm{F}$ is the matrix of basis function evaluations at the training
points given by $F_{ij}=f_{j}(\bm{x}^{(i)}),i=1,\dots,n,j=1\dots p$,
$\bm{r}(\bm{x})$ is the vector of correlations between the prediction point
$\bm{x}$ and the training points $\bm{x}^{(i)}$ given by
$\bm{r}_{i}=\mathcal{R}(\bm{x},\bm{x}^{(i)}|\bm{\theta}),i=1,\dots,n$, and
$\bm{R}$ is the correlation matrix of points in the training set given by
$R_{ij}=\mathcal{R}(\bm{x}^{(i)},\bm{x}^{(j)}|\bm{\theta}),i,j=1,\dots,n$.

This posterior can be estimated by identifying the hyperparameters
$\bm{\theta}$, variance $\sigma_{z}$, and the regression coefficients
$\bm{\beta}$. This is typically done by solving for the parameters that either
maximize the likelihood of the observations or minimize the cross-validation
error. We apply maximum likelihood estimation where the likelihood function
follows from the Gaussian assumption as

$\mathcal{L}(\bm{\theta},\sigma_{z},\bm{\beta}|\bm{Y})=\dfrac{(\det\bm{R})^{-1/2}}{(2\pi\sigma_{z}^{2})^{n/2}}\exp\left[-\dfrac{1}{2\sigma_{z}^{2}}(\bm{Y}-\bm{F}\bm{\beta})^{T}\bm{R}^{-1}(\bm{Y}-\bm{F}\bm{\beta})\right]$
(12)

After estimating the hyperparameters, the mean and the variance of the
Gaussian random variable conditioned on the training data are given by:

$\hat{y}(\bm{x})=\bm{f}(\bm{x})^{T}\bm{\beta}+\bm{r}(\bm{x})^{T}\bm{R}^{-1}(\bm{Y}-\bm{F}\bm{\beta})$
(13)
$\sigma^{2}_{\hat{y}}(\bm{x})=\sigma^{2}_{z}\left(1-\bm{r}(\bm{x})^{T}\bm{R}^{-1}\bm{r}(\bm{x})+\bm{t}(\bm{x})^{T}(\bm{F}^{T}\bm{R}^{-1}\bm{F})^{-1}\bm{t}(\bm{x})\right)$
(14)

where

$\bm{t}(\bm{x})=\bm{F}^{T}\bm{R}^{-1}\bm{r}(\bm{x})-\bm{f}(\bm{x})$ (15)

These provide the Kriging prediction at a new point $\hat{y}(\bm{x})$ along
with a measure of its uncertainty $\sigma_{\hat{y}}(\bm{x})$.

##  3 Estimating Sobol Indices from Gaussian Process Regression

In this section, we review the process through which Sobol sensitivity indices
can be estimated from a GP regression model as originally derived in [48].
Marrel et al. [48] presented two approaches to compute Sobol indices using
GPs. In the first approach, the Sobol indices are computed simply using the GP
predictor, in which case the method provides a deterministic point estimator
based purely on the conditional mean as

$S_{i}=\dfrac{\text{Var}_{X_{i}}(E_{X_{\sim
i}}[E[\mathcal{Y}(X,\omega)]|X_{i}])}{\text{Var}(E[\mathcal{Y}(X,\omega)])}=\dfrac{\text{Var}_{X_{i}}(E_{X_{\sim
i}}[\hat{y}(X)|X_{i}])}{\text{Var}(\hat{y}(X))}$ (16)

where $\hat{y}(X)$ denotes the conditional mean in Eq. (13) and is treated as
a deterministic function in this approach.

The second approach, which we will employ in this study, uses the full GP
model as a stochastic function, which results in random variable sensitivity
indices given by:

$\tilde{S}_{i}(\omega)=\dfrac{\text{Var}_{X_{i}}(E_{X_{\sim
i}}[\mathcal{Y}(X,\omega)|X_{i}])}{E[\text{Var}(\mathcal{Y}(X,\omega))]}$ (17)

having mean and variance given by:

$\mu_{\tilde{S}_{i}}=\dfrac{E[\text{Var}_{X_{i}}(E_{X_{\sim
i}}[\mathcal{Y}(X,\omega)|X_{i}])]}{E[\text{Var}(\mathcal{Y}(X,\omega))]}$
(18)

and

$\sigma^{2}_{\tilde{S}_{i}}=\dfrac{\text{Var}(\text{Var}_{X_{i}}(E_{X_{\sim
i}}[\mathcal{Y}(X,\omega)|X_{i}]))}{(E[\text{Var}(\mathcal{Y}(X,\omega))])^{2}}$
(19)

where $\mu_{\tilde{S}_{i}}$ provides an estimate of the Sobol index and
$\sigma^{2}_{\tilde{S}_{i}}$ gives an estimate of its uncertainty. To compute
the first-order Sobol indices, we start by taking the expectation of the GP
$\mathcal{Y}(X,\omega)$ over all the inputs except ‘$X_{i}$’ in the numerator
of Eq. (17) and denote this as

$A_{i}(X_{i},\omega)=E_{X_{\sim i}}[\mathcal{Y}(X,\omega)|X_{i}]$ (20)

Since, $\mathcal{Y}(X,\omega)$ is a Gaussian process and expectation is a
linear operator, $A_{i}(X_{i},\omega)$ is also a GP referred to as the main
effect GP for $i$th input dimension.

The mean and covariance function of the main effect GP can be determined by
integrating the original GP with respect to the joint probability measure over
all inputs except $X_{i}$. Considering independent inputs, the mean function
is given by [49]:

$\displaystyle\mu_{A_{i}}(x_{i})=E[A_{i}(X_{i})]=\int_{\bm{x}_{\sim
i}}\hat{y}(\bm{x})\prod_{j\neq i}p_{X_{j}}(x_{j})dx_{j}$ (21)

which is easily computed as a product of one-dimensional integrals. Again
considering independent inputs, the covariance function of the main effect GP
is given by [49]:

$\displaystyle\text{Cov}(A_{i}(X_{1i}),A_{i}(X_{2i}))$
$\displaystyle=\int_{\bm{x}_{1\sim i}}\int_{\bm{x}_{2\sim
i}}\text{Cov}(\mathcal{Y}(\bm{x}_{1}),\mathcal{Y}(\bm{x}_{2}))\prod_{j\neq
i}p_{X_{j}}(x_{1j})d{x_{1j}}\prod_{k\neq i}p_{X_{k}}(x_{2k})d{x_{2k}}$ (22)

which can be expressed as a product of 2-dimensional integrals. Detailed
derivations for computing the mean and covariance functions of the main effect
GPs in closed form under specific conditions are provided in Appendix A. We
further note that Eqs. (21) and (22) can be generalized for interactions
between variables as shown in Appendix B. Computation of the interaction
sensitivities follows directly from the following for main effect
sensitivities. These interaction sensitivities can, in principle, be used
within the learning although this is outside the scope of the present work.

To compute Sobol indices, we need to compute the variance of the main effect
GP as

$\sigma^{2}_{A_{i}}=\int_{X_{i}}(A_{i}(x_{i})-E[A_{i}(X_{i})])^{2}p_{X_{i}}(x_{i})dx_{i}$
(23)

Here, a Monte Carlo estimate of this integral is computed by defining a random
discretization of the sample space of $X_{i}$ as
$\\{a_{1},a_{2},\dots,a_{n_{j}}\\}$ (note these points will be reused later
for learning function evaluations as well) and constructing the following
Gaussian random vector to approximate the main effect GP:

$\displaystyle V_{i}=[A_{i}(a_{1}),A_{i}(a_{2}),\ldots,A_{i}(a_{n_{j}})]^{T}$

We can then develop the following two Monte Carlo estimators of the Sobol
indices.

###  3.1 SI computation using the mean of the Main Effect GP

In the first approach, the mean of vector $V_{i}$ (i.e.
$\mu_{V_{i}}=E_{X_{i}}[V_{i}]$) is computed using Eqs. (21), where each
element of the mean vector is defined as
$\mu_{V_{i,j}}=\mu_{A_{i}}(a_{j})=E[A_{i}(a_{j})]$ $\forall$
$j\in\\{1,2,\ldots,n_{j}\\}$. Then, the following standard Monte Carlo
estimators are used to estimate the variance of the main effect GP

$\hat{\sigma}_{A_{i}}^{2}=\dfrac{1}{n_{j}-1}\sum_{j=1}^{n_{j}}(\mu_{V_{i,j}}-\bar{\mu}_{V_{i}})^{2}$
(24)

where $\bar{\mu}_{V_{i}}=\dfrac{1}{n_{j}}\sum_{j=1}^{n_{j}}\mu_{V_{i,j}}$ is
the average of the mean vector. The first-order Sobol indices are then
estimated by:

$\displaystyle\mu_{S_{i}}=\frac{\hat{\sigma}_{A_{i}}^{2}}{\hat{\sigma}^{2}}\quad\forall\quad
i\in\\{1,2,\ldots,d\\}$ (25)

where $\hat{\sigma}^{2}$ is the total variance of $\mathcal{Y}(\bm{x})$. This
estimator provides the Sobol index as identified from the mean GP predictor
for the main effect, but does not account for uncertainty/variability in the
main effect GP. In the next section, an estimator for the standard deviation
of the Sobol index is also given.

###  3.2 SI computation using the Main Effect GP

Here, we simulate the main effect GP and use these simulations to estimate
both the mean and the standard deviation of the Sobol index. First, the mean
vector and covariance matrix of $V_{i}$ are computed using Eqs. (21) and (22),
respectively. We then generate $n_{s}$ realizations of the random vector
$V_{i}$ using the Cholesky decomposition of the covariance matrix [46],
$C=LL^{T}$, where the $k^{th}$ realization for $i^{th}$ input (i.e. $V_{ik}$)
is given as follows

$\displaystyle
V_{ik}=E[A_{i}(X_{i})]+L\epsilon_{k}=[A_{ik}(a_{1}),A_{ik}(a_{2})),\ldots,A_{ik}(a_{n_{j}})]^{T}$
(26)

where $\epsilon_{k}$ is a realization of an $n_{j}$ dimensional zero-mean,
uncorrelated Gaussian random vector. Then, the following standard Monte Carlo
estimators are used to estimate the variance of the main effect GP

$\hat{\sigma}_{A_{i}}^{2}=\dfrac{1}{n_{s}}\sum_{k=1}^{n_{s}}\hat{\sigma}_{A_{ik}}^{2}=\dfrac{1}{n_{s}}\sum_{k=1}^{n_{s}}\dfrac{1}{n_{j}-1}\sum_{j=1}^{n_{j}}(A_{ik}(a_{j})-\hat{\mu}_{V_{ik}})^{2}$
(27)

where $\hat{\mu}_{V_{ik}}$ is the sample mean for vector $V_{ik}$,
$\hat{\sigma}_{A_{ik}}^{2}$ is the $k^{th}$ estimate of the variance of main
effect GP ($A_{i}$), and the corresponding standard error is given by

$\hat{s}_{A_{i}}^{2}=\dfrac{1}{n_{s}-1}\sum_{k=1}^{n_{s}}(\hat{\sigma}_{A_{ik}}^{2}-\hat{\sigma}_{A_{i}}^{2})^{2}.$
(28)

From these, the mean and standard deviation of the first-order Sobol indices
are estimated as:

$\displaystyle\mu_{S_{i}}=\frac{\hat{\sigma}_{A_{i}}^{2}}{\hat{\sigma}^{2}}\quad\text{and}\quad\sigma^{2}_{S_{i}}=\frac{\hat{s}_{A_{i}}^{2}}{\hat{\sigma}^{2}}\quad\forall\quad
i\in\\{1,2,\ldots,d\\}$ (29)

where $\hat{\sigma}^{2}$ is the total variance of $\mathcal{Y}(\bm{x})$.

###  3.3 Computational Considerations

The computational cost of applying these two methods for Sobol index
estimation are vastly different and can play a critical role when integrated
into an active learning loop where estimates need to be computed repeatedly.
Moreover, we have to repeat these estimates for each input dimension, which
can be problematic for problems with a large number of inputs. This must be
weighed against the accuracy of the two approaches and the need to estimate
the standard deviation of the Sobol index. The first approach only utilizes
the mean predictor of Main Effect GP, which is computationally inexpensive but
doesn’t estimate the uncertainty in Sobol indices. The second approach
considers the complete probability structure of the Main Effect GP and
therefore provides an error measure for each Sobol index. But it comes at a
high computational cost because it requires a large set of sample points
($n_{j}$ should be large) to compute Sobol indices with high accuracy. This
results in large covariance matrices, which can make the main effect
simulations in Eq. (26) expensive. We therefore employ the first approach to
estimate Sobol indices because we do not require uncertainty estimates on the
Sobol indices at each iteration.

##  4 Active Learning for Sobol Index Estimation

Next, we compare strategies to perform active learning for Sobol Index
estimation. We specifically consider two different learning functions designed
with the GP-based Sobol index estimates from Eq. (17) in mind. We first
recognize that the Sobol indices given by Eq. (17) are expressed as a ratio of
two variances. The numerator is the variance of the main effect GP and can be
expressed as $\sigma_{A_{i}}^{2}=\text{Var}(A_{i}(X_{i}))$, while the
denominator is the expected value of the variance of the complete GP,
$\mathcal{Y}(X,\omega)$. The challenge of active learning in this case is
therefore to minimize uncertainty in both components. To do so is not
straightforward. We therefore propose two learning strategies that aim to
minimize uncertainty in either the numerator or the denominator.

According to the first strategy, we aim to minimize uncertainty in the overall
GP surrogate (Sobol index denominator). This can be achieved using two
existing, well-known learning functions, (a) Expected Improvement for Global
Fit (EIGF - described in section 4.1 – which has been demonstrated to be
highly effective for this task [22]) and (b) Variance Improvement for Global
Fit (VIGF - described in section 4.2). For the second strategy, we develop a
new learning function, referred to as the MUSIC (Minimize Uncertainty in Sobol
Index Computation) learning function, designed to select samples that reduce
uncertainty in the main effect GPs $A_{i}(X_{i})$ and described in Section
4.3.

Figure 1: Flowchart of the active learning strategies for global sensitivity
analysis. Each of the learning functions (and random sampling) shown in green
are studied in this work.

These two learning strategies are integrated into a standard active learning
loop illustrated in Figure 1, which operates as follows. First, a set of
initial training data $(\bm{X},\bm{Y})$ are generated by randomly sampling the
input random vector $\bm{x}$ and evaluating the computational model
$y(\bm{x})$ for each sample. Then, an initial GP surrogate model is fit to the
data and initial estimates of the Sobol indices can be made using the
estimators above. Next, we generate a large set of $n_{j}$ candidate samples
of the input random variables (typically $n_{j}\geq 10000$) and evaluate the
appropriate learning function at each candidate sample. According to the
specified learning function, a candidate sample, say $\bm{x}^{*}$, is selected
and the model is evaluated at this point, $y(\bm{x}^{*})$. A new GP is fit to
the augmented data set and the Sobol indices are updated. A convergence
criterion on the Sobol indices is assessed. If the Sobol indices are deemed
sufficiently accurate, the process stops. If the indices are not considered
converged, then a new set of candidate points are generated and the iterations
continue.

We compare both of these strategies with a simple random sampling approach in
which, rather than employing a learning function in each iteration, one of the
$n_{j}$ candidate samples is selected at random. This provides a baseline to
measure whether the active learning strategies are effective for improving the
convergence of Sobol index estimates.

###  4.1 Expected Improvement for Global Fit (EIGF)

To improve the fit of the global GP, $\mathcal{Y}(X,\omega))$, and
specifically minimize uncertainty in the estimate of its variance, we employ
the EIGF. Proposed by Lam [22], the EIGF defines an improvement function given
by:

$I(\bm{x})=(\mathcal{Y}(\bm{x})-y(\bm{x}^{(j^{*})}))^{2}$ (30)

where $y(\bm{x}^{(j^{*})})$ is the observation of $y$ at the point
$\bm{x}^{(j^{*})}$ that is nearest to the point $\bm{x}$ and
$\mathcal{Y}(\bm{x})$ is the normal random variable resulting from evaluating
the GP at point $\bm{x}$. Taking the expectation of Eq. (30) yields

$\text{EIGF}(\bm{x})=E[I(\bm{x})]=(\hat{y}(\bm{x})-y(\bm{x}^{(j^{*})}))^{2}+\sigma_{\hat{y}}^{2}(\bm{x})$
(31)

where, $\hat{y}(\bm{x})$ and $\sigma_{\hat{y}}^{2}(\bm{x})$ are the GP mean
and variance at point $\bm{x}$ (see Eqs. (13) and (14)). The EIGF learning
function aims to identify the point $\bm{x}^{*}$ that maximizes the
expectation in Eq. (31) as:

$\bm{x}^{*}=\operatorname*{argmax}_{\bm{x}}(\text{EIGF}(\bm{x}))$

This point is selected as the new sample point for model evaluation and
retraining of the GP surrogate model in the active learning loop shown in
Figure 1.

A closer look at the EIGF in Eq. (31) demonstrates a balance of two factors –
often referred to as exploration and exploitation. In the first term,
$(\hat{y}(\bm{x})-y(\bm{x}^{(j^{*})}))^{2}$, we see that the EIGF targets
regions where variability in the function is large, which corresponds to
exploitation. This is balanced by the second term,
$\sigma_{\hat{y}}^{2}(\bm{x})$, which corresponds to regions where the GP
prediction has high uncertainty, likely caused by a lack of training data in
this region and encouraging exploration. However, it is well-known that the
EIGF learning function favors exploitation (i.e. the first term often
dominates) and may not explore new regions of parameter space sufficiently.
Therefore a more balanced approach may be considered.

###  4.2 Variance of Improvement for Global Fit (VIGF)

Mohammadi and Challenor [50] proposed a learning function using same
improvement function as EIGF function, given in Eq. (30). They further
recognized that $I(x)/\sigma^{2}(x)$ can be characterised by a non-central
chi-square distribution ([50]).

$I(x)/\sigma^{2}(x)\sim\chi^{{}^{\prime}2}\bigg{(}\kappa=1,\lambda=\bigg{(}\frac{\mu(x)-y(\bm{x}^{(j^{*})})}{\sigma(x)}\bigg{)}\bigg{)}$
(32)

where the number of degree of freedom ($\kappa$) and noncentrality parameter
($\lambda$) are defined using the GP’s posterior mean and standard deviation.
Thus, the variance of the improvement function is expressed as

$\displaystyle\text{VIGF}(x)=Var\\{I(x)\\}=4\sigma^{2}(x)\left[(\mu(x)-y(\bm{x}^{(j^{*})}))^{2}+2\sigma^{2}(x)\right].$
(33)

In the VIGF, the first term ($4\sigma^{2}(x)(\mu(x)-y(\bm{x}^{(j^{*})}))^{2}$)
is a product of terms that contribute toward local exploitation and global
exploration, thus favoring both. The second term ($8\sigma^{4}(x)$) focuses
exclusively on global exploration. This learning function therefore gives more
weight to global exploration than EIGF and may provide a better balance in the
exploration/exploitation trade off.

###  4.3 MUSIC Learning Function

From Eq. (25) and (29), we see that the quality of the Sobol index estimates
depends strongly on the quality of the main effect GPs. Large uncertainties in
the variance of the main effect GPs will result in associated large
uncertainties in Sobol index estimates. We therefore aim to create a training
set that produces the most accurate main effect GPs. This problem, as we will
see, can be posed in active learning terms as identifying the set of points
that yields the best global fit for the main effect GPs, taking inspiration
from the EIGF [22] and VIGF [50].

Let us first consider that we can construct an improvement function for the
main effect GPs in each dimension, similar to Eq. (30), as follows:

$I_{i}(x_{i})=(A_{i}(x_{i})-E_{X_{\sim
i}}[Y|X_{i}](x^{(j^{*})}_{i}))^{2},\quad i\in\\{1,\dots,d\\}$ (34)

where $x_{i}$, $x^{(j^{*})}_{i}$ are the $i^{th}$ component of points $\bm{x}$
and $\bm{x}^{(j^{*})}$, respectively and $x^{(j^{*})}_{i}$ is the nearest
training sample to $x_{i}$ in $i^{\text{th}}$ dimension. We immediately
recognize that we cannot actually observe the true main effect $E_{X_{\sim
i}}[Y(x^{(j^{*})}_{i})|X_{i}]$. Instead, we must estimate this value using the
main effect GP yielding a component-wise improvement function given by

$I_{A_{i}}(x_{i})=(A_{i}(x_{i})-A_{i}(x^{(j^{*})}_{i}))^{2},\quad
i\in\\{1,\dots,d\\}$ (35)

Taking the expectation and variance of Eq. (35) yields:

$E[I_{A_{i}}(x_{i})]=(\mu_{A_{i}}(x_{i})-\mu_{A_{i}}(x^{(j^{*})}_{i}))^{2}+\sigma^{2}_{A_{i}}(x_{i})$
(36)
$V[I_{A_{i}}(x_{i})]=4\sigma^{2}_{A_{i}}(x_{i})[(\mu_{A_{i}}(x_{i})-\mu_{A_{i}}(x^{(j^{*})}_{i}))^{2}+2\sigma^{2}_{A_{i}}(x_{i})]$
(37)

where $\mu_{A_{i}}(x_{i})$, $\mu_{A_{i}}(x^{(j^{*})}_{i})$ can be computed
using Eq. (21) and $\sigma^{2}_{A_{i}}(x_{i})$ can be computed from the
diagonal terms of Eq. (22). A closer look at $E[I_{A_{i}}(x_{i})]$
demonstrates that it also balances exploration and exploitation specifically
with regard to resolving the main effects. The first term targets regions
where the main effects vary strongly with $X_{i}$ (exploitation) and the
second term targets regions where the main effect GP has high uncertainty due
to a lack of training data (exploration). Eq. (37) has the same interpretation
as the VIGF, balancing exploration and exploitation, but operating on the main
effects.

Notice that the improvement criteria in Eq. (36) and (37) focus on improving
the prediction of individual main effect GPs. These learning criteria can be
collectively represented by a vector as follows:

$\displaystyle\bm{M_{A}}(\bm{x})=\begin{bmatrix}E[I_{A_{1}}(x_{1})]\\\
E[I_{A_{2}}(x_{2})]\\\ \vdots\\\
E[I_{A_{d}}(x_{d})\end{bmatrix}\quad\text{or}\quad\begin{bmatrix}V[I_{A_{1}}(x_{1})]\\\
V[I_{A_{2}}(x_{2})]\\\ \vdots\\\ V[I_{A_{d}}(x_{d})\end{bmatrix}$ (38)

From Eq. (38), we could directly derive improvement-based learning functions
over individual main effects. This might be useful, for example, if we were
interested in resolving a specific term in the HDMR in Eq. (1). However, we
are not interested only in individual main effects. Rather, we are interested
in ensuring convergence of all main effects, to the extent possible. This
means we must somehow balance the $d$ improvement functions that result from
Eq. (38) to establish in single improvement function. We do so by introducing
a distance-based pre-factor that combines individual learning criteria on main
effect GPs. This results in a new learning function, expressed as

$\displaystyle\text{\Large$\eighthnote$}(\bm{x})=\bm{D}^{T}(\bm{x})\bm{M_{A}}(\bm{x})$
(39)

and termed the MUSIC (Minimizing Uncertainty in Sobol Index Convergence)
learning function. We can consider numerous different distance pre-factors.
Here, we specifically propose the following two pre-factors that produce
convex combinations of main effect improvements

$\displaystyle\bm{D}_{1}(\bm{x})=\bm{W}\cdot\begin{bmatrix}|x_{1}-x^{(j^{*})}_{1}|\\\
|x_{2}-x^{(j^{*})}_{2}|\\\ \vdots\\\
|x_{d}-x^{(j^{*})}_{d}|\end{bmatrix}\quad\text{and}\quad\bm{D}_{2}(\bm{x})=\bm{W}\cdot||\bm{x}-\bm{x}^{(j^{*})}||_{2}^{2}\bm{1}_{d}$
(40)

where $||\bm{x}-\bm{x}^{(j^{*})}||_{2}$ is the Euclidean distance between the
sample $\bm{x}$ and nearest point $\bm{x}^{(j^{*})}$, and
$|x_{i}-x^{(j^{*})}_{i}|$ is the distance between the $i^{\text{th}}$
component of $\bm{x}$ and nearest point $\bm{x}^{(j^{*})}$ along
$i^{\text{th}}$ dimension. In either function, given two points whose total
main effect improvements are comparable the algorithm should favor the point
that is furthest from its nearest neighbor either in a Euclidean sense or in a
component-wise sense. The weights ($\bm{W}=[w_{1},w_{2},\ldots,w_{d}]^{T}$)
can be arbitrarily selected such that $\sum_{i}w_{i}=1$ to favor improvement
in particular main effects. For example, these can be set equal to the
existing Sobol index estimates in order to favor exploration in directions
that correspond to higher Sobol indices.

Using this learning function, we then select the next sample as the one that
maximizes it as:

$\displaystyle\bm{x}^{*}=\operatorname*{argmax}_{\bm{x}}\Large\eighthnote(\bm{x})$

from among the set of $n_{j}$ candidate samples. Table 1 summarizes the
learning function for four combinations based on the improvement criteria and
the distance-based pre-factor.

S.No. | Combinations | Learning Function
---|---|---
1 | MUSIC+EIGF D1 | $\Large\eighthnote_{E_{1}}(\bm{x})=\sum_{i=1}^{d}w_{i}|x_{i}-x^{(j^{*})}_{i}|E[I_{A_{i}}(x_{i})]$
2 | MUSIC+EIGF D2 | $\Large\eighthnote_{E_{2}}(\bm{x})=||\bm{x}-\bm{x}^{(j^{*})}||_{2}^{2}\sum_{i=1}^{d}w_{i}E[I_{A_{i}}(x_{i})]$
3 | MUSIC+VIGF D1 | $\Large\eighthnote_{V_{1}}(\bm{x})=\sum_{i=1}^{d}w_{i}|x_{i}-x^{(j^{*})}_{i}|V[I_{A_{i}}(x_{i})]$
4 | MUSIC+VIGF D2 | $\Large\eighthnote_{V_{2}}(\bm{x})=||\bm{x}-\bm{x}^{(j^{*})}||_{2}^{2}\sum_{i=1}^{d}w_{i}V[I_{A_{i}}(x_{i})]$
Table 1: Four different combinations of MUSIC learning strategy based on
improvement and distance-based pre-factor.

##  5 Exploration of Convergence using Analytical Functions

In this section, we study the convergence of the proposed active learning
methods for four benchmark analytical functions with known Sobol indices and
having increasing dimensions ranging from 2 to 15. Prior to presenting these
examples, we provide a simple basis upon which to explain the observed
convergence by considering the ratio of two quantities with some error.

###  5.1 Ratios of Errors

The dependence on the absolute error of a ratio on estimates of the numerator
and denominator is not straightforward. As we will see, reducing error in the
numerator, denominator, or both does not necessarily reduce error in the
ratio. Consider the following ratio:

$\displaystyle S=\frac{A}{Y}$

Next, consider that we have estimates of $S,A$ and $Y$ given by
$\hat{S},\hat{A}$ and $\hat{Y}$. Define the absolute error in the numerator
($A$) and denominator ($Y$) as:

$\displaystyle\Delta{A}=|A-\hat{A}|$ $\displaystyle\Delta{Y}=|Y-\hat{Y}|$

such that

$\displaystyle\hat{A}=A\pm\Delta{A}$ $\displaystyle\hat{Y}=Y\pm\Delta{Y}$

Next, the absolute error in the ratio can be expressed as:

$\displaystyle\Delta{S}=|S-\hat{S}|=\Bigg{|}S-\frac{\hat{A}}{\hat{Y}}\Bigg{|}$
$\displaystyle\Delta{S}=\Bigg{|}S-\frac{A\pm\Delta{A}}{Y\pm\Delta{Y}}\Bigg{|}$
$\displaystyle\Delta{S}=\Bigg{|}\frac{S(Y\pm\Delta{Y})-A\mp\Delta{A}}{Y\pm\Delta{Y}}\Bigg{|}$
$\displaystyle\Delta{S}=\Bigg{|}\frac{SY\pm
S\Delta{Y}-A\mp\Delta{A}}{Y\pm\Delta{Y}}\Bigg{|}$
$\displaystyle\Delta{S}=\Bigg{|}\frac{\pm
S\Delta{Y}\mp\Delta{A}}{Y\pm\Delta{Y}}\Bigg{|},$

which depends, in a nontrivial way, on the errors $\Delta A$ and $\Delta Y$ –
but also depends on the value of $S$ itself.

Recognizing that our GP models for $A$ and $Y$ presented in the formulation of
the GP-based Sobol indices above may be systematically biased (depending on
the training set and learning function being used) and that we may not know
this bias a priori, we explore this error under four specific cases.

  * •

Case 1: $\hat{A}$ and $\hat{Y}$ are both underestimated such that
$\hat{A}=A-\Delta{A}$ and $\hat{Y}=Y-\Delta{Y}$ yielding

$\displaystyle\Delta{S}=\frac{|-S\Delta{Y}+\Delta{A}|}{Y-\Delta{Y}}$

  * •

Case 2: $\hat{A}$ and $\hat{Y}$ are both overestimated such that
$\hat{A}=A+\Delta{A}$ and $\hat{Y}=Y+\Delta{Y}$ yielding

$\displaystyle\Delta{S}=\frac{|S\Delta{Y}-\Delta{A}|}{Y+\Delta{Y}}$

  * •

Case 3: $\hat{A}$ is overestimated ($\hat{A}=A+\Delta{A}$) and $\hat{Y}$ is
underestimated ($\hat{Y}=Y-\Delta{Y}$) yielding

$\displaystyle\Delta{S}=\frac{|-S\Delta{Y}-\Delta{A}|}{Y-\Delta{Y}}$

  * •

Case 4: $\hat{A}$ is underestimated ($\hat{A}=A-\Delta{A}$) and $\hat{Y}$ is
overestimated ($\hat{Y}=Y+\Delta{Y}$) yielding

$\displaystyle\Delta{S}=\frac{|S\Delta{Y}+\Delta{A}|}{Y+\Delta{Y}}$

Figures 2 and 3 show plots of the dependence of $\Delta S$ on $\Delta A$ and
$\Delta Y$ for $S=0.01$ and $S=0.8$, respectively, in each of these four cases
for a hypothetical ratio. Notice that, for $S=0.01$ (Figure 2) the errors
appear to decrease toward zero as both $\Delta A$ and $\Delta Y$ decrease
111In fact, Cases 1 and 2 have the same bias observed in the $S=0.8$ case but
the influence is negligible because $S$ is very small.. That is, if either (or
both) $\Delta A$ and $\Delta Y$ decrease, then $\Delta S$ decreases. However,
for $S=0.8$, a decrease in either $\Delta A$ or $\Delta Y$ (or both) does not
guarantee a decrease in $\Delta S$. In Cases 1 and 2, regardless of $\Delta A$
(even for very high $\Delta A$), it is possible to obtain zero $\Delta S$. In
fact, for a fixed value of $\Delta Y$, as $\Delta A$ decreases $\Delta S$ will
decrease to zero and then begin increasing again. This occurs because the
systematic bias in both Case 1 and Case 2 can cause the ratio $S=A/Y$ to be
correct even if both $A$ and $Y$ are (perhaps severely) incorrect. Indeed, two
points $X_{1}$ and $X_{2}$ are shown in Figure Y where $\Delta
A_{X_{1}}>\Delta A_{X_{1}}$ and $\Delta Y_{X_{1}}>\Delta Y_{X_{2}}$, yet
$\Delta S_{X_{1}}<\Delta S_{X_{2}}$. As we will see, this systematic bias in
$A$ and $Y$ can make accelerated convergence in the Sobol indices difficult.

Figure 2: Influence of errors in the numerator, $\Delta{A}$, and denominator,
$\Delta{Y}$, on error $\Delta{S}$ in the ratio $S=A/Y$ for $S=0.01$ for four
cases of systematic bias (a) Case 1: Both $\Delta{A}$ and $\Delta{Y}$ are
underestimated, (b) Case 2: Both $\Delta{A}$ and $\Delta{Y}$ are
overestimated, (c) $\Delta{A}$ is overestimated and $\Delta{Y}$ is
underestimated, and (d) $\Delta{A}$ is underestimated and $\Delta{Y}$ is
overestimated. Figure 3: Influence of errors in the numerator, $\Delta{A}$,
and denominator, $\Delta{Y}$, on error $\Delta{S}$ in the ratio $S=A/Y$ for
$S=0.8$ for four cases of systematic bias (a) Case 1: Both $\Delta{A}$ and
$\Delta{Y}$ are underestimated, (b) Case 2: Both $\Delta{A}$ and $\Delta{Y}$
are overestimated, (c) $\Delta{A}$ is overestimated and $\Delta{Y}$ is
underestimated, and (d) $\Delta{A}$ is underestimated and $\Delta{Y}$ is
overestimated.

Next, we will study the convergence, and specifically the behavior of errors
in light of the discussion above, for the different active learning schemes
using four analytical test cases of increasing dimension.

###  5.2 Square Exponential Function

First, we consider the following square exponential function of two random
inputs (i.e. $\mathbf{X}=[X_{1},X_{2}]$).

$\displaystyle y=X_{1}\exp{(-X_{1}^{2}-X_{2}^{2})}$

where, $X_{1}$ and $X_{2}$ are uniformly distributed between [a, b]. Two cases
are considered in this study, one with tight bounds around the two peaks
($a=-2$ and $b=2$) and another with larger bounds $b$ ($a=-2$ and $b=6$) such
that a large flat region is observed across much of the domain. Both cases are
illustrated in Figure 4. The presence of this flat region has a significant
impact on the Sobol indices and, moreover, affects the convergence of the
active learning algorithms. In both cases, the Sobol indices can be computed
analytically. For the first case ($b=2$), the first-order Sobol indices are
$S_{1}=0.6208$ and $S_{2}=0$. For the second case, the Sobol indices are
$S_{1}=0.3119$ and $S_{2}=2.30\times 10^{-5}$.

Figure 4: Square exponential functions (a) with tight bounds and no flat
region and (b) with wide bounds and large flat region. Figure 5: Exponential
Function ($b=2$, no flat region) – Mean square convergence with $\sigma$
confidence intervals from 100 repeated trials of (a) output variance
$\sigma_{Y}^{2}$, (b,c) main effect variances $\sigma_{A_{1}}^{2}$ and
$\sigma_{A_{2}}^{2}$, and (d,e) Sobol indices $S_{1}$ and $S_{2}$ for MUSIC
and EIGF active learning schemes compared with random sampling. Figure 6:
Exponential Function ($b=2$, no flat region) – Mean square convergence with
$\sigma$ confidence intervals from 100 repeated trials of (a) output variance
$\sigma_{Y}^{2}$, (b,c) main effect variances $\sigma_{A_{1}}^{2}$ and
$\sigma_{A_{2}}^{2}$, and (d,e) Sobol indices $S_{1}$ and $S_{2}$ for MUSIC
and VIGF active learning schemes compared with random sampling.

Figure 5 shows the mean square error convergence from 100 repeated trials of
the different active learning schemes using the EIGF and MUSIC learning
functions compared with random sampling for the case where $b=2$ (no flat
region). Figure 5(a) shows convergence of the variance of the output,
$\sigma^{2}_{Y}$, as estimated by the GP. We see that the EIGF generally
outperforms the other two approaches, converging to a very accurate estimate
for $\sigma^{2}_{Y}$ with less than $\sim 30$ samples. Figures 5(b) and (c)
show the same convergence plots for the variance of the main effect GPs
$\sigma_{A_{1}}^{2}$ and $\sigma_{A_{2}}^{2}$, respectively. We see that, in
main effects the MUSIC learning function offers comparable convergence to the
the EIGF in $A_{1}$ and superior convergence in $A_{2}$. Both methods
generally outperform random sampling in the main effect variances.

Finally, figures 5(d) and (e) show convergence of the Sobol indices $S_{1}$
and $S_{2}$. For this problem, the Sobol convergence tracks closely with the
main effect convergence, but with slight differences. For example, we recall
that the variance $\sigma_{A_{1}}^{2}$ from Figure 5(b) is the numerator and
$\sigma_{Y}^{2}$ from Figure 5(a) is the denominator of the Sobol index
$S_{1}$. However, despite superior performance in both $\sigma_{A_{1}}^{2}$
and $\sigma_{Y}^{2}$, the MUSIC function performs mildly worse than random
sampling in $S_{1}$ for small samples. This is due to the ‘ratio of errors’
issues discussed in the previous section.

Figure 6 shows the same convergence plots but compares VIGF and MUSIC (based
on VIGF) learning functions against the sequential random sampling. Similar to
EIGF, the VIGF learning function results in the best estimate for output
variance ($\sigma_{Y}^{2}$), as shown in Fig. 6(a). However, VIGF outperforms
EIGF for main effect and sobol estimates ($\sigma_{A_{i}}^{2}$) and yields
similar convergence to the ‘MUSIC+VIGF D2’ combination. This improvement is
due to the better prediction at both functional features (peak and valley).
Consider, the main effect of first dimension $E[Y|X_{1}]$ for a fixed value of
$X_{1}$ (imagine a cross-section at arbitrary $X_{1}$ in Figure 4), the main
effect function depends on either the peak or the valley (depending on the
value of $X_{1}$). But, the main effect of the second dimension ($E[Y|X_{2}]$)
depends on both features for all values of $X_{2}$. Since VIGF gives more
importance to exploration (as compared to EIGF), it resolves both features
simultaneously – it does not get stuck exploiting one feature. This results in
better estimates for $\sigma_{Y}^{2}$ as well as $\sigma_{A_{i}}^{2}$ and
eventually better Sobol estimates as shown in 6(e).

Figure 7 and 8 show the same convergence plots for the case with $b=6$ (having
a flat region). As expected, the EIGF/VIGF learning functions perform much
better than the other methods for $\sigma_{Y}^{2}$, as they both perform well
for tasks where exploitation is paramount [51]. They both also exhibit
superior performance in the main effect variances $\sigma_{A_{1}}^{2}$ and
$\sigma_{A_{2}}^{2}$ and Sobol indices $S_{1}$ and $S_{2}$ for small sample
sizes. Again, this is a result of its strong exploitation properties. The
MUSIC function meanwhile provides superior convergence to random sampling in
all cases as well.

Figure 7: Exponential Function ($b=6$, with flat region) – Mean square
convergence with $\sigma$ confidence intervals from 100 repeated trials of (a)
output variance $\sigma_{Y}^{2}$, (b,c) main effect variances
$\sigma_{A_{1}}^{2}$ and $\sigma_{A_{2}}^{2}$, and (d,e) Sobol indices $S_{1}$
and $S_{2}$ for MUSIC and EIGF active learning schemes compared with random
sampling. Figure 8: Exponential Function ($b=6$, with flat region) – Mean
square convergence with $\sigma$ confidence intervals from 100 repeated trials
of (a) output variance $\sigma_{Y}^{2}$, (b,c) main effect variances
$\sigma_{A_{1}}^{2}$ and $\sigma_{A_{2}}^{2}$, and (d,e) Sobol indices $S_{1}$
and $S_{2}$ for MUSIC and VIGF active learning schemes compared with random
sampling.

###  5.3 Ishigami function

Next, consider the well-known Ishigmai function with three independent input
random variables ($\mathbf{X}=[X_{1},X_{2},X_{3}]$) that are uniformly
distributed over [$-\pi,\pi$], and defined by:

$\displaystyle
y(\mathbf{x})=\sin{x_{1}}+a\sin^{2}{x_{2}}+bx_{3}^{4}\sin{x_{1}}$

where parameters $a=7$, and $b=0.1$. This function has nonlinear and
nonmonotonic behavior, but the third input dimension exhibits a peculiar
behavior. The first-order contribution of the third dimension is zero towards
output variance, but it can’t be completely ignored as an interaction
sensitivity index (i.e. $S_{13}$) has a significant impact on output variance.
The first-order Sobol indices are $S_{1}=0.3139$, $S_{2}=0.4424$, and
$S_{3}=0$, determined analytically by solving for the variance of the main
effect and output as

$\displaystyle
Var\\{A_{1}\\}=\frac{1}{2}\Big{(}1+\frac{b\pi^{4}}{5}\Big{)}^{2},\quad
Var\\{A_{2}\\}=\frac{a^{2}}{8},\quad Var\\{A_{3}\\}=0$
$Var\\{Y\\}=\frac{a^{2}}{8}+\frac{b\pi^{4}}{5}+\frac{b^{2}\pi^{8}}{18}+\frac{1}{2}$

The performance of EIGF and MUSIC (‘EIGF D1’ and ‘EIGF D2’) adaptive sampling
strategies is shown in Figure. 9. The adaptive algorithms start with a small
set of 10 samples generated using Latin Hypercube Sampling and convergence is
observed up to 500 samples. After updating the surrogate model, Sobol indices
are estimated using 25,000 candidate points. The process is repeated for 100
trials to obtain confidence intervals in the convergence. In this example, the
EIGF struggles to accurately estimate the output variance ($\sigma^{2}_{Y}$ in
Figure 9(a)) and the main effect variances ($\sigma_{A_{i}}^{2}$ in Figure
9(b)–(d)). This can be explained by the sinusoidal nature of the Ishigami
function, which requires a balance between exploration and exploitation. Error
convergence of the ‘MUSIC+EIGF D1’ learning function is similar to random
sampling for every estimate. ‘MUSIC+EIGF D2’, on the other hand, leads to
superior convergence in nearly all aspects. The distance-based pre-factor is
the main difference here. The $D_{2}$ function uses the Euclidean distance and
ensures that training samples are globally far away from each other, which
helps the ‘MUSIC+EIGF’ learning function to avoid local exploitation. Because
the $D_{1}$ function utilizes an element-wise distance component, it does not
avoid local exploitation as effectively. Figures 9 (e-g) show the ultimate
convergence of the Sobol index estimates where again ‘MUSIC+EIGF D2’ generally
outperforms the other criteria.

Figure 9: Ishigami Function – Mean square convergence with $\sigma$ confidence
intervals from 100 repeated trials of (a) output variance $\sigma_{Y}^{2}$,
(b-d) main effect variances $\sigma_{A_{i}}^{2}$ $\forall$ $i\in\\{1,2,3\\}$,
and (e-g) Sobol indices $S_{1}$, $S_{2}$, and $S_{3}$ for MUSIC and EIGF
active learning schemes compared with random sampling.

Figure 10 likewise compares the VIGF and MUSIC (+VIGF) learning functions to
random sampling. The main difference here is that the VIGF improves
performance considerably over the EIGF. Given this consistent improvement with
the VIGF, we will compare only the results using VIGF in the remaining
analytical examples.

Figure 10: Ishigami Function – Mean square convergence of (a) output variance
$\sigma_{Y}^{2}$, (b-d) main effect variances $\sigma_{A_{i}}^{2}$ $\forall$
$i\in\\{1,2,3\\}$, and (e-g) Sobol indices $S_{1}$, $S_{2}$, and $S_{3}$ for
MUSIC and VIGF active learning schemes compared with random sampling.

Here, the convergence of the first and second Sobol indices may appear
counterintuitive. The MUSIC function estimates both the total and main effect
variances with high accuracy, but the error in the Sobol estimates is
comparable to random sampling. This is again due to the effect of ratios of
errors. Figure. 11 shows the convergence path of the Sobol indices by plotting
traces of $(\hat{\sigma}^{2}_{Y},\hat{\sigma}^{2}_{A_{i}})$ averaged over 100
trials. The grey surface represents the Sobol index for any specific pair. The
green star is the actual first-order Sobol index ($S_{i}$) value for the $i$th
dimension. The grey line shows all combinations of
$(\hat{\sigma}^{2}_{Y},\hat{\sigma}^{2}_{A_{i}})$, such that their ratio is
equal to the true Sobol index (i.e. $\hat{S}_{i}=S_{i}$). The black, orange,
and blue curves illustrate the convergence of random, EIGF, and MUSIC
sampling, respectively. We notice that, although MUSIC and EIGF perform better
in the variance estimators (see convergence above), the path that their
estimators takes toward the true value does not follow the grey line. That is,
their Sobol index estimates are biased due to the ratio of errors. This is not
the case for random sampling, which converges largely along the gray line.
This counter-intuitive behavior is not observed for the third dimension
because the actual Sobol index is zero (i.e. $S_{3}=0$) and therefore depends
only on the main effect GP. As illustrated in Figure 2, when $S_{i}$ is very
small convergence in the numerator will yield convergence in the ratio. We can
therefore conclude that the MUSIC function will be very effective at
identifying small Sobol indices.

Figure 11: Dependence of Sobol indices on the variance of main effect GPs and
Output GP.

###  5.4 G-function

Next, we consider the 5-dimensional G-function, defined as the following
product of non-linear functions in each dimension,

$\displaystyle
y(\mathbf{x})=\prod_{k=1}^{d}\frac{|4x^{(k)}-2|+a^{(k)}}{a^{(k)}+1}$

where each input feature is uniformly distributed between [0, 1] (i.e.
$X^{(k)}\sim\text{Uniform}(0,1)$). The analytical Sobol indices in five
dimensions (i.e. $d=5$) are 0.48, 0.21, 0.12, 0.08, 0.06, corresponding to
$a^{(k)}=k\text{, }k\in{1,2,3,4,5}$ from the following equations

$\displaystyle Var\\{A_{i}\\}=\frac{1}{3(1+a^{(i)})^{2}}\quad\quad
Var\\{Y\\}=\Big{[}\prod_{k=1}^{d}(1+Var\\{A_{k}\\})\Big{]}-1$ $\displaystyle
S_{i}=\frac{1/(3(1+a^{(i)})^{2})}{\Big{[}\prod_{k=1}^{d}{[1+1/(3(1+a^{(k)})^{2})}]\Big{]}-1}$

Again we compare the three sampling strategies, where each sampling strategy
initiates by training a GP on 30 samples obtained from Latin Hypercube
Sampling and adaptively selects new samples up to 500 samples. The Sobol
indices are computed using 25000 candidate points. Figure 12 shows the
convergence for VIGF and MUSIC (based on VIGF) compared to random sampling.
Here, the ‘MUSIC+VIGF’ strategy and the VIGF show superior performance to
random sampling in the total variance (Figure 12(a)), main effects (Figure
12(b)–(f)), and Sobol indices ((Figure 12(g)–(k)). Results from EIGF learning
are not shown due to poor convergence in comparison to VIGF.

Figure 12: G-Sobol Function – Mean square convergence with $\sigma$ confidence
intervals from 100 repeated trials of (a) output variance $\sigma_{Y}^{2}$,
(b–f) main effect variances $\sigma_{A_{i}}^{2}$ $\forall$
$i\in\\{1,2,3,4,5\\}$, and (g–k) Sobol indices $S_{i}$ for MUSIC and VIGF
active learning schemes compared with random sampling.

###  5.5 Gaussian Function

To conclude, we compare performance for a high input dimensional problem. A
synthetic $15$d function is considered here, such that 75.5% of the
contribution is from the first-order effects of three inputs. The function is
defined as a product of 15 independent square exponential functions, similar
to a Gaussian kernel given by

$\displaystyle
y(\bm{x})=\prod_{i=1}^{d}\exp{\Bigg{(}-\frac{x_{i}^{2}}{a_{i}}\Bigg{)}}$

where, $x_{i}\sim U(-3,3)$ $\forall$ $i\in\\{1,2,\ldots,d\\}$ and

$\bm{a}=[1.45,3.3,15,50,55,58,59,100,102,112.5,150,160,180,190,200]$

Figure 13(a) shows the main effect contribution and Figure 13(b) shows the
true Sobol indices (in log scale).

Figure 13: Gaussian function: (a) First-order main effect functions and (b)
True Sobol Indices.

Figure 14 shows convergence of the Sobol indices for the first three
dimensions using the VIGF and MUSIC with VIGF learning functions compared to
random sampling for 500 samples. We notice that the VIGF and ‘MUSIC+VIGF D2’
perform quite poorly for this example. In particular, the ‘MUSIC+VIGF D2’
exhibits poor convergence because the Euclidean distance is not a good measure
of distance in 15 dimensions [52]. When this is replaced with the component-
wise distance measure, convergence is much better and is comparable with
random sampling. Convergence for dimensions 4–15 are not shown here because
they immediately produce very low errors in all cases given that the
corresponding Sobol indices are near zero. Importantly, we observe that no
learning function is capable of producing superior convergence to random
sampling in high-dimensions.

Figure 14: Gaussian Function – Mean square convergence with $\sigma$
confidence intervals from 50 repeated trials for Sobol indices (a) $S_{1}$,
(b) $S_{2}$, and (c) $S_{3}$ for MUSIC and VIGF active learning schemes
compared with random sampling.

##  6 Application to Boundary Layer Wind Tunnel Experimental Design

The original motivation for developing this active learning scheme for Sobol
index estimation came from an experimental investigation into the influence of
roughness terrain parameters in large-scale Boundary Layer Wind Tunnel (BLWT)
experiments. In this study, experiments were performed using the University of
Florida BLWT (UF-BLWT illustrated in Figure 15), which features a novel
automated roughness element grid called the Terraformer. This Terraformer is
capable of rapidly and automatically changing the roughness terrain at the
push of a button to mimic a wide-range of wind flow conditions [53].

Figure 15: Schematic of the University of Florida Boundary Layer Wind Tunnel
(BLWT). The Terraformer roughness element grid is capable of rapidly and
automatically reconfiguring its 1116 elements to produce a wide-range of wind
flow characteristics.

In this work, we parameterize the Terraformer (which has 1116 independent
elements capable of changing both height and profile orientation) by a random
field in the along-wind direction using a truncated Karhunen-Loeve (KL)
expansion such that the element at a location $x$ along the Terraformer is
given by

$\displaystyle
h(x,\bm{\theta})=\mu_{H}+\sum_{i=1}^{d}\sqrt{\lambda^{(i)}}{\theta}^{(i)}f^{(i)}(x)$

where $\mu_{H}=80$ mm is the mean element height, $\theta^{(i)}$ are standard
normal random variables, and $\lambda^{(i)}$ and $f^{(i)}(x)$ are the
eigenvalues and eigenfunctions of the covariance function given by

$C(x_{1},x_{2})=\sigma^{2}\exp(-a|x_{2}-x_{1}|)\cos(\omega|x_{2}-x_{1}|)$ (41)

Here $\sigma^{2}=100$ mm is the variance of the field, $\omega$ is the wave
number, and $a$ and $\omega$ are selected such that the length scale of the
covariance function is given by $L=a/(a^{2}2+\omega^{2})=3000$ mm, or
approximately 1/6 of the Terraformer length. The KL expansion is truncated to
$d=10$ terms.

The significant expense in this project is due to experimental time, where the
execution time for a single experiment is approximately 20 minutes. Even with
the high degree of automation afforded by the UF-BLWT, we cannot afford to
explore the influence of all $d=10$ random field parameters on the resulting
wind flow. We therefore aim to adaptively perform global sensitivity analysis
to identify which random field parameters have the most influence on the wind
characteristic of interest and reduce the dimension of the random field to
facilitate a feasible set of experiments. Our specific characteristic of
interest is the turbulence intensity profile of the resulting wind flow
defined at a height above the floor ($z$) by

$I_{u}(z)=\dfrac{\sigma_{u}(z)}{\mu_{u}(z)}$ (42)

where $\mu_{u}(z)$ and $\sigma_{u}(z)$ are the mean and standard deviation of
the longitudinal component of the wind velocity $u(z)$ – which is a random
process. More specifically, we aim to quantify the sensitivity of the
following distance measure to the 10 KL random variables:

$d(\bm{\theta})=||\bm{I}_{u}(\bm{\theta})-\bm{I}_{u}^{*}||_{2}$ (43)

where $\bm{I}_{u}(\bm{\theta})$ is the vector of discretized turbulence
intensities observed at points separated by 20 mm in the range $z=[180,500]$
mm for a Terraformer configuration with parameters $\bm{\theta}$ and
$\bm{I}_{u}^{*}$ is the same vector defined for a reference roughness grid
with all element uniformly set to $h=80$ mm. In other words, we aim to
quantify the sensitivity of the deviation of the turbulence intensity profile
from the reference profile to each of the KL random variables.

For this investigation, we applied an earlier component-wise application of
the proposed EIGF-based MUSIC learning function not demonstrated in detail
above. This initial version of MUSIC focuses on only one dimension at a time
and chooses a new sample to minimize the uncertainty in the estimate of the
particular Sobol index associated with the largest MUSIC component. More
specifically, using the expected main effect improvement defined in Eq. (36),
this version selects the dimension with maximum improvement in the individual
main effect GPs as

$\displaystyle x_{i}^{*}=\operatorname*{argmax}_{x_{i}}{E[I_{A_{i}}(x_{i})]}$
(44)

where $x_{i}^{*}$ is the $i^{\text{th}}$ element of the new sample
$\bm{x^{*}}$ and the other elements are selected randomly over the input
marginal distribution. Although this version of the MUSIC function was
ultimately improved and the methods discussed in detail above were explored in
more detail, it was successfully used to adaptively estimate sensitivities for
the BLWT experimental investigation.

Initially, 30 experiments were performed with parameters $\bm{\theta}$
selected by Latin HyperCube sampling. Then, 270 more experiments were
conducted using the MUSIC learning strategy described above. Figure 16 shows
the final estimates (after 300 experiments) of the first-order Sobol indices,
along with the standard error in the estimates. The experiments show that the
distance metric between two roughness configurations is dominated by the
second eigen parameter. Furthermore, except for input dimensions 1 and 2, no
other input parameter has a considerable first-order effect on the distance.
Using the approach described in Appendix B, we were further able to make an
initial assessment of the interaction sensitivities and determine that
significant interactions between inputs 1, 2, and 3 play an important role in
the resulting turbulence intensity profiles. Using these insights, we were
able to reduce the dimension of the KL random field to $d=3$ and conduct a
more rigorous study, further details of which are provided in [1].

Figure 16: Final estimates of the first-order Sobol Indices.

##  7 Conclusions

This study investigates the application of active learning for global
sensitivity analysis. We specifically propose a novel new learning function,
termed the MUSIC (minimize uncertainty in Sobol index convergence) function,
that leverages the main effect variances in Gaussian Process-based Sensitivity
Indices, i.e. the numerator in Sobol index definition. We compare this
strategy with a more traditional active learning strategy aimed at accurately
capturing output variance $\sigma^{2}_{Y}$, i.e. the denominator of the Sobol
index, using the existing EIGF and VIGF learning functions. Both strategies
are further compared with simple sequential random sampling. We provide a
careful analysis of the convergence for four analytical test functions and
study the convergence behavior. An important insight that arises is that
convergence in Sobol indices are difficult for active learning because they
are defined as a ratio and the learning strategies must focus on either the
numerator or the denominator and achieving rapid convergence in either, or
both, is not sufficient to ensure convergence in the ratio.

Summarizing, we can make the following observations:

  * •

In general, the VIGF function outperforms the EIGF, as does the MUSIC with
VIGF when compared to the MUSIC with EIGF.

  * •

The MUSIC learning functions are very successful at rapidly identifying small
Sobol indices.

  * •

For low-dimensional problems, the MUSIC with VIGF and distance pre-factor
based on a Euclidean distance generally outperforms the other methods in terms
of Sobol index convergence.

  * •

The Euclidean distance-based pre-factor ($\bm{D}_{2}$) is effective for low-
dimensional problems but does not work well in high-dimensions.

  * •

At best, the MUSIC learning strategy is comparable in performance to random
sampling for high-dimensional problems.

Finally, we provide a motivating application where active learning is applied
to learn global sensitivity indices for a set of large-scale boundary layer
wind tunnel experiments. This demonstrates a real, practical use case and
illustrates how it helped to significantly reduce experimental effort.

### Acknowledgements

This material is based upon work supported by the National Science Foundation
under Grant Nos. 1930389 and 1930625. Dr. Dimitrios Tsapetis has been
supported by Defense Threat Reduction Agency, Award HDTRA1202000.

## References

  * Shields et al. [2023] Michael D. Shields, Kurtis Gurley, Ryan Catarelli, Mohit Chauhan, Mariel Ojeda-Tuz, and Forrest Masters.  Active learning applied to automated physical systems increases the rate of discovery.  _Sci. Rep._ , 2023.
  * Shahriari et al. [2015] Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P Adams, and Nando De Freitas.  Taking the human out of the loop: A review of bayesian optimization.  _Proceedings of the IEEE_ , 104(1):148–175, 2015\.
  * Frazier [2018] Peter I Frazier.  Bayesian optimization.  In _Recent advances in optimization and modeling of contemporary problems_ , pages 255–278. Informs, 2018.
  * Greenhill et al. [2020] Stewart Greenhill, Santu Rana, Sunil Gupta, Pratibha Vellanki, and Svetha Venkatesh.  Bayesian optimization for adaptive experimental design: A review.  _IEEE access_ , 8:13937–13948, 2020.
  * Močkus [1975] Jonas Močkus.  On bayesian methods for seeking the extremum.  In _Optimization Techniques IFIP Technical Conference: Novosibirsk, July 1–7, 1974_ , pages 400–404. Springer, 1975.
  * Jones et al. [1998] Donald R Jones, Matthias Schonlau, and William J Welch.  Efficient global optimization of expensive black-box functions.  _Journal of Global optimization_ , 13(4):455, 1998\.
  * Huang et al. [2006] Deng Huang, Theodore T Allen, William I Notz, Ning Zeng, et al.  Global optimization of stochastic black-box systems via sequential kriging meta-models.  _Journal of global optimization_ , 34(3):441–466, 2006.
  * Frazier et al. [2009] Peter Frazier, Warren Powell, and Savas Dayanik.  The knowledge-gradient policy for correlated normal beliefs.  _INFORMS journal on Computing_ , 21(4):599–613, 2009.
  * Scott et al. [2011] Warren Scott, Peter Frazier, and Warren Powell.  The correlated knowledge gradient for simulation optimization of continuous parameters using gaussian process regression.  _SIAM Journal on Optimization_ , 21(3):996–1026, 2011.
  * Wu and Frazier [2016] Jian Wu and Peter Frazier.  The parallel knowledge gradient method for batch bayesian optimization.  _Advances in neural information processing systems_ , 29, 2016.
  * Hennig and Schuler [2012] Philipp Hennig and Christian J Schuler.  Entropy search for information-efficient global optimization.  _Journal of Machine Learning Research_ , 13(6), 2012.
  * Hernández-Lobato et al. [2014] José Miguel Hernández-Lobato, Matthew W Hoffman, and Zoubin Ghahramani.  Predictive entropy search for efficient global optimization of black-box functions.  _Advances in neural information processing systems_ , 27, 2014.
  * Bichon et al. [2008] Barron J Bichon, Michael S Eldred, Laura Painton Swiler, Sandaran Mahadevan, and John M McFarland.  Efficient global reliability analysis for nonlinear implicit performance functions.  _AIAA journal_ , 46(10):2459–2468, 2008.
  * Echard et al. [2011] Benjamin Echard, Nicolas Gayton, and Maurice Lemaire.  Ak-mcs: an active learning reliability method combining kriging and monte carlo simulation.  _Structural Safety_ , 33(2):145–154, 2011.
  * Huang et al. [2016] Xiaoxu Huang, Jianqiao Chen, and Hongping Zhu.  Assessing small failure probabilities by ak–ss: An active learning method combining kriging and subset simulation.  _Structural Safety_ , 59:86–95, 2016.
  * Zhang et al. [2019] Jinhao Zhang, Mi Xiao, and Liang Gao.  An active learning reliability method combining kriging constructed with exploration and exploitation of failure region and subset simulation.  _Reliability Engineering & System Safety_, 188:90–102, 2019.
  * Xu et al. [2020] Chunlong Xu, Weidong Chen, Jingxin Ma, Yaqin Shi, and Shengzhuo Lu.  Ak-mss: An adaptation of the ak-mcs method for small failure probabilities.  _Structural Safety_ , 86:101971, 2020.
  * Sundar and Shields [2019] VS Sundar and Michael D Shields.  Reliability analysis using adaptive kriging surrogates with multimodel inference.  _ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: Civil Engineering_ , 5(2):04019004, 2019.
  * Lelièvre et al. [2018] Nicolas Lelièvre, Pierre Beaurepaire, Cécile Mattrand, and Nicolas Gayton.  Ak-mcsi: A kriging-based method to deal with small failure probabilities and time-consuming models.  _Structural Safety_ , 73:1–11, 2018.
  * Razaaly and Congedo [2020] Nassim Razaaly and Pietro Marco Congedo.  Extension of ak-mcs for the efficient computation of very small failure probabilities.  _Reliability Engineering & System Safety_, 203:107084, 2020.
  * El Haj and Soubra [2021] Abdul-Kader El Haj and Abdul-Hamid Soubra.  Improved active learning probabilistic approach for the computation of failure probability.  _Structural Safety_ , 88:102011, 2021.
  * Lam [2008] Chen Quin Lam.  _Sequential adaptive designs in computer experiments for response surface model fit_.  PhD thesis, The Ohio State University, 2008.
  * Marelli and Sudret [2018] Stefano Marelli and Bruno Sudret.  An active-learning algorithm that combines sparse polynomial chaos expansions and bootstrap for structural reliability analysis.  _Structural Safety_ , 75:67–74, 2018.
  * Cheng and Lu [2020] Kai Cheng and Zhenzhou Lu.  Active learning polynomial chaos expansion for reliability analysis by maximizing expected indicator function prediction error.  _International Journal for Numerical Methods in Engineering_ , 121(14):3159–3177, 2020.
  * Novák et al. [2023] Lukáš Novák, Michael D Shields, Václav Sadílek, and Miroslav Vořechovskỳ.  Active learning-based domain adaptive localized polynomial chaos expansion.  _arXiv preprint arXiv:2301.13635_ , 2023.
  * Sener and Savarese [2017] Ozan Sener and Silvio Savarese.  Active learning for convolutional neural networks: A core-set approach.  _arXiv preprint arXiv:1708.00489_ , 2017.
  * Haut et al. [2018] Juan Mario Haut, Mercedes E Paoletti, Javier Plaza, Jun Li, and Antonio Plaza.  Active learning with convolutional neural networks for hyperspectral image classification using a new bayesian approach.  _IEEE Transactions on Geoscience and Remote Sensing_ , 56(11):6440–6461, 2018.
  * Schröder and Niekler [2020] Christopher Schröder and Andreas Niekler.  A survey of active learning for text classification using deep neural networks.  _arXiv preprint arXiv:2008.07267_ , 2020.
  * DeLoach [1998] Richard DeLoach.  Applications of modern experiment design to wind tunnel testing at nasa langley research center.  In _36th AIAA Aerospace Sciences Meeting and Exhibit_ , page 713, 1998\.
  * Hill et al. [2011] Raymond R Hill, Derek A Leggio, Shay R Capehart, and August G Roesener.  Examining improved experimental designs for wind tunnel testing using monte carlo sampling methods.  _Quality and Reliability Engineering International_ , 27(6):795–803, 2011.
  * VanDercreek et al. [2019] Colin P VanDercreek, Alireza Amiri-Simkooei, Mirjam Snellen, and Daniele Ragni.  Experimental design and stochastic modeling of hydrodynamic wave propagation within cavities for wind tunnel acoustic measurements.  _International Journal of Aeroacoustics_ , 18(8):752–779, 2019.
  * Lookman et al. [2019] Turab Lookman, Prasanna V Balachandran, Dezhen Xue, and Ruihao Yuan.  Active learning in materials science with emphasis on adaptive sampling using uncertainties for targeted design.  _npj Computational Materials_ , 5(1):21, 2019\.
  * Kusne et al. [2020] A Gilad Kusne, Heshan Yu, Changming Wu, Huairuo Zhang, Jason Hattrick-Simpers, Brian DeCost, Suchismita Sarker, Corey Oses, Cormac Toher, Stefano Curtarolo, et al.  On-the-fly closed-loop materials discovery via bayesian active learning.  _Nature communications_ , 11(1):5966, 2020.
  * Jablonka et al. [2021] Kevin Maik Jablonka, Giriprasad Melpatti Jothiappan, Shefang Wang, Berend Smit, and Brian Yoo.  Bias free multiobjective active learning for materials design and discovery.  _Nature communications_ , 12(1):2312, 2021.
  * Reker [2019] Daniel Reker.  Practical considerations for active machine learning in drug discovery.  _Drug Discovery Today: Technologies_ , 32-33:73–79, 2019\.  ISSN 1740-6749.  doi: https://doi.org/10.1016/j.ddtec.2020.06.001.  URL https://www.sciencedirect.com/science/article/pii/S1740674920300019.  Artificial Intelligence.
  * Warmuth et al. [2003] Manfred K Warmuth, Jun Liao, Gunnar Rätsch, Michael Mathieson, Santosh Putta, and Christian Lemmen.  Active learning with support vector machines in the drug discovery process.  _Journal of chemical information and computer sciences_ , 43(2):667–673, 2003.
  * Naik et al. [2016] Armaghan W Naik, Joshua D Kangas, Devin P Sullivan, and Robert F Murphy.  Active machine learning-driven experimentation to determine compound effects on protein patterns.  _Elife_ , 5:e10047, 2016.
  * Sverchkov and Craven [2017] Yuriy Sverchkov and Mark Craven.  A review of active learning approaches to experimental design for uncovering biological networks.  _PLoS computational biology_ , 13(6):e1005466, 2017.
  * Osugi et al. [2005] Thomas Osugi, Deng Kim, and Stephen Scott.  Balancing exploration and exploitation: A new algorithm for active machine learning.  In _Fifth IEEE International Conference on Data Mining (ICDM’05)_ , pages 8–pp. IEEE, 2005.
  * Shan and Wang [2010] Songqing Shan and G Gary Wang.  Survey of modeling and optimization strategies to solve high-dimensional design problems with computationally-expensive black-box functions.  _Structural and multidisciplinary optimization_ , 41:219–241, 2010.
  * Sobol [1993] Ilya M Sobol.  Sensitivity analysis for non-linear mathematical models.  _Math. Modeling Comput. Exp._ , 1:407–414, 1993.
  * Saltelli and Sobol [1995] Andrea Saltelli and Ilya M Sobol.  About the use of rank transformation in sensitivity analysis of model output.  _Reliability Engineering & System Safety_, 50(3):225–239, 1995.
  * Archer et al. [1997] GEB Archer, Andrea Saltelli, and Ilya Meyerovich Sobol.  Sensitivity measures, anova-like techniques and the use of bootstrap.  _Journal of Statistical Computation and Simulation_ , 58(2):99–120, 1997.
  * Rasmussen et al. [2006] Carl Edward Rasmussen, Christopher KI Williams, et al.  _Gaussian processes for machine learning_ , volume 1.  Springer, 2006.
  * Matheron [1963] Georges Matheron.  Principles of geostatistics.  _Economic geology_ , 58(8):1246–1266, 1963.
  * Cressie [2015] Noel Cressie.  _Statistics for spatial data_.  John Wiley & Sons, 2015.
  * Martin and Simpson [2005] Jay D Martin and Timothy W Simpson.  Use of kriging models to approximate deterministic computer models.  _AIAA journal_ , 43(4):853–863, 2005.
  * Marrel et al. [2009] Amandine Marrel, Bertrand Iooss, Béatrice Laurent, and Olivier Roustant.  Calculations of sobol indices for the gaussian process metamodel.  _Reliability Engineering & System Safety_, 94(3):742 – 751, 2009.  ISSN 0951-8320.  doi: https://doi.org/10.1016/j.ress.2008.07.008.  URL http://www.sciencedirect.com/science/article/pii/S0951832008001981.
  * Oakley JE [2004] O’Hagan A. Oakley JE.  Probabilistic sensitivity analysis of complex models: a bayesian approach.  _J R Stat Soc Series B_ , 66(3):751 – 769, 2004\.  ISSN 0951-8320.  doi: https://doi.org/10.1111/j.1467-9868.2004.05304.x.  URL https://rss.onlinelibrary.wiley.com/doi/full/10.1111/j.1467-9868.2004.05304.x?casa_token=Y0y5mWztnHkAAAAA%3AjzEXLY5nnEboGT3EcYWcLbT4GzEXIxHmn8vsax6SE62azxwgzNnIePcsFoE-MthGW0WwRr9KcR2z9Q.
  * Mohammadi and Challenor [2022] Hossein Mohammadi and Peter Challenor.  Sequential adaptive design for emulating costly computer codes.  _arXiv preprint arXiv:2206.12113_ , 2022.
  * Beck and Guillas [2016] Joakim Beck and Serge Guillas.  Sequential design with mutual information for computer experiments (mice): Emulation of a tsunami model.  _SIAM/ASA Journal on Uncertainty Quantification_ , 4(1):739–766, 2016.
  * Aggarwal et al. [2001] Charu C Aggarwal, Alexander Hinneburg, and Daniel A Keim.  On the surprising behavior of distance metrics in high dimensional space.  In _Database Theory—ICDT 2001: 8th International Conference London, UK, January 4–6, 2001 Proceedings 8_ , pages 420–434. Springer, 2001\.
  * Catarelli et al. [2020] Ryan A Catarelli, Pedro L Fernández-Cabán, Brian M Phillips, Jennifer A Bridge, Forrest J Masters, Kurtis R Gurley, and David O Prevatt.  Automation and new capabilities in the university of florida nheri boundary layer wind tunnel.  _Frontiers in Built Environment_ , 6:558151, 2020.
  * Lebrun and Dutfoy [2009] Regis Lebrun and Anne Dutfoy.  An innovating analysis of the nataf transformation from the copula viewpoint.  _Probabilistic Engineering Mechanics_ , 24(3):312–320, 2009.

##  Appendix A Appendix 1

In this appendix, we show how to compute the mean and covariance functions for
the main effect GPs under specific conditions. We then demonstrate that under
conditions of widely used correlation models, these functions can be evaluated
analytically without the need for numerical integrations.

###  A.1 Mean function of the main effect GPs

We begin by restating Eq. (21) as:

$\displaystyle E[A(X_{i})]=\int_{\bm{x}_{\sim
i}}\hat{y}({\mathbf{\bm{x}}})\prod_{j\neq i}p_{X_{j}}(x_{j})dx_{j}$

Next, we simply plug in the expression for $\hat{y}$ from Eq. (13), which
yields

$\displaystyle E[A(X_{i})]$ $\displaystyle=\int_{\bm{x}_{\sim
i}}(\bm{f}(\bm{x})^{T}\bm{\beta}+\bm{r}(\bm{x})^{T}\bm{R}^{-1}(\bm{Y}-\bm{F}\bm{\beta}))\prod_{j\neq
i}p_{X_{j}}(x_{j})dx_{j}$ (A.1) $\displaystyle=\bigg{[}\int_{\bm{x}_{\sim
i}}\bm{f}(\bm{x})^{T}\prod_{j\neq
i}p_{X_{j}}(x_{j})dx_{j}\bigg{]}\bm{\beta}+\bigg{[}\int_{\bm{x}_{\sim
i}}\bm{r}(\bm{x})^{T}\prod_{j\neq
i}p_{X_{j}}(x_{j})dx_{j}\bigg{]}\bm{R}^{-1}(\bm{Y}-\bm{F}\bm{\beta})$ (A.2)

Integrating the first term of Equation (A.2) for a linear basis
$\bm{f}(\bm{x})=\begin{bmatrix}1&x_{1}&x_{2}&\ldots&x_{d}\end{bmatrix}^{T}$.

$\displaystyle\int_{\bm{x}_{\sim i}}\bm{f}(\bm{x})^{T}\prod_{j\neq
i}p_{X_{j}}(x_{j})dx_{j}$ $\displaystyle=\int_{\bm{x}_{\sim
i}}\begin{bmatrix}1&x_{1}&x_{2}&\ldots&x_{d}\end{bmatrix}\prod_{j\neq
i}p_{X_{j}}(x_{j})dx_{j}$
$\displaystyle=\begin{bmatrix}1&E[X_{1}]&E[X_{2}]&\ldots&x_{i}&\ldots&E[X_{d}]\end{bmatrix}$

Integrating the second term of Equation (A.2) yields:

$\displaystyle\int_{\bm{x}_{\sim i}}\bm{r}(\bm{x})^{T}\prod_{j\neq
i}p_{X_{j}}(x_{j})dx_{j}=$ $\displaystyle\int_{\bm{x}_{\sim
i}}\begin{bmatrix}\mathcal{R}({\bm{x}},\bm{x}^{(1)})\\\
\mathcal{R}({\bm{x}},\bm{x}^{(2)})\\\ \vdots\\\
\mathcal{R}({\bm{x}},\bm{x}^{(n)})\end{bmatrix}^{T}\prod_{j\neq
i}p_{X_{j}}(x_{j})dx_{j}$ $\displaystyle=$ $\displaystyle\int_{\bm{x}_{\sim
i}}\begin{bmatrix}\prod_{k}r_{k}(x_{k},x_{k}^{(1)})\\\
\prod_{k}r_{k}(x_{k},x_{k}^{(2)})\\\ \vdots\\\
\prod_{k}r_{k}(x_{k},x_{k}^{(n)})\end{bmatrix}^{T}\prod_{j\neq
i}p_{X_{j}}(x_{j})dx_{j}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}r_{i}(x_{i},x_{i}^{(1)})\prod_{j\neq
i}\int_{x_{j}}r_{j}(x_{j},x_{j}^{(1)})p_{X_{j}}(x_{j})dx_{j}\\\
r_{i}(x_{i},x_{i}^{(2)})\prod_{j\neq
i}\int_{x_{j}}r_{j}(x_{j},x_{j}^{(2)})p_{X_{j}}(x_{j})dx_{j}\\\ \vdots\\\
r_{i}(x_{i},x_{i}^{(n)})\prod_{j\neq
i}\int_{x_{j}}r_{j}(x_{j},x_{j}^{(n)})p_{X_{j}}(x_{j})dx_{j}\end{bmatrix}^{T}$

where we recall that $\bm{x}^{(l)},\sim l=1,\dots,n$ are the training samples
and $r_{k}(x_{k},x_{k}^{(l)})$ is the $k^{th}$ univariate correlation function
expressed in Eq. (10). Finally, Eq. (A.2) can be expressed as:

$\displaystyle E[A(X_{i})]=\begin{bmatrix}1\\\ E[X_{1}]\\\ E[X_{2}]\\\
\vdots\\\ x_{i}\\\ \vdots\\\
E[X_{d}]]\end{bmatrix}^{T}\bm{\beta}+\begin{bmatrix}r_{i}(x_{i},x_{i}^{(1)})\prod_{j\neq
i}\int_{x_{j}}r_{j}(x_{j},x_{j}^{(1)})p_{X_{j}}(x_{j})dx_{j}\\\
r_{i}(x_{i},x_{i}^{(2)})\prod_{j\neq
i}\int_{x_{j}}r_{j}(x_{j},x_{j}^{(2)})p_{X_{j}}(x_{j})dx_{j}\\\ \vdots\\\
r_{i}(x_{i},x_{i}^{(n)})\prod_{j\neq
i}\int_{x_{j}}r_{j}(x_{j},x_{j}^{(n)})p_{X_{j}}(x_{j})dx_{j}\end{bmatrix}^{T}\bm{R}^{-1}(\bm{Y}-\bm{F}\bm{\beta})$
(A.3)

Finally, in Section A.3, we show that the product of 1D integrals can be
computed analytically for the Gaussian correlation model, which results in the
following closed form for the mean of the main effect GP for uniformly
distributed $X$:

$\displaystyle E[A$
$\displaystyle(X_{i})]=\begin{bmatrix}1&E[X_{1}]&E[X_{2}]&\ldots&x_{i}&\ldots&E[X_{d}]]\end{bmatrix}\bm{\beta}$
$\displaystyle+\begin{bmatrix}r_{i}(x_{i},x_{i}^{(1)})\prod_{j\neq
i}\frac{1}{b_{j}-a_{j}}\sqrt{\frac{\pi}{\theta_{j}}}[\Phi(\sqrt{2\theta_{j}}(b_{j}-{x_{j}^{(1)}}))-\Phi(\sqrt{2\theta_{j}}(a_{j}-{x_{j}^{(1)}}))]\\\
r_{i}(x_{i},x_{i}^{(2)})\prod_{j\neq
i}\frac{1}{b_{j}-a_{j}}\sqrt{\frac{\pi}{\theta_{j}}}[\Phi(\sqrt{2\theta_{j}}(b_{j}-{x_{j}^{(2)}}))-\Phi(\sqrt{2\theta_{j}}(a_{j}-{x_{j}^{(2)}}))]\\\
\vdots\\\ r_{i}(x_{i},x_{i}^{(n)})\prod_{j\neq
i}\frac{1}{b_{j}-a_{j}}\sqrt{\frac{\pi}{\theta_{j}}}[\Phi(\sqrt{2\theta_{j}}(b_{j}-{x_{j}^{(n)}}))-\Phi(\sqrt{2\theta_{j}}(a_{j}-{x_{j}^{(n)}}))]\end{bmatrix}^{T}\bm{R}^{-1}(\bm{Y}-\bm{F}\bm{\beta})$
(A.4)

Note that in cases where the GP is trained with noisy data the above equation
can be easily modified to incorporate the noise. We simply substitute the
correlation matrix, $\bm{R}$, in the above equation by $\bm{R_{n}}$ defined
as:

$\displaystyle\bm{R}_{n}=((1-\tau)\bm{R}+\tau\bm{I})$

where ${\tau=\sigma_{\epsilon}^{2}}/{\sigma_{z}^{2}}$ and
$\sigma_{\epsilon}^{2}$ is the variance due to the noise in the data and
$\sigma_{z}^{2}$ is the total variance of the GP.

###  A.2 Covariance function of the main effect GPs

We again begin by restating the expression for the covariance function of the
main effect GPs from Eq. (22):

$\displaystyle\text{Cov}(A(X_{1i}),A(X_{2i}))$
$\displaystyle=\int_{\bm{x}_{1\sim i}}\int_{\bm{x}_{2\sim
i}}\text{Cov}(\mathcal{Y}(\bm{x}_{1}),\mathcal{Y}(\bm{x}_{2}))\prod_{j\neq
i}p_{X_{j}}(x_{1j})d{x_{1j}}\prod_{k\neq i}p_{X_{k}}(x_{2k})d{x_{2k}}$

The covariance of the GP $\mathcal{Y}(\bm{x})$ can be expressed as:

$\displaystyle\text{Cov}(\mathcal{Y}(\bm{x}_{1}),\mathcal{Y}(\bm{x}_{2}))=\sigma^{2}_{z}(r(\bm{x}_{1},\bm{x}_{2})-\bm{r}(\bm{x}_{1})^{T}\bm{R}^{-1}\bm{r}(\bm{x}_{2})+\bm{t}(\bm{x}_{1})^{T}(\bm{F}^{T}\bm{R}^{-1}\bm{F})^{-1}\bm{t}(\bm{x}_{2})),$

which yields the following covariance for the main effect GP:

$\text{Cov}(A(X_{1i}),A(X_{2i}))=\int_{\bm{x}_{1\sim i}}\int_{\bm{x}_{2\sim
i}}\bigg{[}\sigma^{2}_{z}(r(\bm{x}_{1},\bm{x}_{2})-\bm{r}(\bm{x}_{1})^{T}\bm{R}^{-1}\bm{r}(\bm{x}_{2})+\bm{t}(\bm{x}_{1})^{T}(\bm{F}^{T}\bm{R}^{-1}\bm{F})^{-1}\bm{t}(\bm{x}_{2}))\bigg{]}\\\
\prod_{j\neq i}p_{X_{j}}(x_{1j})d{x_{1j}}\prod_{k\neq
i}p_{X_{k}}(x_{2k})d{x_{2k}}$ (A.5)

Expanding this, we can express the first term of Eq. (A.5) as:

$\displaystyle\sigma^{2}_{z}\int_{\bm{x}_{1\sim i}}\int_{\bm{x}_{2\sim
i}}r(\bm{x}_{1},\bm{x}_{2})\prod_{j\neq
i}p_{X_{j}}(x_{1j})d{x_{1j}}\prod_{k\neq i}p_{X_{k}}(x_{2k})d{x_{2k}}=$
$\displaystyle\sigma^{2}_{z}\int_{\bm{x}_{1\sim i}}\int_{\bm{x}_{2\sim
i}}\prod_{l=1}^{d}r_{l}(x_{1l},x_{2l})\prod_{j\neq i}$ $\displaystyle
p_{X_{j}}(x_{1j})d{x_{1j}}\prod_{k\neq i}p_{X_{k}}(x_{2k})d{x_{2k}}$
$\displaystyle=\sigma^{2}_{z}r_{i}(x_{1i},x_{2i})\prod_{l\neq
i}^{d}\int_{\bm{x}_{1\sim i}}\int_{\bm{x}_{2\sim i}}r_{l}(x_{1l}$
$\displaystyle,x_{2l})p_{X_{l}}(x_{1l})d{x_{1l}}p_{X_{l}}(x_{2l})d{x_{2l}}$

In Section A.3, we derive the closed form for this integral given a uniform
distribution and Gaussian correlation model, which yields:

$\sigma^{2}_{z}\int_{\bm{x}_{1\sim i}}\int_{\bm{x}_{2\sim
i}}r(\bm{x}_{1},\bm{x}_{2})\prod_{j\neq
i}p_{X_{j}}(x_{1j})d{x_{1j}}\prod_{k\neq i}p_{X_{k}}(x_{2k})d{x_{2k}}=\\\
\sigma^{2}_{z}r_{i}(x_{1i},x_{2i})\prod_{l\neq
i}^{d}\left(\frac{1}{b_{l}-a_{l}}\right)^{2}\sqrt{\frac{\pi}{\theta_{l}}}\times\\\
\left[a_{l}+b_{l}-2[a_{l}\Phi(\sqrt{2\theta_{l}}(b_{l}-a_{l}))+b_{l}\Phi(\sqrt{2\theta_{l}}(a_{l}-b_{l}))]-\frac{1}{\sqrt{\pi\theta_{l}}}[1-\exp\\{-\theta_{l}(b_{l}-a_{l})^{2}\\}]\right]$

Next, the second term of the Eq. (A.5) is:

$\displaystyle-\sigma^{2}_{z}\int_{\bm{x}_{1\sim i}}\int_{\bm{x}_{2\sim
i}}\bm{r}(\bm{x}_{1})^{T}\bm{R}^{-1}\bm{r}(\bm{x}_{2})\prod_{j\neq
i}p_{X_{j}}(x_{1j})d{x_{1j}}\prod_{k\neq i}p_{X_{k}}(x_{2k})d{x_{2k}}=$
$\displaystyle-\sigma^{2}_{z}\bigg{[}\int_{\bm{x}_{1\sim
i}}\bm{r}(\bm{x}_{1})^{T}\prod_{j\neq
i}p_{X_{j}}(x_{1j})d{x_{1j}}\bigg{]}\bm{R}^{-1}\bigg{[}\int_{\bm{x}_{2\sim
i}}\bm{r}(\bm{x}_{2})$ $\displaystyle\prod_{k\neq
i}p_{X_{k}}(x_{2k})d{x_{2k}}\bigg{]}$

where the integrals in brackets can be solved in closed form as derived in
Section A.1. Finally, the third term of Eq. (A.5) is:

$\displaystyle-\sigma^{2}_{z}\int_{\bm{x}_{1\sim i}}\int_{\bm{x}_{2\sim
i}}\bm{t}(\bm{x}_{1})^{T}(\bm{F}^{T}\bm{R}^{-1}\bm{F})^{-1}\bm{t}(\bm{x}_{2})\prod_{j\neq
i}p_{X_{j}}(x_{1j})d{x_{1j}}\prod_{j\neq i}p_{X_{j}}(x_{2j})d{x_{2j}}=$
$\displaystyle-\sigma^{2}_{z}\bigg{[}\int_{\bm{x}_{1\sim
i}}\bm{t}(\bm{x}_{1})^{T}\prod_{j\neq
i}p_{X_{j}}(x_{1j})d{x_{1j}}\bigg{]}(\bm{F}^{T}\bm{R}^{-1}\bm{F})^{-1}\bigg{[}\int_{\bm{x}_{2\sim
i}}\bm{t}(\bm{x}_{2})$ $\displaystyle\prod_{j\neq
i}p_{X_{j}}(x_{2j})d{x_{2j}}\bigg{]}$

where, $\hat{\sigma}^{2}_{z}$ is the variance of the GP, ${\mathcal{Y}}$, and
$\bm{t}(\bm{x})=\bm{F}^{T}\bm{R}^{-1}\bm{r}(\bm{x})-\bm{f}(\bm{x})$ as defined
in Eq. (15). Finally, the integral of $\bm{t}(\bm{x})$ over the joint
distribution in brackets can be easily obtained from the equations in Section
A.1, as follows

$\displaystyle\int_{\bm{x}_{1\sim i}}\bm{t}(\bm{x})^{T}\prod_{j\neq
i}p_{X_{j}}(x_{j})d{x_{j}}$ $\displaystyle=\int_{\bm{x}_{1\sim
i}}\Big{[}\bm{F}^{T}\bm{R}^{-1}\bm{r}(\bm{x})-\bm{f}(\bm{x})\Big{]}\prod_{j\neq
i}p_{X_{j}}(x_{j})d{x_{j}}$
$\displaystyle=\bm{F}^{T}\bm{R}^{-1}\int_{\bm{x}_{1\sim
i}}\bm{r}(\bm{x})\prod_{j\neq i}p_{X_{j}}(x_{j})d{x_{j}}-\int_{\bm{x}_{1\sim
i}}\bm{f}(\bm{x})\prod_{j\neq i}p_{X_{j}}(x_{j})d{x_{j}}$
$\displaystyle=\bm{F}^{T}\bm{R}^{-1}\begin{bmatrix}r_{i}(x_{i},x_{i}^{(1)})\prod_{j\neq
i}\int_{x_{j}}r_{j}(x_{j},x_{j}^{(1)})p_{X_{j}}(x_{j})dx_{j}\\\
r_{i}(x_{i},x_{i}^{(2)})\prod_{j\neq
i}\int_{x_{j}}r_{j}(x_{j},x_{j}^{(2)})p_{X_{j}}(x_{j})dx_{j}\\\ \vdots\\\
r_{i}(x_{i},x_{i}^{(n)})\prod_{j\neq
i}\int_{x_{j}}r_{j}(x_{j},x_{j}^{(n)})p_{X_{j}}(x_{j})dx_{j}\end{bmatrix}-\begin{bmatrix}1\\\
E[X_{1}]\\\ E[X_{2}]\\\ \vdots\\\ x_{i}\\\ \vdots\\\ E[X_{d}]]\end{bmatrix}$

The above expression can be further simplified for a given uniform
distribution and Gaussian correlation model as:

$\displaystyle\int_{\bm{x}_{1\sim i}}\bm{t}(\bm{x})^{T}$
$\displaystyle\prod_{j\neq i}p_{X_{j}}(x_{j})d{x_{j}}$
$\displaystyle=\bm{F}^{T}\bm{R}^{-1}\begin{bmatrix}r_{i}(x_{i},x_{i}^{(1)})\prod_{j\neq
i}\frac{1}{b_{j}-a_{j}}\sqrt{\frac{\pi}{\theta_{j}}}[\Phi(\sqrt{2\theta_{j}}(b_{j}-{x_{j}^{(1)}}))-\Phi(\sqrt{2\theta_{j}}(a_{j}-{x_{j}^{(1)}}))]\\\
r_{i}(x_{i},x_{i}^{(2)})\prod_{j\neq
i}\frac{1}{b_{j}-a_{j}}\sqrt{\frac{\pi}{\theta_{j}}}[\Phi(\sqrt{2\theta_{j}}(b_{j}-{x_{j}^{(2)}}))-\Phi(\sqrt{2\theta_{j}}(a_{j}-{x_{j}^{(2)}}))]\\\
\vdots\\\ r_{i}(x_{i},x_{i}^{(n)})\prod_{j\neq
i}\frac{1}{b_{j}-a_{j}}\sqrt{\frac{\pi}{\theta_{j}}}[\Phi(\sqrt{2\theta_{j}}(b_{j}-{x_{j}^{(n)}}))-\Phi(\sqrt{2\theta_{j}}(a_{j}-{x_{j}^{(n)}}))]\end{bmatrix}-\begin{bmatrix}1\\\
E[X_{1}]\\\ E[X_{2}]\\\ \vdots\\\ x_{i}\\\ \vdots\\\ E[X_{d}]]\end{bmatrix}$

Thus, the main effect GP, $A(X_{i})$, can be completely defined by the mean
function, $E[A(X_{i})]$, and covariance function, $\text{Cov}(A(X_{1i})$,
$A(X_{2i}))$, derived in Sections A.1 and A.2.

Again in case of noisy training data, we simply substitute the correlation
matrix, $\bm{R}$, in the above equation by $\bm{R_{n}}$ as defined below:

$\displaystyle\bm{R}_{n}=((1-\tau)\bm{R}+\tau\bm{I})$

where, ${\tau=\sigma_{\epsilon}^{2}}/{\sigma_{z}^{2}}$ and
$\sigma_{\epsilon}^{2}$ is the variance due to the noise in the data and
$\sigma_{z}^{2}$ is the total variance of the GP.

###  A.3 Closed-form integration of the Gaussian correlation model

Let $x_{1}$ and $x_{2}$ be univariate sample points drawn from a uniform
distribution over the range $(a,b)$, $X_{1},X_{2}\sim U(a,b)$, having Gaussian
correlation (from Eq. (10)) given by:

$\displaystyle r(x_{1},x_{2}|\theta)=\exp\\{-\theta(x_{1}-x_{2})^{2}\\}$

Note that, without loss of generality, the arbitrary sample points can be
transformed to a uniform distribution through a suitable isoprobabilistic
transformation (e.g. Nataf [54]). As discussed in the previous appendices, we
need to integrate this correlation model over the pdf $p_{X_{1}}(x_{1})$ as
follows:

$\displaystyle\int_{-\infty}^{\infty}r(\theta,x_{1},x_{2})p_{{X_{1}}}({x_{1}})dx_{1}=$
$\displaystyle\int_{a}^{b}\exp\\{-\theta({x_{1}}-{x_{2}})^{2}\\}\frac{1}{b-a}d{x_{1}}$
$\displaystyle=$
$\displaystyle\frac{1}{b-a}\int_{a}^{b}\exp\\{-\theta({x_{1}}-{x_{2}})^{2}\\}d{x_{1}}$
$\displaystyle=$
$\displaystyle\frac{1}{b-a}\sqrt{2\pi(1/2\theta)}\int_{a}^{b}\frac{1}{\sqrt{2\pi(1/2\theta)}}\exp\bigg{\\{}-\frac{1}{2}\frac{({x_{1}}-{x_{2}})^{2}}{(1/2\theta)}\bigg{\\}}d{x_{1}}$
$\displaystyle=$
$\displaystyle\frac{1}{b-a}\sqrt{\frac{\pi}{\theta}}\bigg{[}\Phi\bigg{(}\sqrt{2\theta}(b-{x_{2}})\bigg{)}-\Phi\bigg{(}\sqrt{2\theta}(a-{x_{2}})\bigg{)}\bigg{]}$

where, $\Phi(\cdot)$ is the standard normal cumulative distribution function.
Further, we need the integrate the above result with respect to ${x_{2}}$ as

$\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}r(\theta,{x_{1}},{x_{2}})$
$\displaystyle p_{{X_{1}}}({x_{1}})p_{{X_{2}}}({x_{2}})d{x_{1}}d{x_{2}}$ (A.6)
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\bigg{[}\frac{1}{b-a}\sqrt{\frac{\pi}{\theta}}\bigg{[}\Phi\bigg{(}\sqrt{2\theta}(b-{x_{2}})\bigg{)}-\Phi\bigg{(}\sqrt{2\theta}(a-{x_{2}})\bigg{)}\bigg{]}\bigg{]}p_{{X_{2}}}({x_{2}})d{x_{2}}$
$\displaystyle=$
$\displaystyle\bigg{(}\frac{1}{b-a}\bigg{)}^{2}\sqrt{\frac{\pi}{\theta}}\int_{a}^{b}\bigg{[}\Phi\bigg{(}\sqrt{2\theta}(b-{x_{2}})\bigg{)}-\Phi\bigg{(}\sqrt{2\theta}(a-{x_{2}})\bigg{)}\bigg{]}d{x_{2}}$

We see that the above integral has two identical components that differ only
by the constant $a$ or $b$. Let us introduce a dummy variable, $z$, which can
take values $a$ or $b$. Integration by parts yields

$\displaystyle\int_{a}^{b}(1)\times\Phi$
$\displaystyle\bigg{(}\sqrt{2\theta}(z-{x_{2}})\bigg{)}d{x_{2}}=\Phi\bigg{(}\sqrt{2\theta}(z-{s})\bigg{)}x_{2}\bigg{|}_{a}^{b}-\int_{a}^{b}{x_{2}}\phi\bigg{(}\sqrt{2\theta}(z-{x_{2}})\bigg{)}(-\sqrt{2\theta})d{x_{2}}$

where, $\phi(.)$ is the standard normal probability density function. In this
integration by parts and with the use of a change of variables as
$t=2\theta(z-x_{2})^{2}$, the second term can be expressed as:

$\displaystyle\int_{a}^{b}{x_{2}}\phi\bigg{(}\sqrt{2\theta}$
$\displaystyle(z-{x_{2}})\bigg{)}(-\sqrt{2\theta})d{x_{2}}$ $\displaystyle=$
$\displaystyle\int_{a}^{b}[z-(z-{x_{2}})](-\sqrt{2\theta})\frac{1}{\sqrt{2\pi}}\exp\bigg{\\{}-\frac{1}{2}\frac{(z-{x_{2}})^{2}}{(1/2\theta)}\bigg{\\}}d{x_{2}}$
$\displaystyle=$
$\displaystyle-z\int_{a}^{b}\frac{1}{\sqrt{2\pi(1/2\theta)}}\exp\bigg{\\{}-\frac{1}{2}\frac{(z-{x_{2}})^{2}}{(1/2\theta)}\bigg{\\}}d{x_{2}}+\int_{a}^{b}\frac{(z-{x_{2}})}{(1/\sqrt{2\theta})}\frac{1}{\sqrt{2\pi}}\exp\bigg{\\{}-\frac{1}{2}\frac{(z-{x_{2}})^{2}}{(1/2\theta)}\bigg{\\}}d{x_{2}}$
$\displaystyle=$
$\displaystyle-z\Phi\bigg{(}\sqrt{2\theta}({x_{2}}-z)\bigg{)}\bigg{|}_{a}^{b}+\int_{2\theta(z-a)^{2}}^{2\theta(z-b)^{2}}\sqrt{t}\frac{1}{\sqrt{2\pi}}\exp\bigg{\\{}-\frac{1}{2}t\bigg{\\}}\frac{dt}{2\sqrt{t}(-\sqrt{2\theta})}$
$\displaystyle=$
$\displaystyle-z\Phi\bigg{(}\sqrt{2\theta}({x_{2}}-z)\bigg{)}\bigg{|}_{a}^{b}-\frac{1}{4\sqrt{\pi\theta}}\int_{2\theta(z-a)^{2}}^{2\theta(z-b)^{2}}\exp\bigg{\\{}-\frac{1}{2}t\bigg{\\}}dt$
$\displaystyle=$
$\displaystyle-z\Phi\bigg{(}\sqrt{2\theta}({x_{2}}-z)\bigg{)}\bigg{|}_{a}^{b}+\frac{1}{2\sqrt{\pi\theta}}\exp\bigg{\\{}-\frac{1}{2}t\bigg{\\}}\bigg{|}_{2\theta(z-a)^{2}}^{2\theta(z-b)^{2}}$

Using this expression the integral in Eq. (A.6), can be evaluate as:

$\displaystyle\int_{a}^{b}\Phi$
$\displaystyle\bigg{(}\sqrt{2\theta}(b-{x_{2}})\bigg{)}d{x_{2}}-\int_{a}^{b}\Phi\bigg{(}\sqrt{2\theta}(a-{x_{2}})\bigg{)}d{x_{2}}$
$\displaystyle=a+b-2[a\Phi(\sqrt{2\theta}(b-a))+b\Phi(\sqrt{2\theta}(a-b))]-\frac{1}{\sqrt{\pi\theta}}[1-\exp\\{-\theta(b-a)^{2}\\}]$

Finally, we can compute the double integration of the correlation function in
closed form as:

$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}r(\theta,{x_{1}},{x_{2}})p_{{X_{1}}}({x_{1}})p_{{X_{2}}}({x_{2}})d{x_{1}}d{x_{2}}=\\\
\bigg{(}\frac{1}{b-a}\bigg{)}^{2}\sqrt{\frac{\pi}{\theta}}\bigg{[}a+b-2[a\Phi(\sqrt{2\theta}(b-a))+b\Phi(\sqrt{2\theta}(a-b))]-\frac{1}{\sqrt{\pi\theta}}[1-\exp\\{-\theta(b-a)^{2}\\}]\bigg{]}$
(A.7)

##  Appendix B Interaction Sensitivities

The interaction sensitivities are computed using a similar procedure as the
main effect sensitivities. These indices can be expressed using following
equation

$\tilde{S}_{\bm{p}}(\omega)=\dfrac{\text{Var}_{X_{\bm{p}}}(E_{X_{\sim\bm{p}}}[\mathcal{Y}(X,\omega)|X_{\bm{p}}])}{E[\text{Var}(\mathcal{Y}(X,\omega))]}$
(A.8)

where, $\bm{p}$ is a subset of indices of the input dimensions(i.e.
$\bm{p}\in\\{1,2,\ldots,d\\}$).

To compute the interaction Sobol indices, we start by taking the expectation
of the GP $\mathcal{Y}(X,\omega)$ over all the inputs except ‘$X_{\bm{p}}$’
and denote this as

$A(X_{\bm{p}},\omega)=E_{X_{\sim\bm{p}}}[\mathcal{Y}(X,\omega)|X_{\bm{p}}]$
(A.9)

Since, $\mathcal{Y}(X,\omega)$ is a Gaussian process and expectation is a
linear operator, $A(X_{\bm{p}},\omega)$ is also a GP referred to as the
interaction effect GP.

The mean and covariance function of the interaction effect GP can be
determined by integrating the original GP with respect to the joint
probability measure over all inputs except $X_{\bm{p}}$. Considering
independent inputs, the mean function is given by [49]:

$\displaystyle\mu_{A_{\bm{p}}}(x_{\bm{p}})=\mathbb{E}[A(X_{\bm{p}})]=\int_{\bm{x}_{\sim{\bm{p}}}}\hat{y}(\bm{x})\prod_{j\not\in{\bm{p}}}p_{X_{j}}(x_{j})dx_{j}$
(A.10)

which is easily computed as a product of one-dimensional integrals. Again
considering independent inputs, the covariance function of the interaction
effect GP is given by [49]:

$\displaystyle\text{Cov}(A(X_{1{\bm{p}}}),A(X_{2{\bm{p}}}))$
$\displaystyle=\int_{\bm{x}_{1\sim{\bm{p}}}}\int_{\bm{x}_{2\sim{\bm{p}}}}\text{Cov}(\mathcal{Y}(\bm{x}_{1}),\mathcal{Y}(\bm{x}_{2}))\prod_{j\not\in{\bm{p}}}p_{X_{j}}(x_{1j})d{x_{1j}}\prod_{j\not\in{\bm{p}}}p_{X_{j}}(x_{2j})d{x_{2j}}$
(A.11)

which can be expressed as a product of 2-dimensional integrals. After
obtaining the mean and covariance function of the interaction effect GP, the
estimates for sensitivities indices are obtained in a similar procedure as
main effect sensitivities. Note that the contributions of the interaction
sensitivities are not used during the learning process.