###### Abstract

Computer simulations have become essential for analyzing complex systems, but
high-fidelity simulations often come with significant computational costs. To
tackle this challenge, multi-fidelity computer experiments have emerged as a
promising approach that leverages both low-fidelity and high-fidelity
simulations, enhancing both the accuracy and efficiency of the analysis. In
this paper, we introduce a new and flexible statistical model, the Recursive
Non-Additive (RNA) emulator, that integrates the data from multi-fidelity
computer experiments. Unlike conventional multi-fidelity emulation approaches
that rely on an additive auto-regressive structure, the proposed RNA emulator
recursively captures the relationships between multi-fidelity data using
Gaussian process priors without making the additive assumption, allowing the
model to accommodate more complex data patterns. Importantly, we derive the
posterior predictive mean and variance of the emulator, which can be
efficiently computed in a closed-form manner, leading to significant
improvements in computational efficiency. Additionally, based on this
emulator, we introduce three active learning strategies that optimize the
balance between accuracy and simulation costs to guide the selection of the
fidelity level and input locations for the next simulation run. We demonstrate
the effectiveness of the proposed approach in a suite of synthetic examples
and a real-world problem. An R package for the proposed methodology is
provided in an open repository.

Active Learning for a Recursive Non-Additive Emulator for Multi-Fidelity
Computer Experiments
Junoh Heo and Chih-Li Sung***These authors gratefully acknowledge funding from
NSF DMS 2113407.
Michigan State University

Keywords: Surrogate model; Sequential design; Uncertainty quantification;
Gaussian process; Auto-regressive model.

##  1 Introduction

Computer simulations have become a crucial role in engineering and scientific
research, serving as valuable tools for predicting the performance of complex
systems. They have found successful applications in diverse fields, including
aerospace engineering (Mak et al.,, 2018), natural disaster prediction (Ma et
al.,, 2022), and cell biology (Sung et al.,, 2020). However, conducting high-
fidelity computer simulations for parameter space exploration can be quite
demanding due to the prohibitive costs and time-consuming computations
involved. To address this challenge, multi-fidelity emulation has emerged as a
promising alternative, which uses a statistical model to leverage
computationally expensive yet accurate high-fidelity simulations along with
computationally inexpensive but potentially less accurate low-fidelity
simulations, yielding an efficient predictive model emulating the expensive
computer code. By strategically integrating these simulations and carefully
designing the multi-fidelity experiments, we can potentially improve accuracy
without sacrificing excessive computational resources.

The usefulness of the multi-fidelity emulation framework has driven extensive
research in recent years. One popular approach is the Kennedy-O’Hagan (KO)
model (Kennedy and O’Hagan,, 2000), which models a sequence of computer
simulations from lowest to highest fidelity using a sequence of Gaussian
process (GP) models (Gramacy,, 2020; Rasmussen and Williams,, 2006), linked by
a linear auto-regressive framework. This model has made significant
contributions across various fields employing multi-fidelity computer
experiments (see, e.g., Patra et al., (2020), Kuya et al., (2011), and Demeyer
et al., (2017)), and several recent developments, including Qian et al.,
(2006), Le Gratiet, (2013), Le Gratiet and Garnier, (2014), Qian and Wu,
(2008), Perdikaris et al., (2017), and Ji et al., (2022) (among many others),
have investigated modeling strategies for efficient posterior prediction and
Bayesian uncertainty quantification.

Despite this body of work, most of these approaches rely on the assumption of
linear correlation between low-fidelity and high-fidelity data, resulting in
an additive GP structure. With the growing complexity of modern data, such
models face challenges in capturing complex relationships between data with
different fidelity levels. As shown in the left panel of Figure 1, where the
relationship between high-fidelity data and low-fidelity data is nonlinear,
the KO model falls short in providing accurate predictions due to its limited
flexibility resulting from the additive structure it assumes.

In this paper, we propose a new and flexible model that captures the nonlinear
relationships between multi-fidelity data in a recursive manner. This flexible
nonlinear functional form can encompass many existing models, including the KO
model, as a special case. Specifically, instead of assuming a specific form
for the relationship, our model uses GP priors to model the relationship
between multi-fidelity data, resulting in a non-additive GP model. Hence, we
refer to this proposed method as the Recursive Non-Additive (RNA) emulator. As
shown in the right panel of Figure 1, the RNA emulator demonstrates its
superior ability to the KO model by emulating the high-fidelity simulator with
high accuracy and low uncertainty, showcasing its flexibility and power of our
nonlinear modeling approach.

Figure 1: An example adapted from Perdikaris et al., (2017), where $n_{1}=13$
samples (red dots) are collected from the low-fidelity simulator
$f_{1}(x)=\sin(8\pi x)$ (red dashed line), and $n_{2}=8$ samples (black
triangles) are collected from the high-fidelity simulator
$f_{2}(x)=(x-\sqrt{2})f_{1}^{2}(x)$ (black solid line). The KO emulator (left
panel) and the RNA emulator (right panel) are shown as blue lines. Gray shaded
regions represent the 95% pointwise confidence intervals.

The RNA emulator belongs to the emerging field of linked/deep GP models (see,
e.g., Kyzyurova et al., (2018), Ming and Guillas, (2021), Sauer et al., (2023)
and Ming et al., (2023)), where different GPs are connected in a coupled
manner. To the best of our knowledge, there has been limited research on
extending such results for the analysis of multi-fidelity computer
experiments, and we aim to address this gap in our work. Notably, recent work
by Perdikaris et al., (2017) has made progress in this direction, but their
approach assumes an additive structure for the kernel function, and employs
the Monte Carlo integration to handle intractable posterior distributions. In
a similar vein, Meng and Karniadakis, (2020), Li et al., (2020), Meng et al.,
(2021), and Kerleguer et al., (2024) establish connections between different
fidelities using (Bayesian) neural networks. In contrast, our proposed model
not only provides great flexibility using GP priors with commonly used kernel
structures to connect multi-fidelity data, but also provides analytical
expressions for both the posterior mean and variance, enabling the expectation
propagation method to approximate the intractable posterior distribution. This
computational improvement allows for more efficient calculations, facilitating
efficient uncertainty quantification.

Leveraging this newly developed RNA emulator, we introduce three active
learning strategies to achieve enhanced accuracy while carefully managing the
limited simulation resources, particularly crucial for computationally
expensive simulations. Active learning, also known as sequential design,
involves sequentially searching for and acquiring new data points at optimal
locations based on a given sampling criterion, to construct an accurate
surrogate model/emulator. While active learning has been well-established for
single-fidelity GP emulators (Gramacy,, 2020; Rasmussen and Williams,, 2006;
Santner et al.,, 2018), research in the context of multi-fidelity computer
experiments is scarce and more challenging. This is because it requires
simultaneous selection of optimal input locations and fidelity levels,
accounting for their respective simulation costs. Although some recent works
have considered cost in specific cases like single-objective unconstrained
optimization (Huang et al.,, 2006; Swersky et al.,, 2013; He et al.,, 2017)
and global approximation (Xiong et al.,, 2013; Le Gratiet and Cannamela,,
2015; Stroh et al.,, 2022; Ehara and Guillas,, 2023; Sung et al.,, 2022), most
of these works are developed based on the KO model, and the active learning
for non-additive GP models has not been fully explored in the literature. In
addition, popular sampling criteria for global approximation, such as “Active
Learning MacKay” (McKay et al., (2000), ALM) and “Active Learning Cohn” (Cohn,
(1993), ALC), remain largely unexplored in the context of multi-fidelity
computer experiments. Recent successful applications of these sampling
criteria to other learning problems can be found in Binois et al., (2019),
Park et al., (2023), Koermer et al., (2023), and Sauer et al., (2023).

Our main contribution, in addition to the development of the RNA emulator,
lies in advancing active learning with these popular sampling criteria, based
on this newly developed RNA emulator. This advancement enables the optimal
selection of input locations and fidelity levels while considering the
associated costs, resulting in improved predictability for learning the
nonlinear relationships between multi-fidelity simulators. The closed-form
posterior mean and variance of the RNA emulator allow for efficient
computation when computing these sampling criteria, supporting the active
learning efforts. To facilitate broader usage and ensure reproducibility, we
implement an R (R Core Team,, 2018) package called RNAmf, which is available
in a public repository.

The structure of this article is as follows. In Section 2, we provide a brief
review of the KO model. Our proposed RNA emulator is introduced in Section 3.
Section 4 outlines our active learning strategies based on the RNA emulator.
Numerical and real data studies are presented in Sections 5 and 6,
respectively. Lastly, we conclude the paper in Section 7.

##  2 Background

In this section, we provide an overview of the background information
regarding multi-fidelity emulation.

###  2.1 Problem Setup

Let $f_{l}(\mathbf{x})$ represent the scalar simulation output of the computer
code with input parameters $\mathbf{x}\in\Omega\subseteq\mathbb{R}^{d}$ at
fidelity level $l=1,\ldots,L$. Throughout this paper, we assume that $L$
distinct fidelity levels of simulations are conducted for training an
emulator, where a higher fidelity level corresponds to a simulator with more
accurate outputs but also higher simulation costs per run.

Our primary objective is to construct an efficient emulator for the highest-
fidelity and most expensive simulation code, $f_{L}(\mathbf{x})$. For each
fidelity level $l$, we perform simulations at $n_{l}$ design points denoted by
$\mathcal{X}_{l}=\\{\mathbf{x}_{i}^{[l]}\\}_{i=1}^{n_{l}}$. These simulations
yield the corresponding simulation outputs
$\mathbf{y}_{l}:=(f_{l}(\mathbf{x}))_{\mathbf{x}\in\mathcal{X}_{l}}$,
representing the vector of outputs for $f_{l}(\mathbf{x})$ at design points
$\mathbf{x}\in\mathcal{X}_{l}$, and each element of $\mathbf{y}_{l}$ is
denoted by $y^{[l]}_{i}=f_{l}(\mathbf{x}^{[l]}_{i})$. We assume that the
designs $\mathcal{X}_{l}$ are sequentially nested, i.e.,

$\mathcal{X}_{L}\subseteq\mathcal{X}_{L-1}\subseteq\cdots\subseteq\mathcal{X}_{1}\subseteq\Omega,$
(1)

and $\mathbf{x}^{[l]}_{i}=\mathbf{x}^{[l-1]}_{i}$ for $i=1,\ldots,n_{l}$. In
other words, design points run for a higher fidelity simulator are contained
within the design points run for a lower fidelity simulator. This property has
been shown to lead to more efficient inference in various multi-fidelity
emulation approaches (Qian,, 2009; Qian et al.,, 2009; Haaland and Qian,,
2010).

Furthermore, we let $C_{l}$ denote the simulation cost (e.g., in CPU hours)
for a single run of the simulator at fidelity level $l$. Since higher-fidelity
simulators are more computationally demanding, this implies that
$0<C_{1}<C_{2}<\ldots<C_{L}$.

###  2.2 Auto-regressive model

One of the prominent approaches for modeling $f_{L}(\mathbf{x})$ is the auto-
regressive model proposed by Kennedy and O’Hagan, (2000), referred to as the
KO model hereafter, which can be expressed as follows:

$\displaystyle\begin{cases}&f_{1}(\mathbf{x})=Z_{1}(\mathbf{x}),\\\
&f_{l}(\mathbf{x})=\rho_{l-1}f_{l-1}(\mathbf{x})+Z_{l}(\mathbf{x}),\quad\text{for}\quad
2\leq l\leq L,\end{cases}$ (2)

where $\rho_{l-1}$ is an unknown auto-regressive parameter, and
$Z_{l}=(f_{l}-\rho_{l-1}f_{l-1})$ represents the discrepancy between the
$(l-1)$-th and $l$-th code. The KO model considers a probabilistic surrogate
model by assuming that $\\{Z_{l}\\}^{L}_{l=1}$ follow independent zero-mean
Gaussian process (GP) models:

$Z_{l}(\mathbf{x})\overset{indep.}{\sim}\mathcal{GP}\\{\alpha_{l}(\mathbf{x}),\tau^{2}_{l}\Phi_{l}(\mathbf{x},\mathbf{x}^{\prime})\\},\quad
l=1,\ldots,L,$ (3)

where $\alpha_{l}(\mathbf{x})$ is a mean function, $\tau^{2}_{l}$ is a
variance parameter, and $\Phi_{l}(\mathbf{x},\mathbf{x}^{\prime})$ is a
positive-definite kernel function defined on $\Omega\times\Omega$. In the
original KO paper, $\alpha_{l}(\mathbf{x})$ is assumed to be
$h(\mathbf{x})\beta_{l}$, where $h(\mathbf{x})$ is a vector of $d$ regression
functions. Other common choices for $\alpha_{l}(\mathbf{x})$ include
$\alpha_{l}(\mathbf{x})\equiv 0$ or $\alpha_{l}(\mathbf{x})=\mu$. As for the
kernel function, popular choices include the squared exponential kernel and
Matérn kernel (Stein,, 1999). Specifically, the anisotropic squared
exponential kernel takes the form:

$\displaystyle\Phi_{l}(\mathbf{x},\mathbf{x}^{\prime})=\prod^{d}_{j=1}\phi(x_{j},x^{\prime}_{j};\theta_{lj})=\prod^{d}_{j=1}\exp\left(-\frac{\left(x_{j}-x_{j}^{\prime}\right)^{2}}{\theta_{lj}}\right),$
(4)

where $(\theta_{l1},\ldots,\theta_{ld})$ is the lengthscale hyperparameter,
indicating that the correlation decays exponentially fast in the squared
distance between $\mathbf{x}$ and $\mathbf{x}^{\prime}$. The GP model (3),
combined with the auto-regressive model (2), implies that, conditional on the
parameters $\tau_{l}^{2},\theta_{lj}$, and $\mu(\cdot)$, the joint
distribution of $(\mathbf{y}_{1},\ldots,\mathbf{y}_{L})$ follows a
multivariate normal distribution, and these unknown parameters can be
estimated via maximum likelihood estimation or Bayesian inference. Given the
data $(\mathbf{y}_{1},\ldots,\mathbf{y}_{L})$, it can be shown that the
posterior distribution of $f_{L}(\mathbf{x})$ is also a GP. The posterior mean
function can then be used to emulate the expensive simulator, while the
posterior variance function can be employed to quantify the emulation
uncertainty. We refer the details to Kennedy and O’Hagan, (2000).

The auto-regressive framework of the KO model has led to the development of
several variants. For instance, Le Gratiet and Garnier, (2014) introduce a
faster algorithm based on a recursive formulation for computing the posterior
of $f_{l}(\mathbf{x})$ more efficiently. To enhance the model’s flexibility,
Qian et al., (2006), Le Gratiet and Garnier, (2014), and Qian and Wu, (2008)
allow the auto-regressive parameter $\rho_{l}$ to depend on the input
$\mathbf{x}$, that is,
$f_{l}(\mathbf{x})=\rho_{l-1}(\mathbf{x})f_{l-1}(\mathbf{x})+Z_{l}(\mathbf{x})$
for $2\leq l\leq L$, where the first two assume $\rho_{l-1}(\mathbf{x})$ to be
a linear function, while the last one assume $\rho_{l-1}(\mathbf{x})$ to be a
GP.

##  3 Recursive Non-Additive (RNA) emulator

Despite the advantages of the KO model, it results in an additive GP model
based on (2) and (3), which may not adequately capture the nonlinear
relationships between data at different fidelity levels. To overcome this
limitation and achieve a more flexible representation, we propose a novel
Recursive Non-Additive (RNA) emulator, defined in a recursive fashion:

$\displaystyle\begin{cases}&f_{1}(\mathbf{x})=W_{1}(\mathbf{x}),\\\
&f_{l}(\mathbf{x})=W_{l}(\mathbf{x},f_{l-1}(\mathbf{x})),\quad\text{for}\quad
l=2,\ldots,L.\end{cases}$ (5)

The model structure is illustrated in Figure 2. This RNA model offers greater
flexibility and can be viewed as many existing models as a special case. For
instance, the auto-regressive KO model can be represented in the form of
(LABEL:eq:RNA) by setting
$W_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))=\rho_{l-1}f_{l-1}(\mathbf{x})+\delta(\mathbf{x})$,
and it can also be considered as the model in Qian et al., (2006), Le Gratiet
and Garnier, (2014), and Qian and Wu, (2008) by setting
$W_{l}(x,f_{l-1}(\mathbf{x}))=\rho_{l-1}(\mathbf{x})f_{l-1}(\mathbf{x})+\delta(\mathbf{x})$.

$\mathcal{X}_{1}$$W_{1}$$\mathbf{y}_{1}$$\mathcal{X}_{2}$$W_{2}$$\mathbf{y}_{2}$$\mathcal{X}_{3}$$W_{3}$$\mathbf{y}_{3}$⋮⋮$\mathcal{X}_{L}$$W_{L}$$\mathbf{y}_{L}$

Figure 2: An illustration of the recursive structure of the RNA model.

Instead of assuming a specific form for the relationship, we model the
relationship $W_{l}$ using a GP prior, that is,
$W_{1}(\mathbf{x})\sim\mathcal{GP}\\{\alpha_{1}(\mathbf{x}),\tau^{2}_{l}\Phi_{1}(\mathbf{x},\mathbf{x}^{\prime})\\}$
and

$W_{l}(\mathbf{z})\overset{indep.}{\sim}\mathcal{GP}\\{\alpha_{l}(\mathbf{z}),\tau^{2}_{l}K_{l}(\mathbf{z},\mathbf{z}^{\prime})\\},\quad
l=2,\cdots,L,$ (6)

where $\mathbf{z}=(\mathbf{x},y)$, forming a vector of size $(d+1)$, and
$K_{l}(\mathbf{z},\mathbf{z}^{\prime})$ is a positive-definite kernel. We
consider a constant mean, i.e., $\alpha_{1}(\mathbf{x})=\alpha_{1}$ and
$\alpha_{l}(\mathbf{z})=\alpha_{l}$ for $l\geq 2$. We adopt popular kernel
choices for $K_{l}$, such as the squared exponential kernel and Matérn kernel.
In particular, the squared exponential kernel $K_{l}$, following (4), can be
expressed as:

$\displaystyle
K_{l}(\mathbf{z},\mathbf{z}^{\prime})=\phi(y,y^{\prime};\theta_{ly})\prod^{d}_{j=1}\phi(x_{j},x^{\prime}_{j};\theta_{lj})=\exp\left(-\frac{\left(y-y^{\prime}\right)^{2}}{\theta_{ly}}\right)\prod^{d}_{j=1}\exp\left(-\frac{\left(x_{j}-x_{j}^{\prime}\right)^{2}}{\theta_{lj}}\right),$
(7)

where the lengthscale hyperparameter,
$\bm{\theta}_{l}=(\theta_{l1},\ldots,\theta_{ld},\theta_{ly})$, represents a
vector of size $(d+1)$. The Matérn kernel can be similarly constructed, which
is given in the Appendix A.

Combining the GP model (6) with the recursive formulation (LABEL:eq:RNA) and
assuming the nested design as defined in (1), the observed simulations
$\mathbf{y}_{l}$ follow a multivariate normal distribution:

$\mathbf{y}_{l}\overset{indep.}{\sim}\mathcal{N}_{n_{l}}(\alpha_{l}\mathbf{1}_{n_{l}},\tau^{2}_{l}\mathbf{K}_{l})\quad\text{for}\quad
l=1,\ldots,L,$

where $\mathbf{1}_{n_{l}}$ is a unit vector of size $n_{l}$, $\mathbf{K}_{l}$
is an $n_{l}\times n_{l}$ matrix with each element
$(\mathbf{K}_{l})_{ij}=\Phi_{1}(\mathbf{x}^{[1]}_{i},\mathbf{x}^{[1]}_{j})$
for $l=1$, and
$(\mathbf{K}_{l})_{ij}=K_{l}(\mathbf{z}^{[l]}_{i},\mathbf{z}^{[l]}_{j})$ for
$l\geq 2$, where
$\mathbf{z}^{[l]}_{i}=(\mathbf{x}^{[l]}_{i},f_{l-1}(\mathbf{x}^{[l]}_{i}))$.
Note that $f_{l-1}(\mathbf{x}^{[l]}_{i})=y^{[l-1]}_{i}$ due to the nested
design assumption (1), which is the $i$-th simulation output at level $l-1$.
The parameters $\\{\alpha_{l},\tau^{2}_{l},\bm{\theta}_{l}\\}^{L}_{l=1}$ can
be estimated by maximum likelihood estimation. Specifically, the log-
likelihood (up to an additive constant) is

$-\frac{1}{2}\left(n_{l}\log(\tau_{l}^{2})+\log(\det(\mathbf{K}_{l}))+\frac{1}{\tau^{2}_{l}}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}})^{T}\mathbf{K}_{l}^{-1}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}})\right).$

The parameters can then be efficiently estimated by maximizing the log-
likelihood via an optimization algorithm, such as quasi-Newton optimization
method of Byrd et al., (1995).

It is important to acknowledge that the idea of a recursive GP was previously
proposed by Perdikaris et al., (2017), referred to as a nonlinear auto-
regressive model therein. However, there are two key distinctions between
their model and ours. The first distinction lies in the kernel assumption,
where they assume an additive form of the kernel to better capture the auto-
regressive nature. Specifically, they use
$K_{l}(\mathbf{z},\mathbf{z}^{\prime})=\Phi_{l1}(\mathbf{x},\mathbf{x}^{\prime})\Phi_{l2}(f_{l-1}(\mathbf{x}),f_{l-1}(\mathbf{x}^{\prime}))+\Phi_{l3}(\mathbf{x},\mathbf{x}^{\prime})$
with valid kernel functions $\Phi_{l1}$, $\Phi_{l2}$, and $\Phi_{l3}$. While
our kernel function shares some similarities, particularly the first component
$\Phi_{l1}(\mathbf{x},\mathbf{x}^{\prime})\Phi_{l2}(f_{l-1}(\mathbf{x}),f_{l-1}(\mathbf{x}^{\prime}))$,
the role of the second component $\Phi_{l3}$ in predictions remains unclear.
The inclusion of this component introduces $d$ hyperparameters for an
anisotropic kernel, making the estimation more challenging with limited
samples, especially for high-dimensional problems. In contrast, we adopt the
natural form of popular kernel choices, such as the squared exponential kernel
in (7) and the Matérn kernel, placing our model within the emerging field of
linked/deep GP models, which has shown promising results in the computer
experiment literature (Kyzyurova et al.,, 2018; Ming and Guillas,, 2021; Sauer
et al.,, 2023). The second distinction is in the computation for the posterior
of $f_{L}(\mathbf{x})$. Specifically, their model relies on Monte Carlo (MC)
integration for their computation, which can be quite computationally
demanding, especially in this recursive formulation. In contrast, with these
popular kernel choices, we can derive the posterior mean and variance of
$f_{L}(\mathbf{x})$ in a closed form, which is presented in following
proposition, enabling more efficient predictions and uncertainty
quantification.

The derivation of the posterior follows these steps. Starting with the GP
assumption (6), and utilizing the properties of conditional multivariate
normal distribution, the posterior distribution of $f_{l}$ given
$\mathbf{y}_{l}$ and $f_{l-1}$ at a new input location $\mathbf{x}$ is
normally distributed, namely:
$f_{1}(\mathbf{x})|\mathbf{y}_{1}\sim\mathcal{N}(\mu_{1}(\mathbf{x}),\sigma_{1}^{2}(\mathbf{x}))$
with

$\displaystyle\mu_{1}(\mathbf{x})$
$\displaystyle=\alpha_{1}\mathbf{1}_{n_{1}}+\mathbf{k}_{1}^{T}(\mathbf{x})\mathbf{K}^{-1}_{1}(\mathbf{y}_{1}-\alpha_{l}\mathbf{1}_{n_{1}}),\quad\text{and}$
(8)
$\displaystyle\sigma^{2}_{1}(\mathbf{x})=\tau^{2}_{1}(1-\mathbf{k}_{1}(\mathbf{x})^{T}\mathbf{K}^{-1}_{1}\mathbf{k}_{1}(\mathbf{x})),$
(9)

and
$f_{l}(\mathbf{x})|\mathbf{y}_{l},f_{l-1}(\mathbf{x})\sim\mathcal{N}(\mu_{l}(\mathbf{x},f_{l-1}(\mathbf{x})),\sigma_{l}^{2}(\mathbf{x},f_{l-1}(\mathbf{x})))$
for $l=2,\ldots,L$ with

$\displaystyle\mu_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))$
$\displaystyle=\alpha_{l}\mathbf{1}_{n_{l}}+\mathbf{k}_{l}^{T}(\mathbf{x},f_{l-1}(\mathbf{x}))\mathbf{K}^{-1}_{l}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}}),\quad\text{and}$
(10) $\displaystyle\sigma^{2}_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))$
$\displaystyle=\tau^{2}_{l}(1-\mathbf{k}_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))^{T}\mathbf{K}^{-1}_{l}\mathbf{k}_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))),$
(11)

where $\mathbf{k}_{1}(\mathbf{x})$ and
$\mathbf{k}_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))$ are an $n_{l}\times 1$ matrix
with each element
$(\mathbf{k}_{1}(\mathbf{x}))_{i,1}=\Phi_{l}(\mathbf{x},\mathbf{x}^{[l]}_{i})$
and
$(\mathbf{k}_{l}(\mathbf{x},f_{l-1}(\mathbf{x})))_{i,1}=K_{l}((\mathbf{x},f_{l-1}(\mathbf{x})),(\mathbf{x}_{i}^{[l]},y^{[l-1]}_{i}))$
for $l\geq 2$, respectively. The posterior distribution of $f_{l}$ can then be
obtained by

$\displaystyle p(f_{l}(\mathbf{x})|\mathbf{y}_{1},\ldots,\mathbf{y}_{l})$
$\displaystyle=$ $\displaystyle\int\cdots\int
p(f_{l}(\mathbf{x})|\mathbf{y}_{l},f_{l-1}(\mathbf{x}))p(f_{l-1}(\mathbf{x})|\mathbf{y}_{l-1},f_{l-2}(\mathbf{x}))\cdots
p(f_{1}(\mathbf{x})|\mathbf{y}_{1})\text{d}(f_{l-1}(\mathbf{x}))\ldots\text{d}(f_{1}(\mathbf{x})).$

This posterior is analytically intractable but can be numerically approximated
using MC integration, as done in Perdikaris et al., (2017), which involves
sequential sampling from the normal distribution
$p(f_{l}(\mathbf{x})|\mathbf{y}_{l},f_{l-1}(\mathbf{x}))$ from $l=1$ to $l=L$.
However, this method can be computationally demanding and can become
impractical for uncertainty quantification, especially when the dimensions of
$\mathbf{x}$ and the number of fidelity levels increase. To address this, we
derive recursive closed-form expressions for the posterior mean and variance
under popular kernel choices as follows.

######  Proposition 3.1.

Under the squared exponential kernel function (7), the posterior mean and
variance of $f_{l}(\mathbf{x})$ given the data
$\mathbf{y}_{1},\ldots,\mathbf{y}_{l}$ for $l\geq 2$ can be expressed in a
recursive fashion as follows,

$\displaystyle\mu^{*}_{l}(\mathbf{x}):$
$\displaystyle=\mathbb{E}[f_{l}(\mathbf{x})|\mathbf{y}_{1},\ldots,\mathbf{y}_{l}]$
$\displaystyle=\alpha_{l}+\sum^{n_{l}}_{i=1}r_{i}\prod_{j=1}^{d}\exp\left(-\frac{(x_{j}-x^{[l]}_{ij})^{2}}{\theta_{lj}}\right)\frac{1}{\sqrt{1+2\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}}}\exp{\left(-\frac{(y^{[l-1]}_{i}-\mu^{*}_{l-1}(\mathbf{x}))^{2}}{\theta_{ly}+2\sigma^{*2}_{l-1}(\mathbf{x})}\right)},$

and

$\displaystyle\sigma^{*2}_{l}(\mathbf{x}):=\mathbb{V}$
$\displaystyle[f_{l}(x)|\mathbf{y}_{1},\ldots,\mathbf{y}_{l}]=\tau^{2}_{l}-(\mu^{*}_{l}(\mathbf{x})-\alpha_{l})^{2}+$
$\displaystyle\left(\sum_{i,k=1}^{n_{l}}\zeta_{ik}\left(r_{i}r_{k}-\tau^{2}_{l}(\mathbf{K}^{-1}_{l})_{ik}\right)\prod_{j=1}^{d}\exp{\left(-\frac{(x_{j}-x^{[l]}_{ij})^{2}+(x_{j}-x^{[l]}_{kj})^{2}}{\theta_{lj}}\right)}\right),$
(12)

where
$r_{i}=(\mathbf{K}^{-1}_{l}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}}))_{i}$,
and

$\displaystyle\zeta_{ik}=\frac{1}{\sqrt{1+4\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}}}\exp{\left(-\frac{(\frac{y^{[l-1]}_{i}+y^{[l-1]}_{k}}{2}-\mu^{*}_{l-1}(\mathbf{x}))^{2}}{\frac{\theta_{ly}}{2}+2\sigma^{*2}_{l-1}(\mathbf{x})}-\frac{(y^{[l-1]}_{i}+y^{[l-1]}_{k})^{2}}{2\theta_{ly}}\right)}.$
(13)

For $l=1$, it follows that $\mu^{*}_{1}(\mathbf{x})=\mu_{1}(\mathbf{x})$ and
$\sigma^{*2}_{1}(\mathbf{x})=\sigma^{2}_{1}(\mathbf{x})$ as in (8) and (9),
respectively.

The posterior mean and variance under a Matérn kernel with the smoothness
parameter $\nu=1.5$ and $\nu=2.5$ are provided in the Supplementary Materials
S2, and the detailed derivations for Proposition 3.1 are provided in
Supplementary Materials S1, which follow the proof of Kyzyurova et al., (2018)
and Ming and Guillas, (2021). With this proposition, the posterior mean and
variance can be efficiently computed in a recursive fashion. Similar to
Kyzyurova et al., (2018) and Ming and Guillas, (2021), we adopt the moment
matching method, using a Gaussian distribution to approximate the posterior
distribution with the mean and variance presented in the proposition. The
parameters in the posterior distribution, including
$\\{\alpha_{l},\tau^{2}_{l},\bm{\theta}_{l}\\}^{L}_{l=1}$, can be plugged in
by their estimates.

Similar to conventional GP emulators for single-fidelity deterministic
computer models, the proposed RNA emulator also exhibits the interpolation
property, which is described in the following proposition. The proof is
provided in Appendix B.

######  Proposition 3.2.

The RNA emulator satisfies interpolation property, that is,
$\mu^{*}_{l}(\mathbf{x}_{i}^{[l]})=y_{i}^{[l]}$, and
$\sigma^{*2}_{l}(\mathbf{x}_{i}^{[l]})=0$, where
$\\{(\mathbf{x}_{i}^{[l]},y_{i}^{[l]})\\}_{i=1,\ldots,n_{l}}$ are the training
samples.

An example of this posterior distribution is presented in the right panel of
Figure 1, illustrating that the posterior mean closely aligns with the true
function, and the confidence intervals contructed by the posterior variance
cover the true function. We have developed an R package called RNAmf, which
implements the parameter estimation and computations for the closed-form
posterior mean and variance using a Gaussian kernel and a Matérn kernel with
smoothness parameters of $1.5$ and $2.5$. The closed-form posterior mean and
variance are particularly useful for efficient computation of the criteria
employed in active learning, which is introduced in the next section.

##  4 Active learning for RNA emulator

In this section, we present three active learning (AL) strategies aimed at
enhancing the predictive capabilities of the proposed model through careful
design of computer experiments. Different from traditional AL approaches,
these strategies encompass the dual task of not only identifying the optimal
input locations but also determining the most appropriate fidelity level.
These strategies are inspired from popular sampling criteria in AL for GPs:
“Active Learning MacKay (ALM)” and “Active Learning Cohn (ALC)” (Seo et al.,,
2000).

We suppose that an initial experiment of sample size $n_{l}$ for each fidelity
level $l$, following a nested design
$\mathcal{X}_{L}\subseteq\cdots\subseteq\mathcal{X}_{1}$, is conducted, for
which a space-filling design is often considered, such as the nested Latin
hypercube design (Qian,, 2009). AL seeks to optimize a selection criterion for
choosing the next point $\mathbf{x}^{[l]}_{n_{l}+1}$ at fidelity level $l$,
carrying out its corresponding simulation
$y^{[l]}_{n_{l}+1}=f_{l}(\mathbf{x}^{[l]}_{n_{l}+1})$, and thus augmenting the
dataset.

###  4.1 Active Learning MacKay (ALM)

A straightforward but commonly used sampling criterion in AL is to select the
next point that maximizes the posterior predictive variance (MacKay,, 1992).
Extending this concept to our scenario by accounting for the simulation cost
$C_{l}$, we choose the next point $\mathbf{x}^{[l]}_{n_{l}+1}$ at fidelity
level $l$ by maximizing the ALM citerion:

$\displaystyle(l^{*},\mathbf{x}_{n_{l^{*}}+1}^{[l^{*}]})=\operatorname*{argmax}_{l\in\\{1,\ldots,L\\};\mathbf{x}\in\Omega}\frac{\sigma^{*2}_{l}(\mathbf{x})}{\sum^{l}_{j=1}C_{j}},$
(14)

which aims to maximize the ratio between the posterior predictive variance
$\sigma^{*2}_{l}(\mathbf{x})$ and the simulation cost $\sum^{l}_{j=1}C_{j}$ at
each fidelity level $l$. Note that after running the simulation at the optimal
input location $\mathbf{x}^{[l^{*}]}_{n_{l^{*}}+1}$ at level $l^{*}$, the
posterior predictive variance
$\sigma^{*2}_{l}(\mathbf{x}^{[l^{*}]}_{n_{l^{*}}+1})$ becomes zero (see
Proposition 3.2). In other words, our selection of the optimal level hinges on
achieving the highest ratio of uncertainty reduction at
$\mathbf{x}^{[l^{*}]}_{n_{l^{*}}+1}$ to the simulation cost.

The reason behind assigning a simulation cost of $\sum^{l}_{j=1}C_{j}$ to the
selection of $\mathbf{x}_{n_{l}+1}^{[l]}$ is because of the nested structure
assumption as in (1) for fitting the RNA emulator, implying that we also need
to conduct simulations with the same optimal input location at lower fidelity
levels. That is, to run the simulation
$f_{l^{*}}(\mathbf{x}_{n_{l^{*}}+1}^{[l^{*}]})$, we also need to run
$f_{l}(\mathbf{x}_{n_{l}+1}^{[l]})$ with
$\mathbf{x}^{[l]}_{n_{l}+1}=\mathbf{x}^{[l^{*}]}_{n_{l^{*}}+1}$ for all $1\leq
l<l^{*}$. It is also worth mentioning that the cost can be tailored to depend
on the input $\mathbf{x}$, as done in He et al., (2017) and Stroh et al.,
(2022).

The computation of this ALM criterion is facilitated by the availability of
the closed-form expression of the posterior predictive variance as in (3.1),
which in turn simplifies the optimization process of (14). In particular, the
optimal input location $\mathbf{x}^{[l]}_{n_{l}+1}$ for each $l$ can be
efficiently obtained through the optim library in R, using the method=L-BFGS-B
option, which performs a quasi-Newton optimization approach of Byrd et al.,
(1995).

###  4.2 Active Learning Cohn (ALC)

Another widely employed, more aggregate criterion is Active Learning Cohn
(ALC) (Cohn,, 1993; Seo et al.,, 2000). In contrast to ALM, ALC selects an
input location that maximizes the reduction in posterior variances across the
entire input space after running this selected simulation. Extending the
concept to our scenario by accounting for the simulation cost $C_{l}$, we
choose the next point $\mathbf{x}_{n_{l}+1}$ at fidelity level $l$ by
maximizing the ALC criterion:

$\displaystyle(l^{*},\mathbf{x}_{n_{l^{*}}+1}^{[l^{*}]})=\operatorname*{argmax}_{l\in\\{1,\ldots,L\\};\mathbf{x}\in\Omega}\frac{\Delta\sigma_{L}^{2}(l,\mathbf{x})}{\sum^{l}_{j=1}C_{j}},$
(15)

where $\Delta\sigma_{L}^{2}(l,\mathbf{x})$ is the average reduction in
variance (of the highest-fidelity emulator) from the current design measured
through a choice of the fidelity level $l$ and the input location
$\mathbf{x}$, augmenting the design, that is,

$\displaystyle\Delta\sigma_{L}^{2}(l,\mathbf{x})=\int_{\Omega}\sigma_{L}^{*2}(\bm{\xi})-\tilde{\sigma}_{L}^{*2}(\bm{\xi};l,\mathbf{x}){\rm{d}}\bm{\xi},$
(16)

where $\sigma_{L}^{*2}(\bm{\xi})$ is the posterior variance of
$f_{L}(\bm{\xi})$ based on the current design
$\\{\mathcal{X}_{l}\\}^{L}_{l=1}$, and
$\tilde{\sigma}_{L}^{*2}(\bm{\xi};l,\mathbf{x})$ is the posterior variance
based on the augmented design combining the current design and a new input
location $\mathbf{x}$ at each fidelity level lower than or equal to $l$, i.e.,
$\\{(\mathcal{X}_{1}\cup\mathbf{x}_{n_{1}+1}^{[1]}),\ldots,(\mathcal{X}_{l}\cup\mathbf{x}_{n_{l}+1}^{[l]}),\mathcal{X}_{l+1},\ldots,\mathcal{X}_{L}\\}$
with
$\mathbf{x}_{n_{1}+1}^{[1]}=\cdots=\mathbf{x}_{n_{l}+1}^{[l]}=\mathbf{x}$.
Once again, the incorporation of the new input location $\mathbf{x}$ at each
fidelity level lower than $l$ is due to the nested structure assumption. In
other words, our selection of the optimal level involves maximizing the ratio
of average reduction in variance of the highest-fidelity emulator to the
associated simulation cost. In practice, the integration in (16) can be
approximated by numerical methods, such as MC integration.

It is worth noting that, unlike ALM where the influence of design augmentation
on the variance of the highest-fidelity emulator is unclear, ALC is
specifically designed to maximize the reduction in variance of the highest-
fidelity emulator. However, the ALC strategy involves requiring knowledge of
future outputs $y^{[s]}_{n_{s}+1}=f_{s}(\mathbf{x}^{[s]}_{n_{s}+1})$ for all
$1\leq s\leq l$, as they are involved in
$\tilde{\sigma}_{L}^{*2}(\bm{\xi};l,\mathbf{x})$ (as seen in (13)). These
outputs are not available prior to conducting the simulations, rendering the
posterior variance inaccessible. A possible approach to address this issue is
through MC approximation to impute the outputs. Specifically, we can impute
$y^{[s]}_{n_{s}+1}$ for each $1\leq s\leq l$ by drawing samples from the
posterior distribution of $f_{s}(\mathbf{x}^{[s]}_{n_{s}+1})$ based on the
current design, which is a normal distribution with the posterior mean and
variance presented in Proposition 3.1. This allows us to repeatedly compute
$\tilde{\sigma}_{L}^{*2}(\bm{\xi};l,\mathbf{x})$ using the imputations and
average the results to approximate the variance. Notably, with the imputed
output $y^{[s]}_{n_{s}+1}$, the variance
$\tilde{\sigma}_{L}^{*2}(\bm{\xi};l,\mathbf{x})$ can be efficiently computed
using the Sherman–Morrison formula (Harville,, 1998) for updating the
covariance matrix’s inverse, $\mathbf{K}_{l}^{-1}$, from
$\sigma_{L}^{*2}(\bm{\xi})$. The details regarding this efficient inverse
covariance matrix update can be found in Chapter 6.3 of Gramacy, (2020).

In contrast to ALM, searching for the optimal $l$ and $\mathbf{x}$ that
maximize the ALC criterion (15) can be quite computationally expensive due to
the costly Monte Carlo approximation to compute (16). To this end, an
alternative strategy is proposed to strike a compromise by combining the two
criteria.

###  4.3 Two-step approach: ALMC

Given the distinct advantages and limitations of both ALM and ALC criteria
(details of which are referred to Chapter 6 of Gramacy, (2020)), for a
comprehensive exploration, we can contemplate their combination. Inspired by
Le Gratiet and Cannamela, (2015), we introduce a hybrid approach, which we
refer to as ALMC. First, the optimal input location is selected by maximizing
the posterior predictive variance of the highest fidelity emulator:

$\displaystyle\mathbf{x}^{*}=\operatorname*{argmax}_{\mathbf{x}\in\Omega}\sigma^{*2}_{L}(\mathbf{x}).$

Then, the ALC criterion determines the fidelity level with the identified
input location:

$\displaystyle
l^{*}=\operatorname*{argmax}_{l\in\\{1,\ldots,L\\}}\frac{\Delta\sigma_{L}^{2}(l,\mathbf{x}^{*})}{\sum^{l}_{j=1}C_{j}},$

which aims to maximize the ratio between the variance reduction and the
associated simulation cost. This synergistic approach is not only expected to
capture the advantages of both ALM and ALC, but also offers computational
efficiency advantages compared to the ALC method in the previous subsection.
This is due to the fact that the optimization for $\mathbf{x}^{*}$ by
maximizing the closed-form posterior variance is computationally much cheaper,
as discussed in Section 4.1.

Figure 3: Demonstration of the three active learning strategies for selecting
the next data point to augment the dataset, using the example in the right
panel of Figure 1. In the bottom panel, we visualize the criteria of the three
strategies: ALM (left), ALC (middle), and ALMC (right), plotted across a dense
grid of input locations, where the dots represent the optimal input locations
for each of the simulators. Notably, ALMC employs the gray line to determine
the optimal input location and then utilizes the red and black lines (which
are identical to ALC) to decide the fidelity level. The upper panels show the
corresponding fits after adding the selected points to the training dataset,
where the solid dots represent the chosen samples, all of which select the
low-fidelity simulator.

Figure 3 demonstrates the effectiveness of these three strategies, showing
that how they are employed to determine the next input location with a
suitable fidelity level to run the simulation and augment the dataset,
particularly for the example in right panel of Figure 1 where the RNA emulator
is employed. Consider the simulation costs: $C_{1}=1$ and $C_{2}=3$ for the
two simulators. It shows that, for all the three criteria, the choice is
consistently in favor of selecting the low-fidelity simulator to augment the
dataset. While the selected locations differ, they all fall within the range
of $[0.18,0.25]$, which, as per the current design (prior to running this
simulation), holds the highest uncertainty, as seen in the right panel of
Figure 1. All these selection outcomes contribute to an overall improvement in
emulation accuracy, while simultaneously reducing global uncertainty, even
when opting for low-fidelity data alone, highlighting the effectiveness of the
three strategies.

##  5 Numerical Studies

In this section, we conduct a suite of numerical experiments to examine the
performance of the proposed approach. The experiments encompass two main
aspects. In Section 5.1, we assess the predictive capabilities of the proposed
RNA emulator, while Section 5.2 delves into the evaluation of the performance
of the proposed AL strategies.

We consider the anisotropic squared exponential kernel as in (7) for the
proposed model, a choice that is also shared by our competing methods. All
experiments are performed on a MacBook Pro laptop with 2.9 GHz 6-Core Intel
Core i9 and 16Gb of RAM.

###  5.1 Emulation performance

We begin by comparing the predictive performance of the proposed RNA emulator
(labeled RNAmf) with two other methods in the numerical experiments: the co-
kriging model (labeled CoKriging) by Le Gratiet and Garnier, (2014), and the
nonlinear auto-regressive multi-fidelity GP (labeled NARGP) by Perdikaris et
al., (2017). The two methods are readily available through open repositories,
specifically the R package MuFiCokriging (Le Gratiet,, 2012) and the Python
package on the GitHub repository (Perdikaris,, 2016), respectively.

Six synthetic examples are considered, including the two-level Perdikaris
function (Perdikaris et al.,, 2017),

$\displaystyle\begin{cases}&f_{1}(x)=\sin(8\pi x)\\\
&f_{2}(x)=(x-\sqrt{2})f_{1}^{2}(x)\end{cases}\text{for }x\in[0,1],$

the three-level Branin function (Sobester et al.,, 2008),

$\displaystyle\begin{cases}&f_{1}(\mathbf{x})=f_{2}(1.2(\mathbf{x}+2))-3x_{2}+1\\\
&f_{2}(\mathbf{x})=10\sqrt{f_{3}(\mathbf{x})}+2(x_{1}-0.5)-3(3x_{2}-1)-1\\\
&f_{3}(\mathbf{x})=\left(\frac{-1.275x_{1}^{2}}{\pi^{2}}+\frac{5x_{1}}{\pi}+x_{2}-6\right)^{2}+\left(10-\frac{5}{4\pi}\right)\cos(x_{1})+10\end{cases}\text{for
}\mathbf{x}\in[-5,10]\times[0,15],$

the two-level Park function (Park,, 1991; Xiong et al.,, 2013),

$\displaystyle\begin{cases}&f_{1}(\mathbf{x})=f_{2}(\mathbf{x})+\frac{\sin(x_{1})}{10}f_{2}(\mathbf{x})-2x_{1}+x_{2}^{2}+x_{3}^{2}+0.5\\\
&f_{2}(\mathbf{x})=\frac{x_{1}}{2}\left[\sqrt{1+(x_{2}+x_{3}^{2})\frac{x_{4}}{x_{1}^{2}}}-1\right]+(x_{1}+3x_{4})\exp{(1+\sin(x_{3}))}\end{cases}\text{for
}\mathbf{x}\in[0,1]^{4}.$

Furthermore, we explore the two-level Borehole function (Morris et al.,, 1993;
Xiong et al.,, 2013) and Currin function (Xiong et al.,, 2013), and introduce
a three-level function using the Franke function (Franke,, 1979). The detailed
forms of these three functions are available in the Appendix C.

The data are generated by evaluating these function at input locations
obtained from the nested space filling design introduced by Le Gratiet and
Garnier, (2014) with sample sizes $\\{n_{l}\\}^{L}_{l=1}$.The sample sizes and
input dimension for each example are outlined in Table 1. To examine the
prediction performance, $n_{\rm{test}}=1000$ random test input locations are
generated from the same input space. We evaluate the prediction performance
based on two criteria: the root-mean-squre error (RMSE) and continuous rank
probability score (CRPS) (Gneiting and Raftery,, 2007), which are defined as
follows:

RMSE
$\displaystyle=\left(\sum_{i=1}^{n_{\rm{test}}}\frac{\left(f_{L}(\mathbf{x}^{\rm{test}}_{i})-\mu^{*}_{L}(\mathbf{x}^{\rm{test}}_{i})\right)^{2}}{n_{\rm{test}}}\right)^{1/2},\quad\text{and}$
CRPS
$\displaystyle=\frac{1}{n_{\rm{test}}}\sum^{n_{\rm{test}}}_{i=1}\sigma_{L}^{*}(\mathbf{x}^{\rm{test}}_{i})\left[\frac{1}{\sqrt{\pi}}-2\psi\left(z(\mathbf{x}^{\rm{test}}_{i})\right)-z(\mathbf{x}^{\rm{test}}_{i})\left(2\Psi(z(\mathbf{x}^{\rm{test}}_{i}))-1\right)\right],$

with
$z(\mathbf{x}^{\rm{test}}_{i})=(f_{L}(\mathbf{x}^{\rm{test}}_{i})-\mu^{*}_{L}(\mathbf{x}^{\rm{test}}_{i}))/\sigma_{L}^{*}(\mathbf{x}^{\rm{test}}_{i})$,
where $\mathbf{x}^{\rm{test}}_{i}$ is the $i$-th test input location, $\psi$
and $\Psi$ denote the probability density function and cumulative distribution
function of a standard normal distribution, respectively. Note that CRPS
serves as a performance metric for the posterior predictive distribution of a
scalar observation. Lower values for the RMSE and CRPS indicate better model
accuracy. Additionally, we assess the computational efficiency by comparing
the computation time needed for model fitting and predictions.

|  Perdikaris |  Branin |  Park |  Borehole |  Currin |  Franke
---|---|---|---|---|---|---
$d$ |  1 |  2 |  4 |  8 |  2 |  2
$n_{1}$ |  13 |  20 |  40 |  60 |  20 |  20
$n_{2}$ |  8 |  15 |  20 |  30 |  10 |  15
$n_{3}$ |  |  10 |  |  |  |  10
Table 1: Dimension and samples sizes for each of the synthetic examples.
Figure 4: RMSEs of six synthetic examples across 100 repetitions.

Figures 4 and 5 show the results of RMSE and CRPS metrics across 100
repetitions, each employing a different random nested design for the training
input locations. The proposed RNAmf consistently outperforms CoKriging by both
metrics, particularly for examples exhibiting nonlinear relationships between
simulators, such as the Perdikaris, Borehole, Currin, and Franke functions.
For instances where simulators follow a linear (or nearly linear) auto-
regressive model, like the Brainin and Park functions, the proposed RNAmf
remains competitive with CoKriging, which is designed to excel in such
scenarios. On the other hand, NARGP, another approach modeling nonlinear
relationships, outperforms CoKriging in most of the examples and is
competitive with RNAmf, except in the Perdikaris and Franke examples, where
RNAmf exhibits superior performance. However, it comes with significantly
higher computational costs, as shown in Figure 6, due to its expensive MC
approximation, being roughly fifty times slower than both RNAmf and CoKriging
on average. Notably, in scenarios involving three fidelities, including the
Brainin and Franke examples, the computational time for NARGP exceeds that of
RNAmf by more than 200 times. This shows that NARGP can suffer from intensive
computation as the number of fidelity levels increases, while our method
remains competitive in this regard. In summary, the performance across these
synthetic examples underscores the capability of the proposed method in
providing an accurate emulator at a reasonable computational time, facilitated
by the closed-form posterior mean and variance derived in Proposition 3.1.

Figure 5: CRPSs of six synthetic examples across 100 repetitions. Figure 6:
Computational time of six synthetic functions across 100 repetitions.

###  5.2 Active learning performance

With the accurate RNA emulator in place, we now investigate on the performance
of AL strategies for the emulator using the proposed criteria. We compare with
two existing methods: CoKriging-CV, a cokriging-based sequential design
utilizing cross-validation techniques (Le Gratiet and Cannamela,, 2015), and
MR-SUR, a sequential design maximizing the rate of stepwise uncertainty
reduction using the KO model (Stroh et al.,, 2022). As for implementing
CoKriging-CV, we utilized the code provided in the Supplementary Materials of
Le Gratiet and Cannamela, (2015). Notably, both of these methods employed the
(linear) autoregressive model as in (2) in their implementations. To maintain
a consistent comparison, we use the one-dimensional Perdikaris function
(nonlinear) and the 4-dimensional Park function (linear autoregressive) in
Section 5.1, to illustrate the performance of these methods.

In this experiment, we suppose that the simulation costs associated with the
low- and high-fidelity simulators are $C_{1}=1$ and $C_{2}=3$, respectively.
The initial data is established similar to Section 5.1, with sample sizes
specified in Table 1. We consider a total simulation budget of
$C_{\rm{total}}=80$ for the Perdikaris function and $C_{\rm{total}}=130$ for
the Park function. For ALC and ALMC acquisitions which involve the computation
of the average reduction in variance as in (16), 1000 and 100 uniform samples
are respectively generated from the input space to approximate the integral
and impute the future outputs.

Figure 7 shows the results of RMSE and CRPS metrics for the Perdikaris
function, with respect to the total simulation costs accrued after each sample
selection. Additionally, the left panel of Figure 8 displays a boxplot
depicting the final RMSE at a total simulation cost of $C_{\rm{total}}=80$
across the 10 repetitions.

Figure 7: RMSE and CRPS for the Perdikaris function with respect to the
simulation cost. Solid lines represent the average over 10 repetitions and
shaded regions represent the ranges. Figure 8: Final RMSE (left) and
proportion of AL acquisitions choosing low-fidelity data (right) for the
Perdikaris function. Boxplots indicate spread over 10 repetitions.

The results show that the proposed AL strategies dramatically outperform the
two competing methods, CoKriging-CV and MR-SUR, in terms of both accuracy and
stability, considering the same costs. As the cost increases, MR-SUR begins to
close the gap, while CoKriging-CV lags behind the other methods. Among the
three proposed AL strategies, the distinctions are minimal, although ALC
appears to be less accurate than ALM and ALMC as the cost increases. This
implies that, in practice, it may be beneficial to choose between ALM and
ALMC, particularly since ALC is computationally more intensive. As noted in
Section 4, ALC acquisitions involve intricate numerical integration
approximations and data imputation, taking approximately 400 seconds for each
acquisition in this example. In contrast, ALM and ALMC are significantly more
computationally efficient due to the closed-form nature of the posterior
variance, requiring only around 0.8 seconds and 10 seconds per acquisition,
respectively.

From the right panel of Figure 8, it can be seen that the proposed AL
strategies tend to select low-fidelity simulators more frequently than the
other two comparative methods, notably MR-SUR, which consistently chooses
samples exclusively from the high-fidelity simulator. This suggests that the
proposed RNA model can effectively infer the high-fidelity simulation using
primarily low-fidelity data for the nonlinear Perdikaris function, while the
other two KO-based methods (CoKriging-CV and MR-SUR) require more high-
fidelity data to reduce the uncertainty.

Figures 9 and 10 present the results for the Park function. As expected, the
distinctions between these strategies are not as significant because the
function aligns more closely with the KO model (linear autoregressive).
Nonetheless, our proposed AL strategies still exhibit better average
performance. At the final cost budget of $C_{\rm{total}}=130$, ALMC performs
the best, collecting a larger portion of high-fidelity data, as indicated in
Figure 10. In contrast, the KO-based strategies collect more low-fidelity
data, which is again expected because KO-based models are efficient at
leveraging low-fidelity data to infer the high-fidelity simulator. In these
scenarios, our strategies efficiently prioritize the selection of high-
fidelity data to minimize the uncertainty, resulting in superior prediction
accuracy at the same cost.

Figure 9: RMSE and CRPS for the Park function with respect to the simulation
cost. Solid lines represent the average over 10 repetitions and shaded regions
represent the ranges. Figure 10: Final RMSE (left) and proportion of AL
acquisitions choosing low-fidelity data (right) for the Park function.
Boxplots indicate spread over 10 repetitions.

##  6 Thermal Stress Analysis of Jet Engine Turbine Blade

We leverage our proposed method for a real application involving the analysis
of thermal stress in a jet turbine engine blade under steady-state operating
conditions. The turbine blade, which forms part of the jet engine, is
constructed from nickel alloys capable of withstanding extremely high
temperatures. It is crucial for the blade’s design to ensure that it can
endure stress and deformations while avoiding mechanical failure and friction
between the blade tip and the turbine casing. We refer more details to Carter,
(2005), Wright and Han, (2006), Sung et al., (2022) and Sung et al., (2023).

This problem can be treated as a static structural model and can be solved
numerically using finite element methods. There are two input variables
denoted as $x_{1}$ and $x_{2}$, which represent the pressure load on the
pressure and suction sides of the blade, both of which fall within the range
of 0.25 to 0.75 MPa, i.e.,
$\mathbf{x}=(x_{1},x_{2})\in\Omega=[0.25,0.75]^{2}$. The response of interest
is the maximum value over the thermal stress profile, which is a critical
parameter used to assess the structural stability of the turbine blade. We
perform finite element simulations using the Partial Differential Equation
Toolbox in MATLAB (MATLAB,, 2021).

The simulations are conducted at two fidelity levels, each using different
mesh densities for finite element methods. A denser mesh provides higher
fidelity and more accurate results but demands greater computational
resources. Conversely, a coarser mesh sacrifices some accuracy for reduced
computational cost. Figure 11 visually presents the turbine blade structure
and thermal stress profiles obtained at these two fidelity levels for the
input location $\mathbf{x}=(0.5,0.45)$.

Figure 11:  Illustration of low-fidelity (left) and high-fidelity (right)
finite element simulations at the input setting $\mathbf{x}=(0.5,0.45)$ in the
turbine blade application.

We perform the finite element simulations with sample sizes of $n_{1}=20$ and
$n_{2}=10$ to examine the emulation performance. Similar to Section 5.1, we
use the nested space-filling designs of Le Gratiet and Garnier, (2014) to
generate the input locations of the two-fidelity computer experiments. We
record the simulation time of the finite element simulations, which are
respectively $C_{1}=2.25$ and $C_{2}=6.85$ (seconds) and will be used later
for comparing AL strategies. To examine the performance, we conduct the high-
fidelity simulations (i.e. $f_{2}(\mathbf{x})$) at the test input locations of
size $n_{\text{test}}=100$ generated from a set of Latin hypercube samples
from the same design space. The experiment is repeated 10 times, each time
considering different nested space-filling designs for the training input
locations.

Figure 12 presents a comparison of emulation performance with the other two
competing methods, CoKriging and NARGP. Our proposed method, RNAmf,
outperforms the other two methods in terms of RMSE and CRPS. While NARGP
delivers competitive prediction performance, it comes at a significantly
higher computational cost compared to RNAmf.

Figure 12: RMSE, CRPS, and computation time across 10 repetitions in the
turbine blade application.

Figures 13 and 14 present a comparison of the AL strategies with a fixed cost
budget of $C_{\rm{total}}=160$ seconds. The right panel of Figure 14 reveals
that these strategies collect a similar number of low-fidelity data points.
Notably, CoKriging-CV exhibits significant variability across the 10
repetitions and yields poorer prediction performance compared to the other
strategies. Another KO-based strategy, MR-SUR, performs better but still falls
short of our proposed AL strategies at any given simulation cost. Conversely,
our proposed AL strategies demonstrate remarkably stable results and
outperform the others, with ALC being slightly superior to the other two at
the final cost budget.

Figure 13: RMSE and CRPS for the turbine blade application with respect to the
cost. Solid lines represent the average over 10 repetitions and shaded regions
represent the ranges. Figure 14: Final RMSE (left) and proportion of AL
acquisitions choosing low-fidelity data (right) for the turbine blade
application. Boxplots indicate spread over 10 repetitions.

##  7 Conclusion

Multi-fidelity computer experiments have become an essential tool in
simulating complex scientific problems. This paper introduces a new emulator,
the RNA emulator, tailored for multi-fidelity simulations, which proves
effective in producing accurate, efficient predictions for high-fidelity
simulations, especially when dealing with nonlinear relationships between
simulators. Building upon this new emulator, we present three AL strategies
designed to select optimal input locations and fidelity levels to augment
data, thereby enhancing emulation performance. Our numerical studies,
including a jet turbine blade analysis, underscores the effectiveness of our
approach in efficiently emulating multi-fidelity computer experiments and
ensuring accurate emulation through our AL strategies.

With the RNA emulator’s success, it is worthwhile to explore emulators and AL
strategies built upon similar principles for addressing multi-fidelity
problems with tunable fidelity parameters, such as mesh density (Picheny et
al.,, 2013; Tuo et al.,, 2014). Designing experiments for such scenarios
presents intriguing challenges, as shown in recent studies (see, e.g., Shaowu
Yuchi et al., (2023) and Sung et al., (2022)). Furthermore, considering the
increasing prevalence of stochastic computer models (Baker et al.,, 2022),
extending the proposed RNA emulator to accommodate noisy data would
significantly enhance its relevance in real-world applications. While this
article assumes noise-free data, introducing noise into the model is a
feasible endeavor, a task we leave for our future research..

Supplemental Materials Additional supporting materials can be found in
Supplemental Materials, including the closed-form posterior mean and variance
under a squared exponential kernel and a Matérn kernel with the smoothness
parameter $\nu=1.5$ and $\nu=2.5$. The R code and package for reproducing the
results in Sections 5 and 6 are also provided.

##  Appendix A Matérn kernel functions

Here are the Matérn kernels with smoothness parameters of 1.5 and 2.5, both of
which are widely used and come with simpler expressions as follows:

$K_{l}(\mathbf{z},\mathbf{z}^{\prime})=\phi(y,y^{\prime};\theta_{ly})\prod^{d}_{j=1}\phi(x_{j},x^{\prime}_{j};\theta_{lj})$

with

$\phi(x,x^{\prime};\theta)=\left(1+\frac{\sqrt{3}|x-x^{\prime}|}{\theta}\right)\exp\left(-\frac{\sqrt{3}|x-x^{\prime}|}{\theta}\right)$

for smoothness parameter of 1.5 and

$\phi(x,x^{\prime};\theta)=\left(1+\frac{\sqrt{5}|x-x^{\prime}|}{\theta}+\frac{5(x-x^{\prime})^{2}}{3\theta^{2}}\right)\exp\left(-\frac{\sqrt{5}|x-x^{\prime}|}{\theta}\right)$

for smoothness parameter of 2.5.

##  Appendix B Proof of Proposition 3.2

For $l=1$, since $\mu^{*}_{1}(\mathbf{x})=\mu_{1}(\mathbf{x})$ and
$\sigma^{*2}_{1}(\mathbf{x})=\sigma^{2}_{1}(\mathbf{x})$ are the posterior
mean and variance of a conventional GP, it can be shown that
$\mu^{*}_{1}(\mathbf{x}^{[1]}_{i})=y^{[1]}_{i}$ and
$\sigma^{*2}_{1}(\mathbf{x}^{[1]}_{i})=0$ (Santner et al.,, 2018). For
$l=k-1$, suppose that $\mu^{*}_{k-1}(\mathbf{x}^{[k-1]}_{i})=y^{[k-1]}_{i}$
and $\sigma^{*2}_{k-1}(\mathbf{x}^{[k-1]}_{i})=0$, implying that
$f_{k-1}(\mathbf{x}^{[k-1]}_{i})|\mathbf{y}_{1},\ldots,\mathbf{y}_{k-1}$
remains constant at the value $y^{[k-1]}_{i}$. Then,

$\displaystyle\mu^{*}_{k}(\mathbf{x}_{i}^{[k]})=$
$\displaystyle\mathbb{E}[f_{k}(x_{i}^{[k]})|\mathbf{y}_{1},\ldots,\mathbf{y}_{k}]=\mathbb{E}[\mathbb{E}[f_{k}(\mathbf{x}_{i}^{[k]})|f_{k-1}(\mathbf{x}_{i}^{[k]}),\mathbf{y}_{1},\ldots,\mathbf{y}_{k}]]$
$\displaystyle=$
$\displaystyle\mathbb{E}[\mu_{k}(\mathbf{x}^{[k]}_{i},f_{k-1}(\mathbf{x}_{i}^{[k]}))|\mathbf{y}_{1},\ldots,\mathbf{y}_{k}]$
$\displaystyle=$
$\displaystyle\mu_{k}(\mathbf{x}^{[k]}_{i},y^{[k-1]}_{i})=y^{[k]}_{i},$

and

$\displaystyle\sigma^{*2}_{k}(\mathbf{x}_{i}^{[k]})$
$\displaystyle=\mathbb{V}[f_{k}(\mathbf{x}_{i}^{[k]})|\mathbf{y}_{1},\ldots,\mathbf{y}_{k}]$
$\displaystyle=\mathbb{V}[\mathbb{E}[f_{k}(\mathbf{x}_{i}^{[k]})|f_{k-1}(\mathbf{x}_{i}^{[k]}),\mathbf{y}_{1},\ldots,\mathbf{y}_{k}]]+\mathbb{E}[\mathbb{V}[f_{k}(\mathbf{x}_{i}^{[k]})|f_{k-1}(\mathbf{x}_{i}^{[k]}),\mathbf{y}_{1},\ldots,\mathbf{y}_{k}]]$
$\displaystyle=\mathbb{V}[\mu_{k}(\mathbf{x}^{[k]}_{i},f_{k-1}(\mathbf{x}_{i}^{[k]}))|\mathbf{y}_{1},\ldots,\mathbf{y}_{k}]+\mathbb{E}[\sigma^{2}_{k}(\mathbf{x}^{[k]}_{i},f_{k-1}(\mathbf{x}_{i}^{[k]}))|\mathbf{y}_{1},\ldots,\mathbf{y}_{k}]$
$\displaystyle=0+\sigma^{2}_{k}(\mathbf{x}^{[k]}_{i},y^{[k-1]}_{i})=0,$

where $\mu^{2}_{k}(\mathbf{x},f_{k-1}(\mathbf{x}))$ and
$\sigma^{2}_{k}(\mathbf{x},f_{k-1}(\mathbf{x}))$ are defined in (10) and (11).
In both of the derivations, the second to last equation holds because of
$\mathbf{x}^{[k]}_{i}=\mathbf{x}^{[k-1]}_{i}$ and
$f_{k-1}(\mathbf{x}^{[k-1]}_{i})|\mathbf{y}_{1},\ldots,\mathbf{y}_{k-1}:=y^{[k]}_{i}$.
The last equation holds because of the interpolation property of conventional
GPs. By induction, this finishes the proof.

##  Appendix C Synthetic functions in Section 5.1

###  C.1 Two-level borehole function

$\displaystyle\begin{cases}&f_{1}(\mathbf{x})=\frac{2\pi
T_{u}(H_{u}-H_{l})}{\ln(r/r_{w})(1+\frac{2LT_{u}}{\ln(r/r_{w})r^{2}_{w}K_{w}}+\frac{T_{u}}{T_{l}})}\\\
&f_{2}(\mathbf{x})=\frac{5T_{u}(H_{u}-H_{l})}{\ln(r/r_{w})(1.5+\frac{2LT_{u}}{\ln(r/r_{w})r^{2}_{w}K_{w}}+\frac{T_{u}}{T_{l}})},\end{cases}$

where
$r_{w}\in[0.05,0.15],r\in[100,50000],T_{u}\in[63070,115600],H_{u}\in[990,1110],T_{l}\in[63.1,116],H_{l}\in[700,820],L\in[1120,1680]$
and $K_{w}\in[9855,12045]$.

###  C.2 Two-level Currin function

$\displaystyle\begin{cases}&f_{1}(\mathbf{x})=\frac{1}{4}\left[f_{2}(x_{1}+0.05,x_{2}+0.05)+f_{2}(x_{1}+0.05,\max(0,x_{2}-0.05))\right]\\\
&\quad\quad\quad+\frac{1}{4}\left[f_{2}(x_{1}-0.05,x_{2}+0.05)+f_{2}(x_{1}-0.05,\max(0,x_{2}-0.05))\right]\\\
&f_{2}(\mathbf{x})=\left[1-\exp\left(-\frac{1}{2x_{2}}\right)\right]\frac{2300x_{1}^{3}+1900x_{1}^{2}+2092x_{1}+60}{100x_{1}^{3}+500x_{1}^{2}+4x_{1}+20}\end{cases}\text{for
}\mathbf{x}\in[0,1]^{2}.$

###  C.3 Three-level Franke function

$\displaystyle\begin{cases}&f_{1}(\mathbf{x})=0.75\exp\left(-\frac{(9x_{1}-2)^{2}}{4}-\frac{(9x_{2}-2)^{2}}{4}\right)+0.75\exp\left(-\frac{(9x_{1}+1)^{2}}{49}-\frac{9x_{2}+1}{10}\right)\\\
&\quad\quad+0.5\exp\left(-\frac{(9x_{1}-7)^{2}}{4}-\frac{(9x_{2}-3)^{2}}{4}\right)-0.2\exp\left(-(9x_{1}-4)^{2}-(9x_{2}-7)^{2}\right)\\\
&f_{2}(\mathbf{x})=\exp(-1.4f_{1}(\mathbf{x}))\cos(3.5\pi
f_{1}(\mathbf{x}))\\\
&f_{3}(\mathbf{x})=\sin(2\pi(f_{2}(\mathbf{x})-1))\end{cases}\text{for
}\mathbf{x}\in[0,1]^{2}.$

## References

  * Baker et al., (2022) Baker, E., Barbillon, P., Fadikar, A., Gramacy, R. B., Herbei, R., Higdon, D., Huang, J., Johnson, L. R., Ma, P., Mondal, A., Pires, B., Sacks, J., and Sokolov, V. (2022).  Analyzing stochastic computer models: A review with opportunities.  Statistical Science, 37(1):64–89.
  * Binois et al., (2019) Binois, M., Huang, J., Gramacy, R. B., and Ludkovski, M. (2019).  Replication or exploration? sequential design for stochastic simulation experiments.  Technometrics, 61(1):7–23.
  * Byrd et al., (1995) Byrd, R. H., Lu, P., Nocedal, J., and Zhu, C. (1995).  A limited memory algorithm for bound constrained optimization.  SIAM Journal on Scientific Computing, 16(5):1190–1208.
  * Carter, (2005) Carter, T. J. (2005).  Common failures in gas turbine blades.  Engineering Failure Analysis, 12(2):237–247.
  * Cohn, (1993) Cohn, D. (1993).  Neural network exploration using optimal experiment design.  Advances in Neural Information Processing Systems, 6:1071–1083.
  * Demeyer et al., (2017) Demeyer, S., Fischer, N., and Marquis, D. (2017).  Surrogate model based sequential sampling estimation of conformance probability for computationally expensive systems: application to fire safety science.  Journal de la société française de statistique, 158(1):111–138.
  * Ehara and Guillas, (2023) Ehara, A. and Guillas, S. (2023).  An adaptive strategy for sequential designs of multilevel computer experiments.  International Journal for Uncertainty Quantification, 13(4):61–98.
  * Franke, (1979) Franke, R. (1979).  A critical comparison of some methods for interpolation of scattered data.  Technical report, Monterey, California: Naval Postgraduate School.
  * Gneiting and Raftery, (2007) Gneiting, T. and Raftery, A. E. (2007).  Strictly proper scoring rules, prediction, and estimation.  Journal of the American Statistical Association, 102(477):359–378.
  * Gramacy, (2020) Gramacy, R. B. (2020).  Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences.  CRC press.
  * Haaland and Qian, (2010) Haaland, B. and Qian, P. Z. G. (2010).  An approach to constructing nested space-filling designs for multi-fidelity computer experiments.  Statistica Sinica, 20(3):1063–1075.
  * Harville, (1998) Harville, D. A. (1998).  Matrix Algebra from a Statistician’s Perspective.  New York: Springer-Verlag.
  * He et al., (2017) He, X., Tuo, R., and Wu, C. F. J. (2017).  Optimization of multi-fidelity computer experiments via the EQIE criterion.  Technometrics, 59(1):58–68.
  * Huang et al., (2006) Huang, D., Allen, T. T., Notz, W. I., and Miller, R. A. (2006).  Sequential kriging optimization using multiple-fidelity evaluations.  Structural and Multidisciplinary Optimization, 32:369–382.
  * Ji et al., (2022) Ji, Y., Mak, S., Soeder, D., Paquet, J.-F., and Bass, S. A. (2022).  A graphical multi-fidelity Gaussian process model, with application to emulation of expensive computer simulations.  arXiv preprint arXiv:2108.00306.
  * Kennedy and O’Hagan, (2000) Kennedy, M. C. and O’Hagan, A. (2000).  Predicting the output from a complex computer code when fast approximations are available.  Biometrika, 87(1):1–13.
  * Kerleguer et al., (2024) Kerleguer, B., Cannamela, C., and Garnier, J. (2024).  A bayesian neural network approach to multi-fidelity surrogate modelling.  International Journal for Uncertainty Quantification, 14(1):43–60.
  * Koermer et al., (2023) Koermer, S., Loda, J., Noble, A., and Gramacy, R. B. (2023).  Active learning for simulator calibration.  arXiv preprint arXiv:2301.10228.
  * Kuya et al., (2011) Kuya, Y., Takeda, K., Zhang, X., and Forrester, A. I. (2011).  Multifidelity surrogate modeling of experimental and computational aerodynamic data sets.  AIAA Journal, 49(2):289–298.
  * Kyzyurova et al., (2018) Kyzyurova, K. N., Berger, J. O., and Wolpert, R. L. (2018).  Coupling computer models through linking their statistical emulators.  SIAM/ASA Journal on Uncertainty Quantification, 6(3):1151–1171.
  * Le Gratiet, (2012) Le Gratiet, L. (2012).  MuFiCokriging: Multi-Fidelity Cokriging models.  R package version 1.2.
  * Le Gratiet, (2013) Le Gratiet, L. (2013).  Bayesian analysis of hierarchical multifidelity codes.  SIAM/ASA Journal on Uncertainty Quantification, 1(1):244–269.
  * Le Gratiet and Cannamela, (2015) Le Gratiet, L. and Cannamela, C. (2015).  Cokriging-based sequential design strategies using fast cross-validation techniques for multi-fidelity computer codes.  Technometrics, 57(3):418–427.
  * Le Gratiet and Garnier, (2014) Le Gratiet, L. and Garnier, J. (2014).  Recursive co-kriging model for design of computer experiments with multiple levels of fidelity.  International Journal for Uncertainty Quantification, 4(5):365–386.
  * Li et al., (2020) Li, S., Xing, W., Kirby, R., and Zhe, S. (2020).  Multi-fidelity Bayesian optimization via deep neural networks.  Advances in Neural Information Processing Systems, 33:8521–8531.
  * Ma et al., (2022) Ma, P., Karagiannis, G., Konomi, B. A., Asher, T. G., Toro, G. R., and Cox, A. T. (2022).  Multifidelity computer model emulation with high-dimensional output: An application to storm surge.  Journal of the Royal Statistical Society Series C: Applied Statistics, 71(4):861–883.
  * MacKay, (1992) MacKay, D. J. C. (1992).  Information-based objective functions for active data selection.  Neural Computation, 4(4):590–604.
  * Mak et al., (2018) Mak, S., Sung, C.-L., Wang, X., Yeh, S.-T., Chang, Y.-H., Joseph, V. R., Yang, V., and Wu, C. F. J. (2018).  An efficient surrogate model for emulation and physics extraction of large eddy simulations.  Journal of the American Statistical Association, 113(524):1443–1456.
  * MATLAB, (2021) MATLAB (2021).  MATLAB version 9.11.0.1769968 (R2021b).  The Mathworks, Inc., Natick, Massachusetts.
  * McKay et al., (2000) McKay, M. D., Beckman, R. J., and Conover, W. J. (2000).  A comparison of three methods for selecting values of input variables in the analysis of output from a computer code.  Technometrics, 42(1):55–61.
  * Meng et al., (2021) Meng, X., Babaee, H., and Karniadakis, G. E. (2021).  Multi-fidelity Bayesian neural networks: Algorithms and applications.  Journal of Computational Physics, 438:110361.
  * Meng and Karniadakis, (2020) Meng, X. and Karniadakis, G. E. (2020).  A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse PDE problems.  Journal of Computational Physics, 401:109020.
  * Ming and Guillas, (2021) Ming, D. and Guillas, S. (2021).  Linked Gaussian process emulation for systems of computer models using Matérn kernels and adaptive design.  SIAM/ASA Journal on Uncertainty Quantification, 9(4):1615–1642.
  * Ming et al., (2023) Ming, D., Williamson, D., and Guillas, S. (2023).  Deep Gaussian process emulation using stochastic imputation.  Technometrics, 65(2):150–161.
  * Morris et al., (1993) Morris, M. D., Mitchell, T. J., and Ylvisaker, D. (1993).  Bayesian design and analysis of computer experiments: use of derivatives in surface prediction.  Technometrics, 35(3):243–255.
  * Park et al., (2023) Park, C., Waelder, R., Kang, B., Maruyama, B., Hong, S., and Gramacy, R. B. (2023).  Active learning of piecewise Gaussian process surrogates.  arXiv preprint arXiv:2301.08789.
  * Park, (1991) Park, J. S. (1991).  Tuning Complex Computer Codes to Data and Optimal Designs.  University of Illinois at Urbana-Champaign.
  * Patra et al., (2020) Patra, A., Batra, R., Chandrasekaran, A., Kim, C., Huan, T. D., and Ramprasad, R. (2020).  A multi-fidelity information-fusion approach to machine learn and predict polymer bandgap.  Computational Materials Science, 172:109286.
  * Perdikaris, (2016) Perdikaris, P. (2016).  Multi-fidelity modeling using Gaussian processes and nonlinear auto-regressive schemes.  https://github.com/paraklas/NARGP.
  * Perdikaris et al., (2017) Perdikaris, P., Raissi, M., Damianou, A., Lawrence, N. D., and Karniadakis, G. E. (2017).  Nonlinear information fusion algorithms for data-efficient multi-fidelity modelling.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2198):20160751.
  * Picheny et al., (2013) Picheny, V., Ginsbourger, D., Richet, Y., and Caplin, G. (2013).  Quantile-based optimization of noisy computer experiments with tunable precision.  Technometrics, 55(1):2–13.
  * Qian, (2009) Qian, P. Z. G. (2009).  Nested latin hypercube designs.  Biometrika, 96(4):957–970.
  * Qian et al., (2009) Qian, P. Z. G., Ai, M., and Wu, C. F. J. (2009).  Construction of nested space-filling designs.  Annals of Statistics, 37(6A):3616–3643.
  * Qian and Wu, (2008) Qian, P. Z. G. and Wu, C. F. J. (2008).  Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments.  Technometrics, 50(2):192–204.
  * Qian et al., (2006) Qian, Z., Seepersad, C. C., Joseph, V. R., Allen, J. K., and Wu, C. F. J. (2006).  Building surrogate models based on detailed and approximate simulations.  Journal of Mechanical Design, 128(4):668–677.
  * R Core Team, (2018) R Core Team (2018).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Rasmussen and Williams, (2006) Rasmussen, C. E. and Williams, C. K. (2006).  Gaussian Processes for Machine Learning.  Cambridge, MA: MIT Press.
  * Santner et al., (2018) Santner, T. J., Williams, B. J., and Notz, W. I. (2018).  The Design and Analysis of Computer Experiments.  Springer New York.
  * Sauer et al., (2023) Sauer, A., Gramacy, R. B., and Higdon, D. (2023).  Active learning for deep Gaussian process surrogates.  Technometrics, 65(1):4–18.
  * Seo et al., (2000) Seo, S., Wallat, M., Graepel, T., and Obermayer, K. (2000).  Gaussian process regression: Active data selection and test point rejection.  In Mustererkennung 2000: 22. DAGM-Symposium. Kiel, 13.–15. September 2000, pages 27–34. Springer.
  * Shaowu Yuchi et al., (2023) Shaowu Yuchi, H., Roshan Joseph, V., and Wu, C. F. J. (2023).  Design and analysis of multifidelity finite element simulations.  Journal of Mechanical Design, 145(6):061703.
  * Sobester et al., (2008) Sobester, A., Forrester, A., and Keane, A. (2008).  Engineering Design via Surrogate Modelling: A Practical Guide.  John Wiley & Sons.
  * Stein, (1999) Stein, M. L. (1999).  Interpolation of Spatial Data: Some Theory for Kriging.  Springer Science & Business Media.
  * Stroh et al., (2022) Stroh, R., Bect, J., Demeyer, S., Fischer, N., Marquis, D., and Vazquez, E. (2022).  Sequential design of multi-fidelity computer experiments: maximizing the rate of stepwise uncertainty reduction.  Technometrics, 64(2):199–209.
  * Sung et al., (2020) Sung, C.-L., Hung, Y., Rittase, W., Zhu, C., and Wu, C. F. J. (2020).  Calibration for computer experiments with binary responses and application to cell adhesion study.  Journal of the American Statistical Association, 115(532):1664–1674.
  * Sung et al., (2022) Sung, C.-L., Ji, Y., Tang, T., and Mak, S. (2022).  Stacking designs: designing multi-fidelity computer experiments with target predictive accuracy.  arXiv preprint arXiv:2211.00268.
  * Sung et al., (2023) Sung, C.-L., Wang, W., Ding, L., and Wang, X. (2023).  Mesh-clustered Gaussian process emulator for partial differential equation systems.  arXiv preprint arXiv:2301.10387.
  * Swersky et al., (2013) Swersky, K., Snoek, J., and Adams, R. P. (2013).  Multi-task Bayesian optimization.  Advances in Neural Information Processing Systems, 26:2004–2012.
  * Tuo et al., (2014) Tuo, R., Wu, C. F. J., and Yu, D. (2014).  Surrogate modeling of computer experiments with different mesh densities.  Technometrics, 56(3):372–380.
  * Wright and Han, (2006) Wright, L. M. and Han, J.-C. (2006).  Enhanced internal cooling of turbine blades and vanes.  The Gas Turbine Handbook, 4:1–5.
  * Xiong et al., (2013) Xiong, S., Qian, P. Z. G., and Wu, C. F. J. (2013).  Sequential design and analysis of high-accuracy and low-accuracy computer codes.  Technometrics, 55(1):37–46.

Supplementary Materials for “Active Learning for a Recursive Non-Additive
Emulator for Multi-Fidelity Computer Experiments”

##  S1 Proof of Proposition 3.1

For notational simplicity, we denote
$\mathbf{Y}_{l}=\\{\mathbf{y}_{1},\ldots,\mathbf{y}_{l}\\}$. Based on the
squared exponential kernel, the posterior mean and variance at the input
$\mathbf{x}$ can be derived as follows,

$\displaystyle\mu_{l}^{*}(\mathbf{x})$
$\displaystyle=\mathbb{E}[f_{l}(\mathbf{x})|\mathbf{Y}_{l}]$
$\displaystyle=\mathbb{E}[\mathbb{E}[f_{l}(\mathbf{x})|f_{l-1}(\mathbf{x}),\mathbf{Y}_{l}]]$
$\displaystyle=\alpha_{l}+\mathbb{E}[\mathbf{k}^{T}_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))|\mathbf{Y}_{l}]\mathbf{K}^{-1}_{l}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}})$
$\displaystyle=\alpha_{l}+\sum^{n_{l}}_{i=1}r_{i}\prod_{j=1}^{d}\exp\left(-\frac{(x_{j}-x^{[l]}_{ij})^{2}}{\theta_{lj}}\right)\mathbb{E}\left[\exp{\left\\{-\frac{(y_{i}^{[l-1]}-f_{l-1}(\mathbf{x}))^{2}}{\theta_{ly}}\right\\}}\Big{|}\mathbf{Y}_{l}\right]$
$\displaystyle=\alpha_{l}+\sum^{n_{l}}_{i=1}r_{i}\prod_{j=1}^{d}\exp\left(-\frac{(x_{j}-x^{[l]}_{ij})^{2}}{\theta_{lj}}\right)\frac{1}{\sqrt{1+2\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}}}\exp{\left(-\frac{(y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x}))^{2}}{\theta_{ly}+2\sigma^{*2}_{l-1}(\mathbf{x})}\right)},$
and $\displaystyle\sigma_{l}^{*2}(\mathbf{x})$
$\displaystyle=\mathbb{V}[f_{l}(\mathbf{x})|\mathbf{Y}_{l}]$
$\displaystyle=\mathbb{V}\left[\mathbb{E}[f_{l}(\mathbf{x})|f_{l-1}(\mathbf{x}),\mathbf{Y}_{l}]\right]+\mathbb{E}\left[\mathbb{V}[f_{l}(\mathbf{x})|f_{l-1}(\mathbf{x}),\mathbf{Y}_{l}]\right]$
$\displaystyle=\mathbb{E}\left[\left\\{\mathbb{E}\left[f_{l}(\mathbf{x})|f_{l-1}(\mathbf{x}),\mathbf{Y}_{l}\right]\right\\}^{2}\right]-\mu^{*}_{l}(\mathbf{x})^{2}+\mathbb{E}\left[\mathbb{V}\left[f_{l}(\mathbf{x})|f_{l-1}(\mathbf{x}),\mathbf{Y}_{l}\right]\right]$
$\displaystyle=\mathbb{E}\left[\left\\{\alpha_{l}+\mathbf{k}_{l}^{T}(\mathbf{x},f_{l-1}(\mathbf{x}))\mathbf{K}^{-1}_{l}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}})\right\\}^{2}|\mathbf{Y}_{l}\right]-\mu^{*}_{l}(\mathbf{x})^{2}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad+\mathbb{E}\left[\tau_{l}^{2}\left\\{1-\mathbf{k}_{l}^{T}(\mathbf{x},f_{l-1}(\mathbf{x}))\mathbf{K}^{-1}_{l}\mathbf{k}_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))\right\\}|\mathbf{Y}_{l}\right]$
$\displaystyle=\alpha_{l}^{2}+2\alpha_{l}(\mu^{*}_{l}(\mathbf{x})-\alpha_{l})+\mathbb{E}\left[\left\\{\mathbf{k}_{l}^{T}(\mathbf{x},f_{l-1}(\mathbf{x}))\mathbf{K}^{-1}_{l}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}})\right\\}^{2}|\mathbf{Y}_{l}\right]$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad-\mu^{*}_{l}(\mathbf{x})^{2}+\tau^{2}_{l}-\tau^{2}_{l}\mathbb{E}\left[\left\\{\mathbf{k}_{l}^{T}(\mathbf{x},f_{l-1}(\mathbf{x}))\mathbf{K}^{-1}_{l}\mathbf{k}_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))\right\\}|\mathbf{Y}_{l}\right]$
$\displaystyle=\tau^{2}_{l}-(\mu^{*}_{l}(\mathbf{x})-\alpha_{l})^{2}+\left(\sum_{i,k=1}^{n_{l}}\zeta_{ik}\left(r_{i}r_{k}-\tau^{2}_{l}(\mathbf{K}_{l}^{-1})_{ik}\right)\prod_{j=1}^{d}\exp{\left(-\frac{(x_{j}-{x}^{[l]}_{ij})^{2}+(x_{j}-{x}^{[l]}_{kj})^{2}}{\theta_{lj}}\right)}\right),$

where

$\displaystyle r_{i}$
$\displaystyle=(\mathbf{K}^{-1}_{l}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}}))_{i},$
$\displaystyle\zeta_{ik}$
$\displaystyle=\frac{1}{\sqrt{1+4\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}}}\exp{\left(-\frac{(\frac{y_{i}^{[l-1]}+y_{k}^{[l-1]}}{2}-\mu^{*}_{l-1}(\mathbf{x}))^{2}}{\frac{\theta_{ly}}{2}+2\sigma^{*2}_{l-1}(\mathbf{x})}-\frac{(y_{i}^{[l-1]}-y_{k}^{[l-1]})^{2}}{2\theta_{ly}}\right)}.$

##  S2 Posterior mean and variance under Matérn kernel

The posterior mean and variance at the input $\mathbf{x}$ can be derived as
follows,

$\displaystyle\mu^{*}_{l}(\mathbf{x})$
$\displaystyle=\mathbb{E}[f_{l}(\mathbf{x})|\mathbf{Y}_{l}]$
$\displaystyle=\alpha_{l}+\mathbb{E}[\mathbf{k}^{T}_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))|\mathbf{Y}_{l}]\mathbf{K}^{-1}_{l}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}})$
$\displaystyle=\alpha_{l}+\sum^{n_{l}}_{i=1}r_{i}\xi_{i}\prod^{d}_{j=1}\phi(x_{j},x^{[l]}_{ij};\theta_{lj}),$
$\displaystyle\sigma_{l}^{*2}(\mathbf{x})$
$\displaystyle=\mathbb{V}[f_{l}(\mathbf{x})|\mathbf{Y}_{l}]$
$\displaystyle=\alpha_{l}^{2}+2\alpha_{l}(\mu^{*}_{l}(\mathbf{x})-\alpha_{l})+\mathbb{E}\left[\left\\{\mathbf{k}_{l}^{T}(\mathbf{x},f_{l-1}(\mathbf{x}))\mathbf{K}^{-1}_{l}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}})\right\\}^{2}|\mathbf{Y}_{l}\right]$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad-\mu^{*}_{l}(\mathbf{x})^{2}+\tau^{2}_{l}-\tau^{2}_{l}\mathbb{E}\left[\left\\{\mathbf{k}_{l}^{T}(\mathbf{x},f_{l-1}(\mathbf{x}))\mathbf{K}^{-1}_{l}\mathbf{k}_{l}(\mathbf{x},f_{l-1}(\mathbf{x}))\right\\}|\mathbf{Y}_{l}\right]$
$\displaystyle=\tau^{2}_{l}-(\mu^{*}_{l}(\mathbf{x})-\alpha_{l})^{2}+\left(\sum_{i,k=1}^{n_{l}}\zeta_{ik}\left(r_{i}r_{k}-\tau^{2}_{l}(\mathbf{K}_{l}^{-1})_{ik}\right)\prod_{j=1}^{d}\phi(x_{j},x^{[l]}_{ij};\theta_{lj})\phi(x_{j},x^{[l]}_{kj};\theta_{lj})\right)$

where $\phi$ is defined in Appendix A,
$(r_{1},\ldots,r_{n_{l}})=\mathbf{K}^{-1}_{l}(\mathbf{y}_{l}-\alpha_{l}\mathbf{1}_{n_{l}}),\xi_{i}=\mathbb{E}\left[\phi(f_{l-1}(\mathbf{x}),y_{i}^{[l-1]};\theta_{ly})\Big{|}\mathbf{Y}_{l}\right]$,
and
$\zeta_{ik}=\mathbb{E}\left[\phi(f_{l-1}(\mathbf{x}),y_{i}^{[l-1]};\theta_{ly})\phi(f_{l-1}(\mathbf{x}),y_{k}^{[l-1]};\theta_{ly})\Big{|}\mathbf{Y}_{l}\right]$.
The closed-form expressions of $\xi_{i}$ and $\zeta_{ik}$ are provided in the
following subsections.

###  S2.1 Matérn-1.5 kernel

For Matérn-1.5 kernel,
$\Phi_{l}(\mathbf{x},\mathbf{x}^{\prime})=\prod_{j=1}^{d}\phi(x_{j},x^{\prime}_{j};\theta_{lj})=\prod_{j=1}^{d}\left(1+\frac{\sqrt{3}|x_{j}-x_{j}^{\prime}|}{\theta_{lj}}\right)\exp\left(-\frac{\sqrt{3}|x_{j}-x_{j}^{\prime}|}{\theta_{lj}}\right)$,
$\xi_{i}$ and $\zeta_{ik}$ are provided as follows,

$\displaystyle\xi_{i}$
$\displaystyle=\exp\left(\frac{3\sigma^{*2}_{l-1}(\mathbf{x})+2\sqrt{3}\theta_{ly}(y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x}))}{2\theta_{ly}^{2}}\right)\times\left[E_{1}^{\prime}\Lambda_{11}\Psi\left(\frac{\mu^{*}_{l-1}(\mathbf{x})-y_{i}^{[l-1]}-\frac{\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}}{\sigma^{*}_{l-1}(\mathbf{x})}\right)\right.$
$\displaystyle+\left.E_{1}^{\prime}\Lambda_{12}\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})+\frac{\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\right]$
$\displaystyle+\exp\left(\frac{3\sigma^{*2}_{l-1}(\mathbf{x})-2\sqrt{3}\theta_{ly}(y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x}))}{2\theta_{ly}^{2}}\right)\times\left[E_{2}^{\prime}\Lambda_{21}\Psi\left(\frac{-\mu^{*}_{l-1}(\mathbf{x})+y_{i}^{[l-1]}-\frac{\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}}{\sigma^{*}_{l-1}(\mathbf{x})}\right)\right.$
$\displaystyle+\left.E_{2}^{\prime}\Lambda_{12}\cdot\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})-\frac{\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\right],$
$\displaystyle\zeta_{ik}$
$\displaystyle=\exp\left\\{\frac{6\sigma^{*2}_{l-1}(\mathbf{x})+\sqrt{3}\theta_{ly}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}-2\mu^{*}_{l-1}(\mathbf{x})\right)}{\theta_{ly}^{2}}\right\\}$
$\displaystyle\times\left[E_{3}^{\prime}\Lambda_{31}\Psi\left\\{\frac{\left(\mu^{*}_{l-1}(\mathbf{x})-y_{k}^{[l-1]}-2\sqrt{3}\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)}{\sigma^{*}_{l-1}(\mathbf{x})}\right.\right\\}$
$\displaystyle+E_{3}^{\prime}\Lambda_{32}\left.\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(y_{k}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})+2\sqrt{3}\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\right]+\exp\left\\{-\frac{\sqrt{3}\left(y_{k}^{[l-1]}-y_{i}^{[l-1]}\right)}{\theta_{ly}}\right\\}$
$\displaystyle\times\left[E_{4}^{\prime}\Lambda_{41}\left(\Psi\left\\{\frac{y_{k}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})}{\sigma^{*}_{l-1}(\mathbf{x})}\right\\}-\Psi\left\\{\frac{y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})}{\sigma^{*}_{l-1}(\mathbf{x})}\right\\}\right)\right.$
$\displaystyle+E_{4}^{\prime}\Lambda_{42}\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\left.-E_{4}^{\prime}\Lambda_{43}\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(y_{k}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\right]$
$\displaystyle+\exp\left\\{\frac{6\sigma^{*2}_{l-1}(\mathbf{x})-\sqrt{3}\theta_{ly}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}-2\mu^{*}_{l-1}(\mathbf{x})\right)}{\theta_{ly}^{2}}\right\\}$
$\displaystyle\times\left[E_{5}^{\prime}\Lambda_{51}\Psi\left\\{\frac{\left(-\mu^{*}_{l-1}(\mathbf{x})+y_{k}^{[l-1]}-2\sqrt{3}\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)}{\sigma^{*}_{l-1}(\mathbf{x})}\right\\}\right.$

$\displaystyle+E_{5}^{\prime}\Lambda_{52}\left.\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})-2\sqrt{3}\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\right],$
$\displaystyle\Lambda_{11}$ $\displaystyle=\begin{pmatrix}1\\\
\mu^{*}_{l-1}(\mathbf{x})-\frac{\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\end{pmatrix},\Lambda_{12}=\begin{pmatrix}0\\\
1\end{pmatrix},\Lambda_{21}=\begin{pmatrix}1\\\
-\mu^{*}_{l-1}(\mathbf{x})-\frac{\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\end{pmatrix},$
$\displaystyle\Lambda_{31}$ $\displaystyle=\begin{pmatrix}1\\\
\mu^{*}_{l-1}(\mathbf{x})-\frac{2\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\\\
\left(\mu^{*}_{l-1}(\mathbf{x})-\frac{2\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}+\sigma^{*2}_{l-1}(\mathbf{x})\end{pmatrix},\Lambda_{32}=\begin{pmatrix}0\\\
1\\\
\mu^{*}_{l-1}(\mathbf{x})-\frac{2\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}+y_{k}^{[l-1]}\end{pmatrix},$
$\displaystyle\Lambda_{41}$ $\displaystyle=\begin{pmatrix}1\\\
\mu^{*}_{l-1}(\mathbf{x})\\\
\left(\mu^{*}_{l-1}(\mathbf{x})\right)^{2}+\sigma^{*2}_{l-1}(\mathbf{x})\end{pmatrix},\Lambda_{42}=\begin{pmatrix}0\\\
1\\\ \mu^{*}_{l-1}(\mathbf{x})+y_{i}^{[l-1]}\end{pmatrix},$
$\displaystyle\Lambda_{43}$ $\displaystyle=\begin{pmatrix}0\\\ 1\\\
\mu^{*}_{l-1}(\mathbf{x})+y_{k}^{[l-1]}\end{pmatrix},\Lambda_{51}=\begin{pmatrix}1\\\
-\mu^{*}_{l-1}(\mathbf{x})-\frac{2\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\\\
\left(-\mu^{*}_{l-1}(\mathbf{x})-\frac{2\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}+\sigma^{*2}_{l-1}(\mathbf{x})\end{pmatrix},$
$\displaystyle\Lambda_{52}$ $\displaystyle=\begin{pmatrix}0\\\ 1\\\
-\mu^{*}_{l-1}(\mathbf{x})-\frac{2\sqrt{3}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}-y_{i}^{[l-1]}\end{pmatrix},E_{1}=\frac{1}{\theta_{ly}}\begin{pmatrix}\theta_{ly}-\sqrt{3}y_{i}^{[l-1]}\\\
\sqrt{3}\end{pmatrix},$ $\displaystyle E_{2}$
$\displaystyle=\frac{1}{\theta_{ly}}\begin{pmatrix}\theta_{ly}+\sqrt{3}y_{i}^{[l-1]}\\\
\sqrt{3}\end{pmatrix},E_{3}=\frac{1}{\theta_{ly}^{2}}\begin{pmatrix}\theta_{ly}^{2}+3y_{i}^{[l-1]}y_{k}^{[l-1]}-\sqrt{3}\theta_{ly}(y_{i}^{[l-1]}+y_{k}^{[l-1]})\\\
2\sqrt{3}\theta_{ly}-3(y_{i}^{[l-1]}+y_{k}^{[l-1]})\\\ 3\end{pmatrix},$
$\displaystyle E_{4}$
$\displaystyle=\frac{1}{\theta_{ly}^{2}}\begin{pmatrix}\theta_{ly}^{2}-3y_{i}^{[l-1]}y_{k}^{[l-1]}+\sqrt{3}\theta_{ly}(y_{k}^{[l-1]}-y_{i}^{[l-1]})\\\
3(y_{i}^{[l-1]}+y_{k}^{[l-1]})\\\ -3\end{pmatrix},$ $\displaystyle E_{5}$
$\displaystyle=\frac{1}{\theta_{ly}^{2}}\begin{pmatrix}\theta_{ly}^{2}+3y_{i}^{[l-1]}y_{k}^{[l-1]}+\sqrt{3}\theta_{ly}(y_{i}^{[l-1]}+y_{k}^{[l-1]})\\\
2\sqrt{3}\theta_{ly}+3(y_{i}^{[l-1]}+y_{k}^{[l-1]})\\\ 3\end{pmatrix},$

for $y_{i}^{[l-1]}\leq y_{k}^{[l-1]}$. If $y_{i}^{[l-1]}>y_{k}^{[l-1]}$,
interchange $y_{i}^{[l-1]}$ and $y_{k}^{[l-1]}$. $\Psi$ is the cumulative
distribution function of a standard normal distribution.

###  S2.2 Matérn-2.5 kernel

For Matérn-2.5 kernel
$\Phi_{l}(\mathbf{x},\mathbf{x}^{\prime})=\prod_{j=1}^{d}\left(1+\frac{\sqrt{5}|x_{j}-x_{j}^{\prime}|}{\theta_{lj}}+\frac{5(x_{j}-x_{j}^{\prime})^{2}}{3\theta_{lj}^{2}}\right)\exp\left(-\frac{\sqrt{5}|x_{j}-x_{j}^{\prime}|}{\theta_{lj}}\right)$,
$\xi_{i}$ and $\zeta_{ik}$ are provided as follows,

$\displaystyle\xi_{i}$
$\displaystyle=\exp\left(\frac{5\sigma^{*2}_{l-1}(\mathbf{x})+2\sqrt{5}\theta_{ly}(y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x}))}{2\theta_{ly}^{2}}\right)\times\left[E_{1}^{\prime}\Lambda_{11}\Psi\left(\frac{\mu^{*}_{l-1}(\mathbf{x})-y_{i}^{[l-1]}-\frac{\sqrt{5}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}}{\sigma^{*}_{l-1}(\mathbf{x})}\right)\right.$
$\displaystyle+\left.E_{1}^{\prime}\Lambda_{12}\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(\mu^{*}_{l-1}(\mathbf{x})-y_{i}^{[l-1]}-\frac{\sqrt{5}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\right]$
$\displaystyle+\exp\left(\frac{5\sigma^{*2}_{l-1}(\mathbf{x})-2\sqrt{5}\theta_{ly}(y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x}))}{2\theta_{ly}^{2}}\right)\times\left[E_{2}^{\prime}\Lambda_{21}\Psi\left(\frac{-\mu^{*}_{l-1}(\mathbf{x})+y_{i}^{[l-1]}-\frac{\sqrt{5}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}}{\sigma^{*}_{l-1}(\mathbf{x})}\right)\right.$
$\displaystyle+\left.E_{2}^{\prime}\Lambda_{22}\cdot\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(-\mu^{*}_{l-1}(\mathbf{x})+y_{i}^{[l-1]}-\frac{\sqrt{5}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\right],$
$\displaystyle\zeta_{ik}$
$\displaystyle=\exp\left\\{\frac{10\sigma^{*2}_{l-1}(\mathbf{x})+\sqrt{5}\theta_{ly}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}-2\mu^{*}_{l-1}(\mathbf{x})\right)}{\theta_{ly}^{2}}\right\\}$
$\displaystyle\times\left[E_{3}^{\prime}\Lambda_{31}\Psi\left\\{\frac{\left(\mu^{*}_{l-1}(\mathbf{x})-y_{k}^{[l-1]}-2\sqrt{5}\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)}{\sigma^{*}_{l-1}(\mathbf{x})}\right.\right\\}$
$\displaystyle+E_{3}^{\prime}\Lambda_{32}\left.\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(\mu^{*}_{l-1}(\mathbf{x})-y_{k}^{[l-1]}-2\sqrt{5}\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\right]+\exp\left\\{-\frac{\sqrt{5}\left(y_{k}^{[l-1]}-y_{i}^{[l-1]}\right)}{\theta_{ly}}\right\\}$
$\displaystyle\times\left[E_{4}^{\prime}\Lambda_{41}\left(\Psi\left\\{\frac{y_{k}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})}{\sigma^{*}_{l-1}(\mathbf{x})}\right\\}-\Psi\left\\{\frac{y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})}{\sigma^{*}_{l-1}(\mathbf{x})}\right\\}\right)\right.$
$\displaystyle+E_{4}^{\prime}\Lambda_{42}\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(y_{i}^{[l-1]}-\mu^{*}_{l-1}(\mathbf{x})\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\left.-E_{4}^{\prime}\Lambda_{43}\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(y_{k}^{[l-1]}-\mu^{*}_{t-1}(\mathbf{x})\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\right]$

$\displaystyle+\exp\left\\{\frac{10\sigma^{*2}_{l-1}(\mathbf{x})-\sqrt{5}\theta_{ly}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}-2\mu^{*}_{l-1}(\mathbf{x})\right)}{\theta_{ly}^{2}}\right\\}$
$\displaystyle\times\left[E_{5}^{\prime}\Lambda_{51}\Psi\left\\{\frac{\left(-\mu^{*}_{l-1}(\mathbf{x})+y_{k}^{[l-1]}-2\sqrt{5}\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)}{\sigma^{*}_{l-1}(\mathbf{x})}\right\\}\right.$
$\displaystyle+E_{5}^{\prime}\Lambda_{52}\left.\frac{\sigma^{*}_{l-1}(\mathbf{x})}{\sqrt{2\pi}}\exp\left(-\frac{\left(-\mu^{*}_{l-1}(\mathbf{x})+y_{i}^{[l-1]}-2\sqrt{5}\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}}{2\sigma^{*2}_{l-1}(\mathbf{x})}\right)\right],$
$\displaystyle\Lambda_{11}$ $\displaystyle=\begin{pmatrix}1\\\
\mu^{*}_{t-1}(\mathbf{x})-\frac{\sqrt{5}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\\\
\left(\mu^{*}_{l-1}(\mathbf{x})-\frac{\sqrt{5}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}+\sigma^{*2}_{l-1}(\mathbf{x})\end{pmatrix},\Lambda_{12}=\begin{pmatrix}0\\\
1\\\
\mu^{*}_{l-1}(\mathbf{x})-\frac{\sqrt{5}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}+y_{i}^{[l-1]}\end{pmatrix},$
$\displaystyle\Lambda_{21}$ $\displaystyle=\begin{pmatrix}1\\\
-\mu^{*}_{l-1}(\mathbf{x})-\frac{\sqrt{5}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\\\
\left(\mu^{*}_{l-1}(\mathbf{x})+\frac{\sqrt{5}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}\right)^{2}+y_{i}^{[l-1]}\end{pmatrix},\Lambda_{22}=\begin{pmatrix}0\\\
1\\\
-\mu^{*}_{l-1}(\mathbf{x})-\frac{\sqrt{5}\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}}-y_{i}^{[l-1]}\end{pmatrix},$
$\displaystyle\Lambda_{31}$ $\displaystyle=\begin{pmatrix}1\\\ \mu_{c}\\\
\mu_{c}^{2}+\sigma^{*2}_{l-1}(\mathbf{x})\\\
\mu_{c}\left(\mu_{c}^{2}+3\sigma^{*2}_{l-1}(\mathbf{x})\right)\\\
\mu_{c}^{4}+6\mu_{c}^{2}\sigma^{*2}_{l-1}(\mathbf{x})+3\sigma^{*4}_{l-1}(\mathbf{x})\end{pmatrix},\mu_{c}=\mu^{*}_{l-1}(\mathbf{x})-2\sqrt{5}\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}},$
$\displaystyle\Lambda_{32}$ $\displaystyle=\begin{pmatrix}0\\\ 1\\\
\mu_{c}+y_{k}^{[l-1]}\\\
\mu_{c}^{2}+2\sigma^{*2}_{l-1}(\mathbf{x})+\left(y_{k}^{[l-1]}\right)^{2}+\mu_{c}y_{k}^{[l-1]}\\\
\mu_{c}^{3}+\left(y_{k}^{[l-1]}\right)^{3}+y_{k}^{[l-1]}\mu_{c}\left(\mu_{c}+y_{k}^{[l-1]}\right)+\sigma^{*2}_{l-1}(\mathbf{x})\left(5\mu_{c}+3y_{k}^{[l-1]}\right)\end{pmatrix},$

$\displaystyle\Lambda_{41}$ $\displaystyle=\begin{pmatrix}1\\\ \mu^{*}\\\
\mu^{*2}+\sigma^{*2}_{l-1}(\mathbf{x})\\\
\mu^{*}\left(\mu^{*2}+3\sigma^{*2}_{l-1}(\mathbf{x})\right)\\\
\mu^{*4}+6\mu_{c}^{2}\sigma^{*2}_{l-1}(\mathbf{x})+3\sigma^{*4}_{l-1}(\mathbf{x})\end{pmatrix},\mu^{*}=\mu^{*}_{l-1}(\mathbf{x}),$
$\displaystyle\Lambda_{42}$ $\displaystyle=\begin{pmatrix}0\\\ 1\\\
\mu^{*}+y_{i}^{[l-1]}\\\
\mu^{*2}+2\sigma^{*2}_{l-1}(\mathbf{x})+\left(y_{i}^{[l-1]}\right)^{2}+\mu^{*}y_{i}^{[l-1]}\\\
\mu^{*3}+\left(y_{i}^{[l-1]}\right)^{3}+y_{i}^{[l-1]}\mu^{*}\left(\mu^{*}+y_{i}^{[l-1]}\right)+\sigma^{*2}_{l-1}(\mathbf{x})\left(5\mu^{*}+3y_{i}^{[l-1]}\right)\end{pmatrix},$
$\displaystyle\Lambda_{43}$ $\displaystyle=\begin{pmatrix}0\\\ 1\\\
\mu^{*}+y_{k}^{[l-1]}\\\
\mu^{*2}+2\sigma^{*2}_{l-1}(\mathbf{x})+\left(y_{k}^{[l-1]}\right)^{2}+\mu^{*}y_{k}^{[l-1]}\\\
\mu^{*3}+\left(y_{k}^{[l-1]}\right)^{3}+y_{k}^{[l-1]}\mu^{*}\left(\mu^{*}+y_{k}^{[l-1]}\right)+\sigma^{*2}_{l-1}(\mathbf{x})\left(5\mu^{*}+3y_{k}^{[l-1]}\right)\end{pmatrix},$
$\displaystyle\Lambda_{51}$ $\displaystyle=\begin{pmatrix}1\\\ -\mu_{d}\\\
\mu_{d}^{2}+\sigma^{*2}_{l-1}(\mathbf{x})\\\
-\mu_{d}\left(\mu_{d}^{2}+3\sigma^{*2}_{l-1}(\mathbf{x})\right)\\\
\mu_{d}^{4}+6\mu_{d}^{2}\sigma^{*2}_{l-1}(\mathbf{x})+3\sigma^{*4}_{l-1}(\mathbf{x})\end{pmatrix},\mu_{d}=\mu^{*}_{l-1}(\mathbf{x})+2\sqrt{5}\frac{\sigma^{*2}_{l-1}(\mathbf{x})}{\theta_{ly}},$
$\displaystyle\Lambda_{52}$ $\displaystyle=\begin{pmatrix}0\\\ 1\\\
-\mu_{d}-y_{i}^{[l-1]}\\\
\mu_{d}^{2}+2\sigma^{*2}_{l-1}(\mathbf{x})+\left(y_{i}^{[l-1]}\right)^{2}+\mu_{d}y_{i}^{[l-1]}\\\
-\mu_{d}^{3}-\left(y_{i}^{[l-1]}\right)^{3}-y_{i}^{[l-1]}\mu_{d}\left(\mu_{d}+y_{i}^{[l-1]}\right)-\sigma^{*2}_{l-1}(\mathbf{x})\left(5\mu_{d}+3y_{i}^{[l-1]}\right)\end{pmatrix},$

$\displaystyle E_{1}$
$\displaystyle=\frac{1}{3\theta_{ly}^{2}}\begin{pmatrix}3\theta_{ly}^{2}-3\sqrt{5}\theta_{ly}y_{i}^{[l-1]}+5\left(y_{i}^{[l-1]}\right)^{2}\\\
3\sqrt{5}\theta_{ly}-10y_{i}^{[l-1]}\\\
5\end{pmatrix},E_{2}=\frac{1}{3\theta_{ly}^{2}}\begin{pmatrix}3\theta_{ly}^{2}+3\sqrt{5}\theta_{ly}y_{i}^{[l-1]}+5\left(y_{i}^{[l-1]}\right)^{2}\\\
3\sqrt{5}\theta_{ly}+10y_{i}^{[l-1]}\\\ 5\end{pmatrix},$ $\displaystyle E_{3}$
$\displaystyle=\frac{1}{9\theta_{ly}^{4}}\begin{pmatrix}E_{31}&E_{32}&E_{33}&E_{34}&E_{35}\end{pmatrix}^{\top},$
$\displaystyle E_{31}$
$\displaystyle=9\theta_{ly}^{4}+25\left(y_{i}^{[l-1]}\right)^{2}\left(y_{k}^{[l-1]}\right)^{2}-3\sqrt{5}\theta_{ly}\left(3\theta_{ly}^{2}+5y_{i}^{[l-1]}y_{k}^{[l-1]}\right)\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right)$
$\displaystyle+15\theta_{ly}^{2}\left(\left(y_{i}^{[l-1]}\right)^{2}+\left(y_{k}^{[l-1]}\right)^{2}+3y_{i}^{[l-1]}y_{k}^{[l-1]}\right),$
$\displaystyle E_{32}$
$\displaystyle=18\sqrt{5}\theta_{ly}^{3}+15\sqrt{5}\theta_{ly}\left(\left(y_{i}^{[l-1]}\right)^{2}+\left(y_{k}^{[l-1]}\right)^{2}\right)-75\theta_{ly}^{2}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right)$
$\displaystyle-50y_{i}^{[l-1]}y_{k}^{[l-1]}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right)+60\sqrt{5}\theta_{ly}y_{i}^{[l-1]}y_{k}^{[l-1]},$
$\displaystyle E_{33}$
$\displaystyle=5\left\\{5\left(y_{i}^{[l-1]}\right)^{2}+5\left(y_{k}^{[l-1]}\right)^{2}+15\theta_{ly}^{2}-9\sqrt{5}\theta_{ly}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right)+20\left(y_{i}^{[l-1]}y_{k}^{[l-1]}\right)\right\\}$
$\displaystyle E_{34}$
$\displaystyle=10\left(3\sqrt{5}\theta_{ly}-5y_{i}^{[l-1]}-5y_{k}^{[l-1]}\right),E_{35}=25,$
$\displaystyle E_{4}$
$\displaystyle=\frac{1}{9\theta_{ly}^{4}}\begin{pmatrix}E_{41}&E_{42}&E_{43}&E_{44}&E_{45}\end{pmatrix}^{\top},$
$\displaystyle E_{41}$
$\displaystyle=9\theta_{ly}^{4}+25\left(y_{i}^{[l-1]}\right)^{2}\left(y_{k}^{[l-1]}\right)^{2}+3\sqrt{5}\theta_{ly}\left(3\theta_{ly}^{2}-5y_{i}^{[l-1]}y_{k}^{[l-1]}\right)\left(y_{i}^{[l-1]}-y_{k}^{[l-1]}\right)$
$\displaystyle+15\theta_{ly}^{2}\left(\left(y_{i}^{[l-1]}\right)^{2}+\left(y_{k}^{[l-1]}\right)^{2}-3y_{i}^{[l-1]}y_{k}^{[l-1]}\right),$
$\displaystyle E_{42}$
$\displaystyle=5\left\\{3\sqrt{5}\theta_{ly}\left(\left(y_{k}^{[l-1]}\right)^{2}-\left(y_{i}^{[l-1]}\right)^{2}\right)+3\theta_{ly}^{2}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right)-10y_{i}^{[l-1]}y_{k}^{[l-1]}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right)\right\\},$
$\displaystyle E_{43}$
$\displaystyle=5\left\\{5\left(y_{i}^{[l-1]}\right)^{2}+5\left(y_{k}^{[l-1]}\right)^{2}-3\theta_{ly}^{2}-3\sqrt{5}\theta_{ly}\left(y_{k}^{[l-1]}-y_{i}^{[l-1]}\right)+20\left(y_{i}^{[l-1]}y_{k}^{[l-1]}\right)\right\\}$
$\displaystyle E_{44}$
$\displaystyle=-50\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right),E_{45}=25,$
$\displaystyle E_{5}$
$\displaystyle=\frac{1}{9\theta_{ly}^{4}}\begin{pmatrix}E_{51}&E_{52}&E_{53}&E_{54}&E_{55}\end{pmatrix}^{\top},$
$\displaystyle E_{51}$
$\displaystyle=9\theta_{ly}^{4}+25\left(y_{i}^{[l-1]}\right)^{2}\left(y_{k}^{[l-1]}\right)^{2}+3\sqrt{5}\theta_{ly}\left(3\theta_{ly}^{2}+5y_{i}^{[l-1]}y_{k}^{[l-1]}\right)\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right)$
$\displaystyle+15\theta_{ly}^{2}\left(\left(y_{i}^{[l-1]}\right)^{2}+\left(y_{k}^{[l-1]}\right)^{2}+3y_{i}^{[l-1]}y_{k}^{[l-1]}\right),$
$\displaystyle E_{52}$
$\displaystyle=18\sqrt{5}\theta_{ly}^{3}+15\sqrt{5}\theta_{ly}\left(\left(y_{i}^{[l-1]}\right)^{2}+\left(y_{k}^{[l-1]}\right)^{2}\right)+75\theta_{ly}^{2}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right)$
$\displaystyle+50y_{i}^{[l-1]}y_{k}^{[l-1]}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right)+60\sqrt{5}\theta_{ly}y_{i}^{[l-1]}y_{k}^{[l-1]},$

$\displaystyle E_{53}$
$\displaystyle=5\left\\{5\left(y_{i}^{[l-1]}\right)^{2}+5\left(y_{k}^{[l-1]}\right)^{2}+15\theta_{ly}^{2}+9\sqrt{5}\theta_{ly}\left(y_{i}^{[l-1]}+y_{k}^{[l-1]}\right)+20\left(y_{i}^{[l-1]}y_{k}^{[l-1]}\right)\right\\}$
$\displaystyle E_{54}$
$\displaystyle=10\left(3\sqrt{5}\theta_{ly}+5y_{i}^{[l-1]}+5y_{k}^{[l-1]}\right),E_{55}=25,$

for $y_{i}^{[l-1]}\leq y_{k}^{[l-1]}$. If $y_{i}^{[l-1]}>y_{k}^{[l-1]}$,
interchange $y_{i}^{[l-1]}$ and $y_{k}^{[l-1]}$.