# An adaptive strategy for sequential designs of multilevel computer
experiments

Ayao Ehara
Department of Statistical Science
University College London
London, UK
&Serge Guillas
Department of Statistical Science
University College London
London, UK
Corresponding author<EMAIL_ADDRESS>
###### Abstract

Investigating uncertainties in computer simulations can be prohibitive in
terms of computational costs, since the simulator needs to be run over a large
number of input values. Building an emulator, i.e. a statistical surrogate
model of the simulator, using a small design of experiments, greatly
alleviates the computational burden to carry out such investigations.
Nevertheless, this can still be above the computational budget for many
studies. Two major approaches have been used are to reduce the budget needed
to build the emulator: efficient design of experiments, such as sequential
designs, and combining training data of different degrees of sophistication in
a so-called multi-fidelity method, or multilevel in case these fidelities are
ordered typically for increasing resolutions. We present here a novel method
that combines both approaches, the multilevel adaptive sequential design of
computer experiments (MLASCE) in the framework of Gaussian process (GP)
emulators. We make use of reproducing kernel Hilbert spaces as a new tool for
our GP approximations of the increments across levels. This dual strategy
allows us to allocate efficiently limited computational resource over
simulations of different levels of fidelity and build the GP emulator. The
allocation of computational resources is shown to be the solution of a simple
optimization problem in a special case where we theoretically prove the
validity of our approach. Our proposed method is compared with other existing
models of multi-fidelity Gaussian process emulation. Gains of orders of
magnitudes in accuracy for medium-size computing budgets are demonstrated in
numerical examples.

_K_ eywords uncertainty quantification $\cdot$ surrogate models $\cdot$
Gaussian process $\cdot$ multi-fidelity $\cdot$ sequential design $\cdot$
reproducing kernel Hilbert spaces (RKHS)

##  1 Introduction

Computer experiments have been used to study real physical systems in many
areas of scientific research. The parameters which determine the system are
typically treated as input variables of the quantity of interest. Since
running a complex model is computationally expensive, it is hard to understand
fully this input-output structure without relying on another method. In [27],
the fundamental framework of design and analysis of computer experiments
(DACE) was developed. Their approach is to model the deterministic output from
computer simulation as the realization of a stochastic process, thereby
providing a statistical basis for designing experiments for efficient
prediction. Statistical emulators typically based on Gaussian processes are
constructed subsequently for the prediction of input-output structure in their
framework and it allows us to be free from relying only on running a computer
simulation. Gaussian process regression [37, 30] (also called Kriging) are a
commonly used class of surrogate models as a statistical emulator, which
assume that prior beliefs about the input-output structure can be modelled by
a Gaussian process [6, 15].

A computer simulation often can be run at different levels of accuracy, with
versions ranging from the most sophisticated simulator to the most basic. The
sophisticated versions are accurate but so expensive that only a limited
numbers of runs are obtained. In contrast, the low-level simulations return
relatively inaccurate results but are cheaper to implement and we can use a
large number of samples. Then, it is a natural question how the data from
cheap computer simulations can effectively complement the small size of
samples from expensive simulations. In other words, how can we use the data
from several simulators with different levels of accuracy in the most
effective way to construct a statistical emulator? Multi-fidelity statistical
emulators have been employed to address this type of question. The first
multi-fidelity statistical emulator proposed in [7] was based on a linear
regression formulation. Then, this model was improved in [8] by employing a
Bayes linear formulation. The drawback is a lack of accuracy since they are
based on a linear regression formulation. Following the fundamental framework
of [27], multi-fidelity Gaussian process regression, or co-kriging, was
considered in a seminal paper in [14] as a more general methodology of
interpolation, also followed by [11, 26, 25]. The formulation was articulated
[17, 19] and another method for obtaining training data as a strategy of
sequential experimental design was developed in [18]. Recent papers such as
[24, 40, 41] took advantage of this type of method and applied to settings in
physics. The benefit of these models is that the correlation between the data
of high and low fidelity is introduced so that it allows us to efficiently
construct an accurate emulator by performing relatively few evaluations of an
expensive high-fidelity model and more evaluations of a cheaper one. However,
none of these approaches addresses the relationships between the precision of
an emulator and the number of runs from each simulator, which is crucial to
obtain the most accurate surrogate model under limited computational budget.
In other words, the bias between the true response and a statistical emulator
is never accounted hence the emulator would have inherently significant bias.
In addition, choosing the location of data, or design of experiments, directly
affects the quality of an emulator. Most of the related works use space-
filling or Latin hypercube design, which would be inefficient especially in
the case of high dimensions. Moreover, they assume nested design across every
level of fidelity, which means that the simulated points for the more
sophisticated computer simulations are included in the design points for less
accurate ones; this may cause insufficient samplings from the locations where
higher fidelity simulations behaves differently from the ones of low fidelity.

The contributions of this research are two-fold. First our approach is the
first multi-fidelity emulation method that takes advantage of sequential
designs, which have been shown to be superior to one-shot designs for computer
experiments [1]. Secondly, compared to previous multi-fidelity methods, we
allow our choice of design points across fidelities to be fully flexible,
based upon a criterion that seeks to select the fidelity that reduces the most
the difference of mean predicted increments across levels. This dual strategy
provides large gains compared to state-of-the-art multi-fidelity experimental
designs. Defining the accuracy of an emulator and evaluating the effect of
additional expensive computer simulations are necessary steps in optimizing
the usage of computational resources. This is made possible by decomposing the
computer simulations into a sum of incremental functions and assuming a
reproducing kernel Hilbert space (RKHS) as a function space to which each
incremental function belong. Evaluating the accuracy of an approximation in
the RKHS guarantees that our emulator gets closer to the underlining function
as the number of runs increases. Regarding this type of asymptotic theory, a
profound literature exists in scattered data approximation and related fields,
which is often used in Gaussian process approximation. The exhaustive study
was given in [35]. In the context of practical implementation of Gaussian
process regressions in recent papers, the effects of misspecified parameters
in covariance functions were investigated in [33, 34] while convergence rates
with statistically estimated parameters was studied in [32, 39]. The
convergence rates by these methods are based on geometric features of the
domain and design. Especially in the case of high dimension and a tensor grid
design, these convergence rates are quite slow. Therefore, we rely on an
empirical but usually more efficient design strategy, Mutual Information
Computer Experiment (MICE) [1], to obtain faster convergence in a practical
way.

This paper is organized as follows. In section 2 we present a review of
Gaussian process regression for emulation and existing multi-fidelity methods.
In section our strategy of experimental design and the choice of covariance
structure in relation with reproducing kernel Hilbert spaces. It also provides
an overview of multilevel Monte Carlo and an error bound of the kriging
predictor for our analysis as the background theory. Section 3 is the core
part of this research. We show the novel multilevel adaptive sequential design
of computer experiments (MLASCE) in the framework of multi-fidelity GPR and an
algorithm for constructing it. Section 4 details an an independent analysis of
optimal numbers of runs in the special case of tensor grid design. Numerical
examples are given in section 5 showing MLASCE performs better than other
existing multi-fidelity models.

##  2 Gaussian process (GP) emulation: background and multi-fidelity context

We assume that $\mathbb{D}\subseteq\mathbb{R}^{d}$ is the domain of input for
a positive integer $d$ and $y:\mathbb{D}\rightarrow\mathbb{R}$ denotes the
deterministic quantity of interest dependent on the input. Further, let
$\boldsymbol{y}(X_{N})=\left(y(x_{1}),y(x_{2})\ldots,y(x_{N})\right)^{\top}$
denote $N$ observations of $y$ collected at the inputs
$X_{N}=\left(x_{1},x_{2},\ldots,x_{N}\right)^{\top}$ where
$x_{i}\in\mathbb{D}\subseteq\mathbb{R}^{d}$, $y(x_{i})\in\mathbb{R}$. We aim
to predict $\boldsymbol{y}$ at any new input $\boldsymbol{x}\in\mathbb{D}$.
The GP regression (GPR) assumes that the observation vector
$\boldsymbol{y}(X_{N})$ is a realization of the $N$ -dimensional random vector
with multivariate Gaussian distribution

$\boldsymbol{Y}=\left(Y\left(\boldsymbol{x}_{1},\omega\right),Y\left(\boldsymbol{x}_{2},\omega\right),\ldots,Y\left(\boldsymbol{x}_{N},\omega\right)\right)^{\top}$

where $Y(\cdot,\cdot):\mathbb{D}\times\Omega\rightarrow\mathbb{R}$ is a GP
defined on the probability space $(\Omega,\mathcal{F},P)$. For brevity,
$Y(\boldsymbol{x},\cdot)$ is denoted by $Y(\boldsymbol{x})$. The GP $Y$ is
usually represented using GP notation as

$Y(\boldsymbol{x})\sim
GP\left(m(\boldsymbol{x}),K\left(\boldsymbol{x},\boldsymbol{x}^{\prime}\right)\right)$

where $m(\cdot):\mathbb{D}\rightarrow\mathbb{R}$ and
$K(\cdot,\cdot):\mathbb{D}\times\mathbb{D}\rightarrow\mathbb{R}$ are the mean
and covariance functions

$\displaystyle m(\boldsymbol{x})$
$\displaystyle=\mathbb{E}[Y(\boldsymbol{x})]$ $\displaystyle
K\left(\boldsymbol{x},\boldsymbol{x}^{\prime}\right)$
$\displaystyle=\operatorname{Cov}\bigl{(}Y(\boldsymbol{x}),Y\left(\boldsymbol{x}^{\prime}\right)\bigr{)}=\mathbb{E}\bigl{[}\bigl{(}Y(\boldsymbol{x})-m(\boldsymbol{x})\bigr{)}\bigl{(}Y\left(\boldsymbol{x}^{\prime}\right)-m\left(\boldsymbol{x}^{\prime}\right)\bigr{)}\bigr{]}$

The variance of $Y(\boldsymbol{x})$ is $K(\boldsymbol{x},\boldsymbol{x})$ and
its standard deviation is
$\sigma(\boldsymbol{x})=\sqrt{K(\boldsymbol{x},\boldsymbol{x})}.$ The
covariance matrix of the random vector $\boldsymbol{Y}$ is then given by

$K(X_{N},X_{N})=\left(\begin{array}[]{ccc}k\left(\boldsymbol{x}_{1},\boldsymbol{x}_{1}\right)&\cdots&k\left(\boldsymbol{x}_{1},\boldsymbol{x}_{N}\right)\\\
\vdots&\ddots&\vdots\\\
k\left(\boldsymbol{x}_{N},\boldsymbol{x}_{1}\right)&\cdots&k\left(\boldsymbol{x}_{N},\boldsymbol{x}_{N}\right)\end{array}\right)$

The hyperparameters in $m(\boldsymbol{x})$ and
$K\left(\boldsymbol{x},\boldsymbol{x}^{\prime}\right)$ are typically
identified by maximizing the $\log$ marginal likelihood of the observations
(see examples in Section 2.2 in [37]). The GPR prediction at $x$ consists of
the posterior normal distribution
$y\left(\boldsymbol{x}\right)\sim\mathcal{N}\left(m_{N}^{y}\left(x\right),\hat{s}^{2}_{N}\left(\boldsymbol{x}\right)\right),$
with posterior or predictive mean and variance given by

$\displaystyle m_{N}^{y}\left(\boldsymbol{x}\right)$
$\displaystyle=m\left(\boldsymbol{x}\right)+K\left(\boldsymbol{x},X_{N}\right)^{\top}K(X_{N},X_{N})^{-1}(\boldsymbol{y}(X_{N})-\boldsymbol{m}(X_{N}))$
$\displaystyle\hat{\mathrm{s}}^{2}_{N}\left(\boldsymbol{x}\right)$
$\displaystyle=K(\boldsymbol{x},\boldsymbol{x})-K\left(\boldsymbol{x},X_{N}\right)^{\top}K(X_{N},X_{N})^{-1}K\left(\boldsymbol{x},X_{N}\right)$

where
$\boldsymbol{m}(X_{N})=\left(m\left(x_{1}\right),\ldots,m\left(x_{N}\right)\right)^{\top}$
and $K\left(\boldsymbol{x},X_{N}\right)$ is the vector of covariances

$K\left(\boldsymbol{x},X_{N}\right)=\left(K(\boldsymbol{x},\boldsymbol{x}_{1}),K(\boldsymbol{x},\boldsymbol{x}_{2}),\cdots,K(\boldsymbol{x},\boldsymbol{x}_{N})\right)^{\top}.$

We now present a brief overview of the initial approach of [14]. A co-kriging
model is introduced based on a first-order auto-regressive relation between
model outputs of different levels of fidelity.

Suppose we have $L$-levels of variable-fidelity computer simulations
$\left(y_{l}(x)\right)_{l=1}^{L}$.
$\boldsymbol{y}_{l}(X_{l,N_{l}})=\bigl{(}y_{l}(x_{1,l}),\ldots,$
$y_{l}(x_{N_{l},l})\bigr{)}^{\top}$ denotes the observations at the inputs
$X_{l,N_{l}}=\bigl{(}x_{1,l},\ldots,x_{N_{l},l}\bigr{)}$ sorted by increasing
fidelity and modelled as observations of a Gaussian process
$\left(Y_{l}(\mathrm{x})\right)_{l=1}^{L}$ where
$x_{i,l}\in\mathbb{D}_{l}\subseteq\mathbb{R}^{d},y_{l}(x_{i,l})\in\mathbb{R}$.
Then, $y_{L}(x)$ denotes the output of the most accurate and expensive model
while $y_{1}(x)$ is the output of the cheapest and least accurate surrogate at
our disposal. The auto-regressive co-kriging scheme of [14] reads as

$\displaystyle Y_{l}(x)=\rho_{l-1}Y_{l-1}(x)+\delta_{l}(x),\quad l=2,\ldots,L$
(1)

where $\delta_{l}(x)$ is a GP independent of
$\left\\{Y_{l-1}(x),\ldots,Y_{1}(x)\right\\}$. For convention, $Y_{1}(x)$ is a
GP $\delta_{1}(x)$. $\rho$ is a scaling factor and a deterministic scalar
independent of $\mathbf{x}$. $\rho$ also quantifies the correlation between
the model outputs $\left(y_{l}(x),y_{l-1}(x)\right)$. The derivation of this
model is based on the Markov property

$\operatorname{Cov}\bigl{(}Y_{l}(\mathbf{x}),Y_{l-1}\left(\mathbf{x}^{\prime}\right)|Y_{l-1}(\mathbf{x})\bigr{)}=0,\quad\forall\mathbf{x}\neq\mathbf{x}^{\prime}$

which translates into assuming that given $Y_{l-1}(x),$ we can learn nothing
more about $Y_{l}(\mathrm{x})$ from any other model output
$Y_{l-1}\left(\mathrm{x}^{\prime}\right),$ for
$\mathrm{x}\neq\mathrm{x}^{\prime}$.

The extended formulation of the multi-fidelity model was presented in [19] and
its special form and a strategy of sequential design in this special case were
considered in [18]. The unique difference with the previous model is that
$Y_{l}(x)$ is a function of $\tilde{Y}_{l-1}(x)$, which is the Gaussian
process $Y_{l-1}(x)$ conditioned by the values
$\\{\boldsymbol{y}_{i,N_{i}}\\}_{i=1}^{l-1}$ at points in the experimental
design sets $\\{\mathbb{D}_{i}\\}_{i=1,\ldots,l-1}$. They also introduced the
dependency of $\rho_{l}(x)$ on input $x$. As in the previous model, the nested
property for the experimental design sets
$\mathbb{D}_{L}\subseteq\mathbb{D}_{L-1}\subseteq\cdots\subseteq\mathbb{D}_{1}$
is assumed.

They consider the following model for $l=2,\ldots,L:$

$\displaystyle\left\\{\begin{array}[]{l}Y_{l}(x)=\rho_{l-1}(x)\tilde{Y}_{l-1}(x)+\delta_{l}(x)\\\
\tilde{Y}_{l-1}(x)\perp\delta_{l}(x)\\\
\rho_{l-1}(x)=g_{l-1}^{\top{}}(x)\beta_{\rho_{l-1}}\end{array}\right.$ (5)

where $\delta_{l}(x)$ is a Gaussian process
$GP\left(f_{t}^{T}(x)\beta_{t},\sigma_{t}^{2}r_{t}\left(x,x^{\prime}\right)\right)$,
$\perp$ denotes the independence relationship.

The estimation of hyperparameters is implemented by Maximum likelihood method
or minimizing a loss function of a Leave-One-Out Cross-Validation procedure. A
nested space filling design for choosing $\mathbb{D}_{l}$ was considered in
[19] and in [18], assuming constant $\rho$ thus a special case of [19], they
proposed a strategy of sequential design, which combines both the error
evaluation provided by the GPR model and the observed errors of a leave-one-
out cross-validation procedure.

One of the drawbacks of these formulations is that the bias between the true
function and their emulators is not considered. Therefore, the above methods
would have no guarantee of accurate prediction in theory and inherently
incorporate the bias.

###  2.1 Sequential design

An experimental design is a specification of inputs in the domain of our
interest at which we wish to compute the outputs [28]. In our case,
experimental design is about obtaining the training data and specifying the
appropriate location. In the case of GPR, the posterior mean and covariance
functions are used for its prediction. Therefore, the locations of inputs and
corresponding outputs used as the training data as well as the its number of
data points does affect the quality of the prediction. In design of
experiments, it is customary to use space-filling designs, which spread
observations evenly throughout the input region [29, 28], such as uniform
designs, multilayer designs, maximin (Mm)- and minimax (mM)-distance designs,
and Latin hypercube designs (LHD). Space-filling designs treat all parts of
the design space as equally important hence it can be inefficient especially
if the domain of input is high dimensional. A variety of adaptive designs (an
algorithm is called adaptive if it updates its behavior to new data) have been
proposed which can take advantage of information collected during the
experimental design process [28], typically in the form of input-output data
from simulation runs. In the case of GPR, the hyperparamters in kernel
functions are updated recursively whenever new data becomes available. Some
classical adaptive design criteria are the maximum mean squared prediction
error (MMSPE), the integrated MSPE (IMSPE), and the entropy criterion [27].

We employ a useful method called MICE (mutual information computer experiment)
design as a design strategy in this research. Other approaches for sequential
design have been introduced previously such as The active learning MacKay
(ALM) [20] and Active learning Cohn (ALC) [5] algorithm but MICE outperforms
these techniques in general. For our purpose of selecting a very small number
of runs in the design for high fidelity simulators, it is essential to pick
these runs very carefully to optimize the design. MICE is making a bigger
difference with its competitors with smaller designs, and hence is
contributing heavily to the quality of our approach.

For computer experiment, the domain $\mathbb{D}$ is often discretized into
grid $\mathbb{D}_{G}$ with $N_{G}$ points in practice since the optimization
search of training data is a formidable task if $\mathbb{D}$ is not.
Consequently, we replace the search over $\mathbb{D}$ by a search over a set
of candidate points $X_{\text{cand}}\subseteq\mathbb{D}_{G}$. We consider
sequential designs as practical, computationally cheaper alternatives to “one
shot” designs, which determines all input points in one step, albeit often
suboptimal. The sequential design is defined as follows: Suppose that we have
an initial design
$\\{(x_{j},y(x_{j}))\\}_{j=1}^{n}=\\{X_{n},\boldsymbol{y}(X_{n})\\}$, then for
each $k=n,n+1,n+2,\ldots$ one collects an input-output pair
$(x_{k+1},y(x_{k+1}))$ by choosing the input values

$\displaystyle x_{k+1}=\mathop{\rm arg~{}max}\limits_{x\in
X_{\text{cand}}\setminus X_{k}}F_{k}(x).$

where $F_{k}(\cdot)$ is a design criterion to be maximized. The algorithm
iterates until a stopping criterion is met or the computational budget
allocated is exhausted. We denote $\hat{s}^{2}_{k}(x)$ as a predictive
variance of a Gaussian process conditional on $X_{k}$ with $k$ points with
estimated parameters in $K$. In our context, the parameters in the kernel of
GP need to be estimated. Sequential designs allow the estimates to be improved
sequentially with the additional new design points. This is especially
advantageous when some input variables are considerably more influential on
the output than others.

####  2.1.1 Mutual information criteria

Mutual information (MI) is a classical information theoretic measure. The goal
is to place a set of design $\mathcal{X}$ that will give us good predictions
at all unexplored location $\mathbb{D}\setminus\mathcal{X}$. Specifically, we
want to find

$\displaystyle\mathcal{X}^{*}=\mathop{\rm
arg~{}max}\limits_{\mathcal{X}\subseteq\mathbb{D}_{cand}}H(\mathbb{D}_{G}\setminus\mathcal{X})-H(\mathbb{D}_{G}\setminus\mathcal{X}|\mathcal{X})$
(6)

where $H(X)$ denotes the entropy of Gaussian process
$Y(X)=\left(Y(x_{1}),\ldots,Y(x_{n})\right)$ with the location
$X=\left(x_{1},\ldots,x_{n}\right)$ and $H(X|Z)$ is the entropy of a
conditional Gaussian process $Y(X|Z)$ conditioned by the realized values
$Y(Z)$ with the location $Z=\left(z_{1},\ldots,z_{m}\right)$. The set
$\mathcal{X}^{*}$ maximally reduces the entropy over the rest of space
$\mathbb{D}\setminus\mathcal{X}$. The entropy of a Gaussian process is

$H(X)={n\over
2}\log(2\pi)\operatorname{det}\bigl{(}\operatorname{Var}(Y(X))\bigr{)}$

where $\operatorname{det}\left(\cdot\right)$ is the determinant of a matrix
and $\operatorname{Var}(Y(X))$ means the covariance matrix of $Y(X)$. Though
this optimization problem is NP-hard, [16] presented an alternative algorithm,
known as MI algorithm, which avoids directly maximizing the difference of
entropies.

####  2.1.2 Mutual information for computer experiments

MI reduces uncertainties over the domain with relatively cheap computational
costs. However, MI lacks adaptivity, that is, the parameters in the covariance
are assumed to be known. Moreover, computational instability is inherent due
to the inversion of the covariance matrix in the denominator, in Step 2 of MI
algorithm. [1] developed a modified approach known as Mutual information for
computer experiments (MICE). They introduced a different nugget parameter in
diagonal elements of the covariance matrix in the denominator to stabilize the
inversion of a matrix. In addition, they added adaptivity in the sense that
$\mathbb{D}_{G}$ is sampled in the design space $\mathbb{D}$ instead of
keeping it fixed and the parameters in the covariance matrix are estimated at
each step of sampling. The outlined algorithm is as follows [1]. The
covariance $K(\cdot,\cdot)$ is replaced with $K_{\theta}(\cdot,\cdot)$ to
emphasize the hyperparameters $\theta$ in $K$.

Require

: Function $y(x)$, GP emulator $(m(\cdot),K_{\theta}(\cdot,\cdot))$, nugget
parameters $\tau^{2}$ and $\tau^{2}_{s}$, domain of input $\mathbb{D}$,
initial data $\\{X_{k},\boldsymbol{y}(X_{k})\\}$, discrete set size $N_{G}$,
candidate set size $N_{\text{cand}}$, desired size $N$.

Step 1

Implement maximum likelihood estimation to obtain estimates $\hat{\theta}_{k}$
of $\theta$ based on $\\{X_{k},\boldsymbol{y}(X_{k})\\}$ in
$K_{\theta}(\cdot,\cdot)$.

Step 2

Fit GP emulator to data $\\{X_{k},\boldsymbol{y}(X_{k})\\}$.

Step 3

Generate a discrete set $\mathbb{D}_{G}$ of size $N_{G}$, and choose a
candidate set $X_{\text{cand}}\subseteq\mathbb{D}_{G}$.

Step 4

Solve $x_{k+1}=\mathop{\rm arg~{}max}\limits_{x\in
X_{\text{cand}}}\hat{s}^{2}_{k}(x,\hat{\theta}_{k},\tau^{2})/\hat{s}^{2}_{G\setminus(k\cup
x)}(x,\hat{\theta}_{k},\max({\tau^{2},\tau^{2}_{s}}))$

Step 5

Evaluate $y_{k+1}=y(x_{k+1})$, and let $X_{k+1}=X_{k}\cup x_{K+1}$ and
$\mathbb{y}(X_{k+1})=\mathbb{y}(X_{k})\cup y_{k+1}$.

Step 6

If $k+1=N$, then stop; otherwise let $k=k+1$, and go to Step 1.

Output

: The collection of inputs and outputs $\\{X_{N},\boldsymbol{y}(X_{N})\\}$ of
size $N$.

In step $2$, maximum likelihood method is implemented to get the estimation of
the hyperparameters in the covariance. For example, if the covariance kernel
is Gaussian, the shape and scale paramters are estimated while all of the
hyperparameters except the smoothness parameter are estimated in the case of
Matérn kernel, which we will discuss later. In Step $3$, they suggest that
$\mathbb{D}_{G}$ be sampled in the design space $\mathbb{D},$ instead of
keeping $\mathbb{D}_{G}$ fixed throughout. In Step 4 the MICE criterion is
evaluated for all $x\in X_{\text{cand }}$.

###  2.2 Reproducing kernel Hilbert spaces and choice of kernels

To implement the error analysis of an emulator and evaluate the scale of
expensive simulations, a reproducing kernel Hilbert space is a convenient
tool. The definition of RKHS is given in [12].

######  Definition 1 ([12], Definition 2.7.1.).

Let $\mathbb{H}$ be a Hilbert space composed of the real valued functions
defined on $\mathbb{D}$. A bivariate function $K$ on
$\mathbb{D}\times\mathbb{D}$ is said to be a reproducing kernel (RK) for
$\mathbb{H}$ if

  1. 1.

for every $t\in\mathbb{D},K(\cdot,t)\in\mathbb{H}$ and

  2. 2.

$K$ satisfies the reproducing property that $f(t)=\langle
f,K(\cdot,t)\rangle_{\mathbb{H}}$ for every $f\in\mathbb{H}$ and
$t\in\mathbb{D}$ where $\langle\cdot,\cdot\rangle_{\mathbb{H}}$ is the inner
product of $\mathbb{H}$.

When $\mathbb{H}$ possesses an RK it is said to be an RKHS.

In implementing GPR, a specific covariance structure, which exactly
corresponds to the reproducing kernel in an RKHS, should be chosen. The
structure of RKHS, or a variety of functions with which a Gaussian process
emulator can deal, is totally determined by the covariance structure. In other
words, we determine the function space where the true response lies by
choosing the covariance structure. If a reproducing kernel is poorly chosen,
we cannot expect posterior mean function of GP (the prediction) to be
approximated accurately. In this paper, Matérn covariance functions are used
due to clear interpretation of the corresponding function spaces. We start by
reviewing the following form of Matérn covariance functions.

$\displaystyle
K_{\nu,\lambda,\sigma^{2}}(x,x^{\prime})=\sigma^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{\|x-x^{\prime}\|_{2}}{\lambda}\right)^{\nu}J_{\nu}\left(\sqrt{2\nu}\frac{\|x-x^{\prime}\|_{2}}{\lambda}\right)$
(7)

where $\|x-x^{\prime}\|_{2}$ denotes the Euclidean distance between $x$ and
$x^{\prime}$ for $x,x^{\prime}\in\mathbb{R}^{d}$, $\Gamma_{\nu}$ is the Gamma
function for $\nu>0$ and $J_{\nu}$ is a modified Bessel function of order
$\nu$. The parameter $\sigma^{2}$ is usually referred to as the (marginal)
variance, $\lambda$ as the correlation length and $\nu$ as the smoothness
parameter. The expression for the Matérn covariance kernel simplifies for
particular choices of $\nu$. Notable examples include the exponential
covariance kernel
$\sigma^{2}\exp\left(-\left\|x-x^{\prime}\right\|_{2}/\lambda\right)$ with
$\nu=1/2,$ and the Gaussian covariance kernel
$\sigma^{2}\exp\left(-\left\|x-x^{\prime}\right\|_{2}^{2}/\lambda^{2}\right)$
in the $\operatorname{limit}\nu\rightarrow\infty$.

Choosing the hyperparemeters $\\{\nu,\lambda,\sigma^{2}\\}$ in Mat́ern kernel
functions is quite important. We have assumed $\delta_{l}$ belongs to the
corresponding RKHS $\mathbb{H}_{l}$ with RK $K_{l}$, which is crucial in the
approximation, but poorly chosen hyperparameters can ruin this assumption and
hence the approximation. Therefore, we have to grasp how they affects our
approximation. Let us denote $\mathbb{H}(\mathbb{R}^{d})$ be a RKHS defined on
$\mathbb{R}^{d}$ with Matérn kernel ${K_{\nu,\lambda,\sigma^{2}}}$ specified
in $(\ref{mat})$. Moreover, we define the Sobolev spaces
$H^{s}(\mathbb{R}^{d})$ as

$\displaystyle H^{s}(\mathbb{R}^{d})\coloneqq\left\\{f\in
L^{2}\left(\mathbb{R}^{d}\right):(1+\|\xi\|^{2})^{s/2}\widehat{f}(\xi)\in
L^{2}(\mathbb{R}^{d})\right\\}$ (8)

where $\hat{f}(\xi)$ is the Fourier transformation of $f(x)$. Then, we can see
the degree of the smoothness of functions in $\mathbb{H}(\mathbb{R}^{d})$ by
the following proposition.

######  Proposition 1 ([31], Proposition 3.2).

Let $k_{\nu,\lambda,\sigma^{2}}$ be a Matérn covariance kernel as defined in
(7). Then the RKHS $\mathbb{H}(\mathbb{R}^{d})$ with the reproducing kernel
$K_{\nu,\lambda,\sigma^{2}}$ is equal to the Sobolev space
$H^{\nu+d/2}(\mathbb{R}^{d})$ as a vector space, and the norm of
$\mathbb{H}(\mathbb{R}^{d})$ and $H^{\nu+d/2}(\mathbb{R}^{d})$ are equivalent.

It is seen that as a vector space, the RKHS $\mathbb{H}(\mathbb{R}^{d})$ of
the Matérn kernel $K_{\nu,\lambda,\sigma^{2}}$ is fully determined by the
smoothness parameter $\nu$. The parameters $\lambda$ and $\sigma^{2}$ do,
however, influence the constants in the norm equivalence of the RKHS norm and
the standard Sobolev norm. We will see the relationships of these parameters
and error analysis of our GP emulator. As discussed in [2], we treat an RKHS
$\mathbb{H}(\mathbb{D})$ as a restriction of $\mathbb{H}(\mathbb{R}^{d})$ and
this RKHS is equipped with the kernel being the restriction of $K$ to
$\mathbb{D}\times\mathbb{D}$. Therefore, the above argument is true in the
case of $\mathbb{H}(\mathbb{D})$ as well.

###  2.3 Multilevel Monte Carlo (MLMC)

Multilevel Monte Carlo (MLMC) [10], decomposes an expensive function into
increments, as a similar framework of our analysis. As discussed in [23], in
the context of multifidelity, less accurate but cheap simulators are used to
reduce the runtime where possible and more sophisticated and expensive
simulators are employed to preserve the accuracy of the emulator. Then, a
natural question is to know the optimal numbers of runs for each simulator. In
MLMC, the goal is to estimate the expectation $\mathbb{E}[P]$ for some random
variable $P$ using Monte Carlo simulation. Our goal is to build an emulator of
the functional form of the computer model. It yields expectations and
variances but goes much further as the entire distribution of predicted values
is supplied by the emulator.

The random variables $P_{0},\ldots,P_{L}$ are prepared and the goal is to
estimate $\mathbb{E}[P_{L}]$. $P_{0},\ldots,P_{L-1}$ approximate $P_{L}$ with
increasing accuracy, but also increasing cost. Then, we have the simple
identity by decomposing $P_{L}$:
$\mathbb{E}\left[P_{L}\right]=\mathbb{E}\left[P_{0}\right]+\sum_{\ell=1}^{L}\mathbb{E}\left[P_{\ell}-P_{\ell-1}\right]$.
Therefore we can use the following unbiased estimator for
$\mathbb{E}\left[P_{L}\right]:$
$N_{0}^{-1}\sum_{n=1}^{N_{0}}P_{0}^{(0,n)}+\sum_{\ell=1}^{L}\left\\{N_{\ell}^{-1}\sum_{n=1}^{N_{\ell}}\left(P_{\ell}^{(\ell,n)}-P_{\ell-1}^{(\ell,n)}\right)\right\\}$
where the inclusion of the level $\ell$ in the superscript $(\ell,n)$
indicates that independent samples are used at each level of correction. For a
fixed cost, the variance is minimized by choosing $N_{\ell}$ to minimize

$\displaystyle\sum_{\ell=0}^{L}\left(N_{\ell}^{-1}V_{\ell}+\mu^{2}N_{\ell}C_{\ell}\right)$
(9)

for some value of the Lagrange multiplier $\mu^{2}$. This gives
$N_{\ell}=\mu\sqrt{V_{\ell}/C_{\ell}}$. The formal statement of the above
discussion is stated in Theorem 2.1 in [10].

###  2.4 Uniform error bounds in GPR prediction

In the latter section, we implement an independent analysis of optimal numbers
of runs mixed with the concept of MLMC. A clear relationship between the
accuracy of an emulator and number of training data points was provided in
[34]. Matérn kernels (7) are assumed for the covariance structure of GPR and
we study the case of misspecified kernel functions. This addresses a common
question in GPR when the true kernel function is unknown. We denote
$\\{\nu_{0},\lambda_{0},\sigma_{0}^{2}\\}$ as the true parameters in (7) and
we denote the Matérn kernel in (7) as $K_{0}$ if the parameters are equal to
$\\{\nu_{0},\lambda_{0},\sigma_{0}^{2}\\}$. For $x\in\mathbb{D}$, we also
define a Power function $P_{K,X_{N}}$ as

$\displaystyle
P_{K,X_{N}}(x)=1-K\left(\boldsymbol{x},X_{N}\right)^{\top}K(X_{N},X_{N})^{-1}K\left(\boldsymbol{x},X_{N}\right)$
(10)

The Power function $P_{K,X_{N}}$ is the special case of predictive variance
$\hat{s}^{2}_{N}\left(\boldsymbol{x}\right)$ where $\sigma^{2}=1$. We further
define the supremum of the Power function as

$\displaystyle P_{K,X_{N}}=\sup_{x\in\mathbb{D}}P_{K,X_{N}}(x)$ (11)

The following condition for the kernel is necessary to present the convergence
result in [34].

######  Condition 1 ([34]).

The kernels $K$ and $K_{0}$ as radial basis functions are continuous and
integrable on $\mathbb{R}^{d}$, satisfying

$\|\hat{K}_{0}/\hat{K}\|_{L_{\infty}\left(\mathbb{R}^{d}\right)}=:A_{1}^{2}<+\infty$

where $\hat{K}_{0}$ and $\hat{K}$ are Fourier transformation of $K_{0}$ and
$K$. In addition, there exists $\tilde{\alpha}\in(0,1]$, such that

$\int_{\mathrm{R}^{d}}\|\omega\|^{\tilde{\alpha}}\hat{K}(\omega)d\omega=:A_{0}<+\infty$

Then, we can present the convergence result of [34]. It is important to see
that $\sigma$ in the following upper bound represents the scale of target
function $y$.

######  Theorem 1 ([34]).

Suppose Condition 1 holds, and the design set $X_{N}$ is dense enough in the
sense that $P_{K,X_{N}}$ defined in (11) is no more than some given constant
$C$.

For any $u>0$ with probability at least
$1-2\exp\left\\{-u^{2}/\left(2A_{1}^{2}\sigma^{2}P_{K,X_{N}}^{2}\right)\right\\}$,
the prediction error has the upper bound

$\displaystyle\sup_{x\in\mathbb{D}}\left|m_{N}^{y}\left(x\right)-\delta(x)\right|\leq\mathcal{K}\sigma
P_{K,X_{N}}\log^{1/2}\left(\mathrm{e}/P_{K,X_{N}}\right)+u$ (12)

Here the constants $C,\mathcal{K}>0$ depend only on
$\mathbb{D},\tilde{\alpha},A_{0}$, and $A_{1}$.

$K_{0}$ is a Matérn kernel with true smoothness parameter $\nu_{0}$ as
defined. As discussed in [34], it can be verified that Condition 1 holds if
and only if $0<\nu\leq\nu_{0}$. If we want this convergence result to hold
with probability at least $p$, $u$ should be equal to $\sigma
P_{K,X_{N}}\left(2A_{1}\log{2\over 1-p}\right)^{1/2}$.

##  3 Multilevel adaptive sequential design of computer experiments (MLASCE)

In this section, we present our novel multilevel adaptive sequential design of
computer experiments (MLASCE) in the framework of multi-fidelity GPR. The
novelty is an error analysis and associated optimization of the usage in
computational resources. The assumptions and formulation of this model are
dealt in Section 3.1 to 3.2. Section 3.3 presents an error analysis and
Section 3.4 illustrates the algorithm for implementing our method.

###  3.1 Incremental decomposition and reproducing kernel Hilbert spaces
(RKHS)

Although an elegant framework of multi-fidelity GPR or co-kriging was
established in [14], this type of method fails to evaluate the improvement
from simulations of higher fidelity in the sense of optimizing the accuracy of
emulator and usage of computational resources. Moreover, their emulators would
be inherently biased due to the lack of error analysis. Here, we construct the
necessary settings in order to capture the degree of benefit of more
sophisticated simulations and construct an asymptotically unbiased statistical
emulator. One way is to decompose computer simulations into basic
approximation and incremental functions. This is an analogy of the correlation
structure (1) of [14] with $\rho_{l}=1$ but essential for us to implement
analysis of the bias. Decomposing expensive function into increments is not
only necessary to evaluate the additional effects of sophisticated computer
simulations but also a common technique in computational statistics.

We assume $y$ is unknown and aim to approximate it by other deterministic
functions of different fidelity $y_{l}:\mathbb{D}\rightarrow\mathbb{R}$ for
level $l=1,\ldots,L$, which are called computer simulations. Each $y_{l}$ has
its own degree of accuracy denoted as $h_{l}$ and cost $t_{l}$. $h_{l}$
satisfies $1\geq h_{1}>h_{2}>\cdots>h_{L}$ and we assume $y_{l}$ becomes
identical with $y$ as $h_{l}\rightarrow 0$. $h_{l}$ is an ordinal number and
should be decided by the user in proportion to different accuracy of computer
simulations. For example, $h_{l}$ can be determined by the proportion of
different discretization levels in a numerical solver. The actual value of
$h_{l}$ is not used in MLASCE but is present only in the next section. The
cost $t_{l}$ satisfies $t_{L}>t_{L-1}>\cdots>t_{1}$ hence $y_{l_{1}}$ is more
accurate but expensive than $y_{l_{2}}$ for different levels $l_{1}>l_{2}$. We
introduce discrepancy term $\delta_{l}:\mathbb{D}\rightarrow\mathbb{R}$ for
$l=2,\ldots,L$ so that we establish the following relationships

$\displaystyle
y_{l}(\boldsymbol{x})=y_{l-1}(\boldsymbol{x})+\delta_{l}(\boldsymbol{x}).$

where $x$ is an arbitrary element in $\mathbb{D}$. Using this expression,
$y_{L}$ is decomposed into the sum of the basic approximation and increments.

$\displaystyle y_{L}(\boldsymbol{x})=\sum_{l=1}^{L}\delta_{l}(\boldsymbol{x})$
(13)

where $\delta_{1}(\boldsymbol{x})=y_{1}(\boldsymbol{x})$ and
$\delta_{l}(\boldsymbol{x})=y_{l}(\boldsymbol{x})-y_{l-1}(\boldsymbol{x})$ for
$l=2,\ldots,L$. We can see $\delta_{1}$ is the baseline approximation and
$\\{\delta_{l}\\}_{l=2}^{L}$ are the incremental refinements. It is expected
that the value of increment $\delta_{l}(\boldsymbol{x})$ should decrease as
the level $l$ increases thus less contributions come from higher fidelity
simulations. In order to evaluate accurately the degree of improvements, we
first assume each $\delta_{l}$ is an element of a reproducing kernel Hilbert
space (RKHS) $\mathbb{H}_{l}$ with a reproducing kernel $K_{l}$ for
$l=1,\ldots,L$. Let
$\boldsymbol{\delta_{l}}(X_{l,N_{l}})=\left(\delta_{l}(x_{l,1}),\delta_{l}(x_{l,2})\ldots,\delta_{l}(x_{l,N_{l}})\right)^{\top}$
denote $N_{l}$ outputs of $\delta_{l}$ collected at the inputs
$X_{l,N_{l}}=\left(x_{l,1},x_{l,2},\ldots,x_{l,N_{l}}\right)^{\top}$. Next, we
build a new multi-fidelity GPR by assuming that
$\boldsymbol{\delta_{l}}(X_{l,N_{l}})$ is a realization of GP
$\eta_{l}(\boldsymbol{x})\sim
GP\left(0,K_{l}(\boldsymbol{x},\boldsymbol{x}^{\prime})\right)$, where every
$\eta_{l}$ is independent of each other for different $l$. Based on the
experimental design or training data $X_{l,N_{l}}$ and
$\boldsymbol{\delta_{l}}(X_{l,N_{l}})$, we obtain the GPR prediction at $x$
consisting of the posterior distribution
$\eta_{l}\left(\boldsymbol{x}\right)|X_{l,N_{l}},\boldsymbol{\delta_{l}}(X_{l,N_{l}})\sim\mathcal{N}\left(m_{l,N_{l}}^{\delta_{l}}\left(\boldsymbol{x}\right),\hat{s}^{2}_{l,N_{l}}\left(\boldsymbol{x}\right)\right)$
with posterior mean and variance given by

$\begin{aligned}
m_{l,N_{l}}^{\delta_{l}}\left(\boldsymbol{x}\right)&=K_{l}\left(\boldsymbol{x},X_{l,N_{l}}\right)^{\top}K_{l}(X_{l,N_{l}},X_{l,N_{l}})^{-1}\boldsymbol{\delta_{l}}(X_{l,N_{l}})\\\
\hat{s}^{2}_{l,N_{l}}\left(\boldsymbol{x}\right)&=K_{l}(\boldsymbol{x},\boldsymbol{x})-K_{l}\left(\boldsymbol{x},X_{l,N_{l}}\right)^{\top}K_{l}(X_{l,N_{l}},X_{l,N_{l}})^{-1}K_{l}\left(\boldsymbol{x},X_{l,N_{l}}\right)\end{aligned}.$

Therefore, we can construct the prediction for $y_{L}(\boldsymbol{x})$ by
taking note that every $\eta_{l}$ is independent of each other and summing up
$\eta_{l}(\boldsymbol{x})|X_{N_{l}},\boldsymbol{\delta_{l,N_{l}}}$. The
predictive distribution of our multi-fidelity GPR is given by

$\displaystyle\mathcal{N}\left(m^{y_{L}}_{L}(\boldsymbol{x}),s^{2}_{L}(\boldsymbol{x})\right)$
(14)

where
$m^{y_{L}}_{L}(\boldsymbol{x})=\sum_{l=1}^{L}m_{l,N_{l}}^{\delta_{l}}\left(\boldsymbol{x}\right)$
and
$s^{2}_{L}(\boldsymbol{x})=\sum_{l=1}^{L}\hat{s}^{2}_{l,N_{l}}\left(\boldsymbol{x}\right)$.

The posterior mean function
$m_{l,N_{l}}^{\delta_{l}}\left(\boldsymbol{x}\right)$ is used for the
predicting $\delta_{l}(\boldsymbol{x})$ and it is an element of the
corresponding RKHS. We further expect that it would converge to the true
response in some sense if the number of data points tends to infinity. Theorem
2 guarantees this desired property and this is going to be discussed later.

If we use a non-zero mean $m_{l}(\cdot)$ in GP $\eta_{l}$, the formula for the
predictive mean $m_{l,N_{l}}^{\delta_{l}}$ changes to

$m_{l,N_{l}}^{\delta_{l}}\left(\boldsymbol{x}\right)=m_{l}(\boldsymbol{x})+K_{l}\left(\boldsymbol{x},X_{l,N_{l}}\right)^{\top}K_{l}(X_{l,N_{l}},X_{l,N_{l}})^{-1}\left(\boldsymbol{\delta_{l}}(X_{l,N_{l}})-m_{l}\left(X_{l,N_{l}}\right)\right)$

where
$m_{l}\left(X_{N_{l}}\right):=\left[m_{l}\left(x_{l,1}\right),\ldots,m\left(x_{l,N_{l}}\right)\right]\in\mathbb{R}^{N_{l}}$
and he predictive covariance is unchanged. The benefit of specifying a non-
zero mean function is that one could reduce “the amount of learning” by
considering the discrepancy $m_{l}-\delta_{l}$ instead of $\delta_{l}$ itself
if a specific form of $m_{l}$ is promising. On the other hand, parametrizing
$m_{l}$ would require large amount of parameters especially in the case of
high dimensional inputs.

###  3.2 Choice of smoothness parameter

As we can see in the Proposition 1, the smoothness parameter decides the
underlining function space to which the target function is assumed to belong.
However, choosing the true smoothness parameter $\nu_{0}$, which specifies the
smoothness of the target function, may be a difficult task. The
hyperparameters in a covariance function is usually unknown. We might be able
to estimate $\nu$ by maximum likelihood method or other statistical methods
but sufficient amount of data may not be available to obtain the reliable
result, which is a typical case where the computer simulation is expensive to
run. Moreover, evaluating the Matérn kernel with $\nu$ other than
$0.5,1.5,2.5,\infty$ numerous times in the process of maximum likelihood
method or experimental design is considerably expensive and it is common to
fix the smoothness parameter $\nu$ before implementing the computation in the
packages of statistical software. Therefore, we usually have to set the value
for $\nu$ in advance. For the different smooth parameters $\nu$ and
$\nu^{\prime}$ $(\nu<\nu^{\prime})$, the corresponding Sobolev spaces satisfy
$H^{\nu^{\prime}}(\mathbb{R}^{d})\subseteq H^{\nu}(\mathbb{R}^{d})$ hence one
possible strategy of choosing $\nu$ is to set a small value.

Some papers investigate the effect of misspecified smoothness parameters on
the GPR approximation in the context of space filling design. The convergence
rate of the posterior mean of GPR becomes slower if the (fixed) smoothness
parameters are different from the true value [34, 33]. Similar results [32,
39] are available in the case that the estimated smoothness parameter is used.
Although we do not use space filling design for locating the training data,
these results are insightful for choosing the smoothness parameter.

###  3.3 Error analysis

####  3.3.1 Convergence of mean prediction

The estimated parameters of $\lambda,\sigma^{2}$ and a priori fixed $\nu$ are
used in the Matérn function. The estimation of $\lambda,\sigma^{2}$ are based
on the data obtained by design of experiment. In this subsection, we show that
our multifidelity GPR is a reasonable approximation of the target function. In
this framework, the training data is obtained by implementing MICE as the
design of experiment. It is assumed that the observation vector
$\boldsymbol{y}$ is a realization of the Gaussian process
$Y(\boldsymbol{x})\sim
GP\left(0,K_{\nu,\theta}\left(\boldsymbol{x},\boldsymbol{x}^{\prime}\right)\right)$
where $K_{\nu,\theta}=K_{\nu,\lambda,\sigma^{2}}$ is a Matérn kernel defined
in $(\ref{mat})$ and $\theta$ denotes $\\{\lambda,\sigma^{2}\\}$ for
convenience. The GPR prediction at $x$ are obtained in the form of the normal
distribution
$y\left(x\right)\sim\mathcal{N}\left(m_{N,\theta}^{y}\left(x\right),\hat{s}^{2}_{N,\theta}\left(x\right)\right),$
with posterior mean and variance given by

$\begin{aligned}
m_{N,\theta}^{y}\left(\boldsymbol{x}\right)&=K_{\nu,\theta}\left(\boldsymbol{x},X_{N}\right)^{\top}K_{\nu,\theta}(X_{N},X_{N})^{-1}\boldsymbol{y}(X_{N})\\\
\hat{\mathrm{s}}^{2}_{N,\theta}\left(\boldsymbol{x}\right)&=K_{\nu,\theta}(\boldsymbol{x},\boldsymbol{x})-K_{\nu,\theta}\left(\boldsymbol{x},X_{N}\right)^{\top}K_{\nu,\theta}(X_{N},X_{N})^{-1}K_{\nu,\theta}\left(\boldsymbol{x},X_{N}\right)\end{aligned}.$

We are interested in a pointwise prediction of $y\left(\boldsymbol{x}\right)$
and it is a natural question to investigate the error between the
$y\left(\boldsymbol{x}\right)$ and its mean prediction
$m_{N,\theta}^{y}\left(\boldsymbol{x}\right)$. To formulate this analysis, let
$\Theta$ be a compact set in $\mathbb{R}^{2}$ and $y\in\mathbb{H}$ for an
arbitrary $\theta\in\Theta$ where $\mathbb{H}$ is a RKHS with the RK
$K_{\nu,\theta}$. For any $\boldsymbol{x}\in\mathbb{D}$ and $\theta\in\Theta$,
we can easily check the error is given by

$\displaystyle|y\left(\boldsymbol{x}\right)-m_{N,\theta}^{y}\left(\boldsymbol{x}\right)|\leq\|y\|_{\mathbb{H}}\bigl{(}\hat{\mathrm{s}}^{2}_{N}\left(\boldsymbol{x}\right)\bigr{)}^{\frac{1}{2}}$
(15)

Moreover, the (pointwise) error converges to $0$ if enough size of data is
available. Although it may be argued that the amount of data for an expensive
computer simulation is limited in practice, the following theorem still
ensures that the approximation goes in the right direction. The proof is given
in the Appendix.

######  Theorem 2.

Let $\mathbb{H}_{\theta}$ denote an RKHS defined in
$\mathbb{D}\subseteq\mathbb{R}^{d}$ with the reproducing kernel
$K_{\nu,\theta}=K_{\nu,\lambda,\sigma^{2}}$ defined in (7) and
$\mathbb{R}^{2}_{+}$ denote the subset of $\mathbb{R}^{2}$ composed of
positive values. Further, we assume the following conditions.

  1. 1.

$\Theta$ is a compact set in $\mathbb{R}^{2}_{+}$.

  2. 2.

$y\in\mathbb{H}_{\theta}$ for an arbitrary $\theta\in\Theta$ and a fixed
$\nu$.

  3. 3.

$X=\\{x_{i}:1\leq i\leq\infty\\}$ is a set of points in $\mathbb{D}$ such that
$\text{span}\\{K_{\nu,\theta}(\cdot,x_{i}):1\leq i\leq\infty\\}$ is dense in
$\mathbb{H}_{\theta}$ and $y(X)=\\{y(x_{i}):1\leq i\leq\infty\\}$ is the
values of $y$ located at $X$.

  4. 4.

$X_{N}=\left(x_{1},x_{2},\ldots,x_{N}\right)^{\top}$ and
$\boldsymbol{y}(X_{N})=\bigl{(}y(x_{1}),y(x_{2})\ldots,y(x_{N})\bigr{)}$ are
subset of $X$ and $y(X)$.

  5. 5.

$K_{\nu,\theta}(X_{N},X_{N})^{-1}$ exists for any $\theta\in\Theta$ and
$X_{N}$.

Then, the posterior covariance
$K_{\nu,\theta}(x,u)-K_{\nu,\theta}(x,X_{N})K_{\nu,\theta}(X_{N},X_{N})^{-1}K_{\nu,\theta}(X_{N},u)\to
0$ and posterior mean
$K_{\nu,\theta}(x,X_{N})K_{\nu,\theta}(X_{N},X_{N})^{-1}y(X_{N})\to y(x)$ for
any $x,u\in\mathbb{D}$, $\theta\in\Theta$ and a fixed $\nu$ as $N\to\infty$.

In practice, the hyperparameters $\theta$ or $\\{\lambda,\sigma^{2}\\}$ in
$K_{\nu,\theta}$, which are usually unknown, are identified with maximum
likelihood estimator $\hat{\theta}_{N}$ whereas $\nu$ is set to a specific
fixed value in advance. One might argue that the changing value of
$\widehat{\theta}_{N}$ in
$m_{N,\hat{\theta}_{N}}^{y}\left(\boldsymbol{x}\right)$ as $N$ grows might not
be trivial but it can be verified that
$m_{N,\hat{\theta}_{N}}^{y}\left(\boldsymbol{x}\right)$ also converges to
$y(x)$ if we further assume a mild condition on
$m_{N,\theta}^{y}\left(\boldsymbol{x}\right)$ and that
$\\{\hat{\theta}_{N}\\}_{N=1}^{\infty}$ is in $\Theta$ with a fixed limit
$\theta_{0}\in\Theta$ (see Remark 2 in Appendix).

####  3.3.2 Applying convergence of mean to the multilevel setting

We are going to apply the above results to our multi-fidelity GPR. First, we
prepare Matérn kernels $K_{\nu_{l},\lambda_{l},\sigma^{2}_{l}}$ with fixed
$\nu_{l}$ for $l=1,\ldots,L$. We estimate $\\{\lambda_{l},\sigma_{l}^{2}\\}$
based on the design $\\{X_{l,N_{l}},\boldsymbol{y}(X_{l,N_{l}})\\}$
accordingly and set these values to the hyperparameters in
$K_{\nu_{l},\lambda_{l},\sigma^{2}_{l}}$. We plug these Matérn kernels in the
covariance structure of $(\ref{mymodel})$. We also assume that
$\delta_{l}\in\mathbb{H}_{l}$ with the RK
$K_{\nu_{l},\lambda_{l},\sigma^{2}_{l}}$ for a chosen $\nu_{l}$ and any
$\lambda_{l},\sigma^{2}_{l}$ in a compact set
$\Theta_{l}\subseteq\mathrm{R}^{2}$. Further, the estimated hyperparameters of
$\lambda_{l},\sigma^{2}_{l}$ are in $\Theta_{l}$. Then, we derive an error
bound for our mean prediction $m^{y_{L}}_{L}(x)$ in $(\ref{mymodel})$.

$\displaystyle|m^{y_{L}}_{L}(\boldsymbol{x})-y_{L}(\boldsymbol{x})|\leq\sum_{l=1}^{L}\bigl{(}\hat{s}^{2}_{l,N_{l}}\left(\boldsymbol{x}\right)\bigr{)}^{\frac{1}{2}}\|\delta_{l}\|_{\mathbb{H}_{l}}$
(16)

We write $m_{N,\theta}^{y}$ as $m_{N}^{y}$ and $K_{\theta}(\cdot,\cdot)$ as
$K(\cdot,\cdot)$ when there is no need for expressing the dependency on
specific $\theta$.

Further, Theorem 2 applies to each
$|m_{l,N_{l}}^{\delta_{l}}\left(\boldsymbol{x}\right)-\delta_{l}(\boldsymbol{x})|$
hence we are able to make the error as small as possible if sufficient amount
of data can be obtained. In practice, we usually cannot afford to run
expensive simulations to obtain enough data. Therefore, the next task is
optimizing the number of training data points
$\\{X_{l,N_{l}},\delta_{l}(X_{l,N_{l}})\\}$ across every $l$ under a
constraint of computational resources to achieve small upper bound in
$(\ref{myerror})$.

###  3.4 Algorithm for MLASCE and optimal usage of computational resources

Let us go back to the error bound in $(\ref{myerror})$. Our goal is to
implement the method of experimental design so that the above upper bound is
made as small as possible under a restriction of the computational resources.
That is, we prepare the training data
$\\{X_{l,N_{l}},\boldsymbol{\delta}_{l,N_{l}}(X_{l,N_{l}})\\}$ for all levels
$l=1,\ldots,L$ in the way that the error bound becomes as small as possible.
The error bound depends on not only the number of data points $N_{l}$ but also
the locations of the data $X_{N_{l}}$. However, evaluating it in an explicit
form without specific knowledge is a rather formidable task. Thus, an
optimization strategy should be developed. The essence of our strategy is to
choose the most effective level $l$ for reducing the overall error in
$(\ref{myerror})$ and run $\delta_{l}$ to obtain the new training data
$x_{N_{l}+1}$ and $\delta_{l}(x_{N_{l}+1})$ using MICE.

Although
$\hat{s}^{2}_{l,N_{l}}\left(x\right)^{\frac{1}{2}}\|\delta_{l}\|_{\mathbb{H}_{l}}$
tends to 0 for every $l$ as $N_{l}$ goes to infinity, its weight in the
overall upper bound (16) is different because of the scale
$\|\delta_{l}\|_{\mathbb{H}_{l}}$. Thus, it is sensible to get more samples at
the level where the reduction in
$\hat{s}^{2}_{l,N_{l}}\left(x\right)^{\frac{1}{2}}\|\delta_{l}\|_{\mathbb{H}_{l}}$
is large as $N_{l}$ increases. This is to ensure that the upper bound of (16)
becomes smaller quickly. $\|\delta_{l}\|_{\mathbb{H}_{l}}$ is generally
unknown but the convergence of $\|m_{l,N_{l}}^{\delta_{l}}\|_{\mathbb{H}_{l}}$
is equivalent to that of $\hat{s}^{2}_{l,N_{l}}\left(x\right)^{\frac{1}{2}}$
as discussed in Theorem 2 and it also includes the information of the weight
$\|\delta_{l}\|_{\mathbb{H}_{l}}$. Moreover,
$\|m_{l,N_{l}}^{\delta_{l}}\|_{\mathbb{H}_{l}}$ is bounded by
$\|\delta_{l}\|_{\mathbb{H}_{l}}$ as a non-decreasing function of $N_{l}$ (see
Proposition 3.1 in [32]). Therefore, this optimization problem is equivalent
of maximizing
$\sum_{l=1}^{L}\|m_{l,N_{l}}^{\delta_{l}}\|_{\mathbb{H}_{l}}^{2}$ under the
restriction of $\sum_{l=1}^{L}N_{l}t_{l}=T$. It is similar to the context of
nonlinear knapsack problem or resource allocation problems [22, 4] and our
optimization strategy can be based on the greedy methods, which is one of the
common approaches in this discipline (for example, see [21, 3, 9]). In our
case, the change in the value
$\|m_{l,N_{l}}^{\delta_{l}}\|_{\mathbb{H}_{l}}^{2}$ as $N_{l}$ is increased to
$N_{l}+1$ is seen as a criterion or “score”. When increasing the number of
points $N_{l}$, we decide which $l$ should be chosen based on the score of
every level and continue this process recursively. We can compute the actual
value of $\|m_{l,N_{l}}^{\delta_{l}}\|_{\mathbb{H}_{l}}$ by the relationship
$\|m_{l,N_{l}}^{\delta_{l}}\|_{\mathbb{H}_{l}}^{2}=\boldsymbol{\delta}_{l,N_{l}}(X_{l,N_{l}})^{\top}K_{l}(X_{l,N_{l}},X_{l,N_{l}})^{-1}\boldsymbol{\delta}_{l,N_{l}}(X_{l,N_{l}})$.

Let us define the “score” $\gamma_{l,N_{l}}$ by

$\displaystyle\gamma_{l,N_{l}}\coloneqq|\|m_{l,N_{l}}^{\delta_{l}}\|_{\mathbb{H}_{l}}^{2}-\|m_{l,N_{l}-1}^{\delta_{l}}\|_{\mathbb{H}_{l}}^{2}|\times\frac{a_{l}}{t_{l}}$
(17)

where $a_{l}$ is a tuning parameter which decides the degree of penalty of
cost $t_{l}$. It takes positive value up to $t_{l}$ and smaller value means
less emphasis on that level. The idea is to see the score as a benefit of
running $\delta_{l}$ weighted by its cost $t_{l}$. This kind of heuristic
strategy is widely proposed as a feasible solution to the nonlinear knapsack
problem such as in [21, 3, 9].

$a_{l}=t_{l}$ may be recommended as a first attempt since too much emphasis on
cheaper simulations would occur if the computational costs are hugely
different. As we shall see later, another existing strategy of sequential
design [18] yields inefficient emulators for the same reason. Even if the
lower levels of $\delta_{l}(x)$ are not emphasized by $a_{l}$,
$\|m_{l,N_{l}}^{\delta_{l}}\|_{\mathbb{H}_{l}}^{2}$ (hence $\gamma_{l,N_{l}}$)
tends to take smaller values for higher levels $l$ since
$\delta_{l}(x)=y_{l}(x)-y_{l-1}(x)$ becomes closer to $0$. In other words, our
algorithm with $a_{l}=t_{l}$ automatically tends to emphasize lower levels of
simulations since the fact $y_{l}(x)\simeq y_{l-1}(x)$ for a bigger $l$ gives
$\gamma_{l,N_{l}}$ a smaller value. In the case that $y_{l}(x)$ behaves very
differently than $y_{l-1}(x)$ does, which can happen in fluid dynamics
simulations with kinky movements, a wisely chosen weight $a_{l}$ may improve
the algorithm.

The available computational budget is denoted as $T$, often the available
computational time. We have to pay the computational cost $t_{l}$ for running
$\delta_{l}$ and $T$ is updated to $T-t_{l}$ accordingly if it is run. We
repeat this process until the computational budget depletes. The score
$\gamma_{l,N_{l}}$ is used for computing the potential reduction in every
level. We choose the level where $\gamma_{l,N_{l}}$ is the largest. It should
be noticed that $\delta_{l}$ cannot be run if $t_{l}$ becomes larger than $T$.
Thus, we have to grasp which computer simulations are available at each step.
The summary of our strategy, MLASCE, is as follows.

Preparation (1)

For $l=1,\ldots,L$, we prepare the computer simulation $\delta_{l}$ with the
computational cost $t_{l}$ for each run ($t_{1}<\cdots<t_{L}$), discretized
input domain $\mathrm{D}_{G,l}\subseteq\mathbb{R}^{d}$ and candidate set
$X_{\text{cand},l}\subseteq\mathrm{D}_{G,l}$. We denote the total
computational budget $T$ and $X_{l,k_{l}}\subseteq X_{\text{cand},l}$ the
experimental design at level $l$ of size $k_{l}$.

Preparation (2)

We prepare GP$(0,K_{l}(\cdot,\cdot))$ for $l=1,\ldots,L$. $K_{l}$ is a Matérn
kernel with some fixed $\nu_{l}$ and the other hyperparametrs $\sigma_{l}$ and
$\lambda_{l}$ updated during the procedure of MICE. We first set the score
$\gamma_{l,X_{l,0}}=0$ before obtaining any data.

Step 1

By implementing random sampling for every level, we obtain initial design
(training data) with the size $1$ for every $l$. That is, let $N_{l}=1$,
sample $\bigl{(}X_{l,1},\delta_{l}(X_{l,1})\bigr{)}$ and the hyperparameters
$\sigma^{2}_{l}$ and $\lambda_{l}$ of $K_{l}$ are updated based on the
obtained data.

Step 2

Calculate $\gamma_{l,X_{l,k_{l}}}$ based on the current design
$\bigl{(}X_{l,k_{l}},\delta_{l}(X_{l,k_{l}})\bigr{)}$ for every $l$. The
candidate sets are updated accordingly so that $X_{\text{cand},l}\leftarrow
X_{\text{cand},l}\setminus X_{l,k_{l}}$ for every $l$.

Step 3

We choose the level $l$ where the score $\gamma_{l,X_{l,k_{l}}}$ is the
largest under the condition $t_{l}\leq T$.

Step 3

At level $l$, sample a next training point $\\{x_{l},\delta_{l}(x_{l})\\}$
using MICE. Set $X_{l,k_{l}+1}\leftarrow X_{l,k_{l}}\cup x_{l}$. The
hyperparameters $\sigma^{2}_{l}$ and $\lambda_{l}$ are updated by fitting the
new data $\\{X_{l,k_{l}+1},\delta_{l}(X_{l,k_{l}+1})\\}$ to
GP$(0,K_{l}(\cdot,\cdot))$.

Step 4

If $T<t_{1}$, then stop; otherwise go to Step 2.

Output

: The experimental design $\\{X_{l,N^{*}_{l}},\delta_{l}(X_{l,N^{*}_{l}})\\}$
for every $l$ and the predictive distribution (14).

By implementing this algorithm, the Gaussian process emulator is automatically
built since the hyperparameters $\sigma^{2}_{l}$ and $\lambda_{l}$ are updated
whenever new data is sampled.

Figure 1: Workflow of MLASCE. New design point is generated at level 1 first
and next at level $L$.

##  4 The special case of a uniform tensor grid: optimal numbers of runs

Although MLASCE produces an efficient statistical emulator, the number of data
points for each computer simulation is actually unknown until the algorithm
finishes its computation. It would be beneficial to get a sense of optimal
size of experimental design before running computer simulations albeit in the
case of space filling design. This is the case if locations of design are
specified by prior knowledge of a system. In this section, an analysis of
optimal numbers of runs (prior to actual implementation of computer
simulations) is provided and this is independent of the algorithm of MLASCE
presented in previous sections. The scheme presented in this part is similar
to that of MLMC (9): the error of the emulator is expressed as a function of
the number of data points $N_{l}$ and it is minimized under a restriction. The
assumptions and notations are same as in the previous section. The main
principle of this analysis is composed of (i) deriving a convergence rate of
$|m_{l,N_{l}}^{\delta_{l}}\left(x\right)-\delta_{l}(x)|$ in (16) with respect
to the number of runs $N_{l}$ and (ii) minimizing the upper bound of (16) as a
function of $N_{l}$. There is rich literature regarding obtaining the
convergence rate in the field of scattered data approximation (e.g. [35, 36]).
An approximation is constructed using radial basis functions originating from
reproducing kernels and the approximation coincides with the mean of a
conditioned GP. The results from this literature are hence directly
transferable to approximation with GPs. The work of [34] and [32] contains the
most up-to-date application of scattered data approximation results for GPR.
Given the limited size of training data, we are more interested in a pointwise
prediction than an error bound in a Sobolev norm. In [34], an upper bound of
worst-case prediction error is given, which is to be applied to our analysis.

We start from the following error decomposition (18). We use Matérn kernels
(7) with hyperparameters $\\{\nu_{l},\lambda_{l},\sigma_{l}^{2}\\}$ again for
$m_{l,N_{l}}^{\delta_{l}}\left(x\right)$ for every $l$ and assume the case of
misspecified kernel functions. This analysis of optimal numbers of runs is
implemented before obtaining training data and there is no clue to know the
actual values of parameters in the kernels hence this assumption is crucial.
The outline of this study is (i) the error bound in (18) is converted into the
function of Power functions based on Theorem 1, (ii) Power functions are
further converted into the fill distance and finally (iii) the fill distance
is translated into the $N_{l}$.

$\displaystyle|y_{L}(\boldsymbol{x})-m^{y_{L}}_{L}(\boldsymbol{x})|\leq\sum_{l=1}^{L}|m_{l,N_{l}}^{\delta_{l}}\left(\boldsymbol{x}\right)-\delta_{l}(\boldsymbol{x})|$
(18)

Then, we present the theorem as the result of above direction. To simplify the
notation, we write $a\lesssim b$ for two positive quantities $a$ and $b,$ if
$a/b$ is uniformly bounded independent of the number of points $N_{l}$ and the
degree of accuracy $h_{l}$ of $\delta_{l}$. We also write $a\simeq b,$ if
$a\lesssim b$ and $b\lesssim a$. The definition of the fill distance is given
after this theorem.

######  Theorem 3.

We assume the following conditions.

  1. 1.

$\Omega\subseteq\mathbb{R}^{d}$ is compact and convex with a positive Lebesgue
measure

  2. 2.

$X_{l,N_{l}}$ is a uniform tensor grid for $l=1,\ldots,L$.

  3. 3.

Condition 1 is satisfied for Matérn kernels $K_{l}$ and $K_{l,0}$ in (7) for
$l=1,\ldots,L$.

  4. 4.

$\sigma_{l}\lesssim(h_{l}-h_{l-1})^{2\alpha}$ holds for some positive $\alpha$
and $l=1,\ldots,L$ where $h_{0}$ is $0$.

  5. 5.

$\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}$ in Lemma 4 takes value in
$(0,\sqrt{\mathrm{e}}]$ and the fill distance
$h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}$ is smaller than $\frac{1}{e}$ for
$l=1,\ldots,L$.

Then, the following inequality holds for probability at least $p$ where $p$ is
an arbitrary number in $[0,1)$.

$\displaystyle|y_{L}(x)-m^{y_{L}}_{L}(x)|$ $\displaystyle\lesssim$
$\displaystyle\sum_{l=1}^{L}\Bigl{(}(h_{l}-h_{l-1})^{2\alpha}N_{l}^{-{\nu_{l}\over
d}}\log^{1/2}N_{l}\Bigr{)}$ (19)

Before showing the detail of proof, the auxiliary lemma is presented so that
the Power function is into the function of the fill distance, which is a
quantity depending only on the design $X_{l,N_{l}}$. Given the experimental
region $\mathbb{D}$, the fill distance of a design $\mathbf{X}$ is defined as
$h_{\mathrm{X}}\coloneqq\sup_{x\in\Omega}\min_{x_{j}\in\mathrm{X}}\left\|x-x_{j}\right\|.$
Clearly, the fill distance quantifies the space-filling property [28] of a
design. A design having the minimum fill distance among all possible designs
with the same number of points is known as a minimax distance design [13].
Thanks to the following Lemma, the unknown parameter $\lambda_{l}$ in the
kernel is incorporated into the constant $\tilde{c}_{l},h_{l,0}$.

######  Lemma 4 ([38], Theorem 5.14).

Let $\mathbb{D}$ be compact and convex with a positive Lebesgue measure;
$K_{l}(x,x^{\prime})$ be a Matérn kernel given by (7) with the smoothness
parameter $\nu_{l}$. Then there exist constants $\tilde{c}_{l},h_{l,0}$
depending only on $\mathbb{D},\nu_{l}$ and the scale parameter $\lambda_{l}$
in (7), such that
$P_{K_{l},X_{l,N_{l}}}\leq\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}$
provided that $h_{\mathrm{X}_{l,N_{l}}}\leq h_{l,0}$

###### Proof of Theorem 3.

By applying Theorem 1 and Lemma 4 to
$|m_{l,N_{l}}^{\delta_{l}}\left(x\right)-\delta_{l}(x)|$, we obtain the
following inequality with probability at least $p$ :

$\displaystyle\sup_{x\in\mathbb{D}}\left|m_{l,N_{l}}^{\delta_{l}}\left(x\right)-\delta_{l}(x)\right|$
$\displaystyle\leq$
$\displaystyle\mathcal{K}_{l}\sigma_{l}P_{K_{l},X_{l,N_{l}}}\log^{1/2}\left(\mathrm{e}/P_{K_{l},X_{l,N_{l}}}\right)+u$
(20) $\displaystyle\leq$
$\displaystyle\mathcal{K}_{l}\sigma_{l}\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\log^{1/2}\left(\mathrm{e}/\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\right)+u$
$\displaystyle\leq$
$\displaystyle\mathcal{K}_{l}\sigma_{l}\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\log^{1/2}\left(1/h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\right)^{\kappa_{l}}+u$

In the second inequality, we used the fact the function
$x\log^{1/2}(\mathrm{e}/x)=x(1-\log x)^{1/2}$ is monotone increasing in
$(0,\sqrt{\mathrm{e}}]$ and assumption that
$\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}$ falls in this interval.
Besides,
$(\mathrm{e}/\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}})\leq(1/h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}})^{\kappa_{l}}$
holds for $h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}<\frac{1}{e}$ and some positive
$\kappa_{l}$.

For uniform tensor grids $X_{l,N_{l}}$, the fill distance
$h_{\mathrm{X}_{l,N_{l}}}$ is of the order $N_{l}^{-1/d}$ (Proposition 14.1 in
[35]). Further, we choose u to be
$\sigma_{l}P_{K_{l},X_{l,N_{l}}}\left(2A_{1,l}\log{2\over 1-p}\right)^{1/2}$
so that the inequality holds for probability at least $p$. Then, the upper
bound in (20) is written in the following:

$\displaystyle\mathcal{K}_{l}\sigma_{l}\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\log^{1/2}\left(1/h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\right)^{\kappa_{l}}+u$
$\displaystyle=$
$\displaystyle\mathcal{K}_{l}\sigma_{l}\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\log^{1/2}\left(1/h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\right)^{\kappa_{l}}+\sigma_{l}P_{K_{l},X_{l,N_{l}}}\left(2A_{1,l}\log{2\over
1-p}\right)^{1/2}$ $\displaystyle\leq$
$\displaystyle\mathcal{K}_{l}\sigma_{l}\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\log^{1/2}\left(1/h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\right)^{\kappa_{l}}+\sigma_{l}\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\left(2A_{1,l}\log{2\over
1-p}\right)^{1/2}$ $\displaystyle\leq$
$\displaystyle\mathcal{K}_{l}\sigma_{l}\tilde{c}_{l}\kappa_{l}^{1/2}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\log^{1/2}\left(1/h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\right)+\sigma_{l}\left(2A_{1,l}\log{2\over
1-p}\right)^{1/2}\tilde{c}_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\log^{1/2}\left(1/h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\right)$
$\displaystyle\simeq$
$\displaystyle\sigma_{l}h_{\mathrm{X}_{l,N_{l}}}^{\nu_{l}}\log^{1/2}\left(1/h_{\mathrm{X}_{l,N_{l}}}\right)$
$\displaystyle\simeq$ $\displaystyle\sigma_{l}N_{l}^{-{\nu_{l}\over
d}}\log^{1/2}N_{l}$ $\displaystyle\lesssim$
$\displaystyle(h_{l}-h_{l-1})^{2\alpha}N_{l}^{-{\nu_{l}\over
d}}\log^{1/2}N_{l}$

In the fourth line, we used the fact that $x\leq x\log^{1/2}\frac{1}{x}$ for
$x<\frac{1}{e}$. We finally arrive at the conclusion by summing up the above
upper bound :

$\displaystyle|y_{L}(x)-m^{y_{L}}_{L}(x)|$ $\displaystyle\leq$
$\displaystyle\sum_{l=1}^{L}|m_{l,N_{l}}^{\delta_{l}}\left(x\right)-\delta_{l}(x)|$
$\displaystyle\lesssim$
$\displaystyle\sum_{l=1}^{L}\Bigl{(}(h_{l}-h_{l-1})^{2\alpha}N_{l}^{-{\nu_{l}\over
d}}\log^{1/2}N_{l}\Bigr{)}.$

∎

The comments on the assumptions of $p$ and $\sigma_{l}$ are made here. The
probability $p$ is no longer present in the upper bound since the constants
are independent of order of convergence but $p$ does affect the actual
convergence speed. Besides, we cannot know the exact or estimated value before
running computer simulations though $\sigma_{l}$ represents the scale (hence
importance) of the increment $\delta_{l}$. Therefore, we assume $\sigma_{l}$,
which implies the scale of increment $\delta_{l}$, is proportional to the
degree of accuracy $h_{l}$ and $h_{l-1}$ corresponding to the increment
$\delta_{l}$.

Our final goal is to minimize the (equivalent) upper bound of (19) under the
restriction of computational budget $\sum_{l=1}^{L}N_{l}t_{l}=T$ and
$N_{l}\geq 1$ for every $l$. We are now ready to present the following
proposition.

######  Proposition 2.

Under the same conditions in Theorem 3, the optimal numbers of runs $N_{l}$ is
given by the solution of the following optimization problem:

$\displaystyle\mathop{\rm arg~{}min}\limits_{N_{1},\cdots,N_{L}}$
$\displaystyle\
\sum_{l=1}^{L}\Bigl{(}(h_{l}-h_{l-1})^{2\alpha}N_{l}^{-{\nu_{l}\over
d}}\log^{1/2}N_{l}\Bigr{)}$ $\displaystyle s.t.$
$\displaystyle\sum_{l=1}^{L}N_{l}t_{l}=T$ $\displaystyle N_{l}\geq 1\hskip
10.0ptl=1,\ldots,L$

where the optimal numbers of runs mean that the equivalent quantity of the
upper bound $|y_{L}(x)-m^{y_{L}}_{L}(x)|$ is minimized with respect to
$N_{l}$.

The contributions from $\delta_{l}$ for large $l$ is weighted with relatively
small value $\sigma_{l}^{2}$, which is proportional to the accuracy level
$(h_{l}-h_{l-1})^{2}$. Hence we do not have to run expensive computer
simulations many times, which is the same trick in MLMC as well. The
optimization problem would not be solved in a closed form in general but we
can apply an ordinary method of numerical optimization such as Trust-Region
Constrained Algorithm or BFGS Algorithm. The choice of $\alpha$ might be a
challenge while $h_{l}$ can be determined by the user, typically taken to be
proportional to discretization level in a numerical solver. In the similar
frame work of MLMC [10], deciding unknown parameters is also difficult task. A
feasible strategy for choosing $\alpha$ in our analysis is to select a
relatively small value. For $\alpha<\alpha^{\prime}$, if
$\sigma_{l}\lesssim(h_{l}-h_{l-1})^{2\alpha^{\prime}}$ holds then
$\sigma_{l}\lesssim(h_{l}-h_{l-1})^{2\alpha}$ is satisfied since
$h_{i}-h_{i-1}$ is smaller than 1 as defined before. The Figure 2(b) show how
the numbers of data points at different levels change as $\alpha$ takes
different value. The more emphasis is placed on the lower levels as $\alpha$
takes lager values, which is consistent with our expectation. We assume that
the number of levels $L$ is $5$, $d=1$, $\nu_{l}=2.5$ for every $l$ and total
budget $T$ is $800$. The degree of accuracy$h_{l}$ and computational costs
$t_{l}$ are in Table 1. The rate of increase in the cost is exponential and
this assumption is similar to that of the theorem of MLMC (see Theorem 2.1 in
[10]).

Table 1: Computational costs of computer simulations | $l=1$ | $l=2$ | $l=3$ | $l=4$ | $l=5$
---|---|---|---|---|---
$h_{l}$ | $1$ | $1/2$ | $1/4$ | $1/8$ | $1/16$
$t_{l}$ | $0.5$ | $2$ | $8$ | $32$ | $128$

Larger value of $\alpha$ means the more importance is placed on $\delta_{l}$
of lower accuracy hence the smaller number of samples is picked from expensive
simulations. Unless the simulations of low fidelity is known to be trustworthy
in advance, the small $\alpha$ is recommended.

######  Remark 1.

If we assume $\nu_{l}=\nu^{*}$ for every $l$ and $\nu^{*}>d/2e$, a closed form
of the numbers of sample points can be derived. By taking note that $\log
x\leq x^{a}$ for $x\geq 1$ and $a\geq 1/e$, the upper bound in (19) is
replaced with the more loose bound
$\sum_{l=1}^{L}\Bigl{(}(h_{l}-h_{l-1})^{2\alpha}N_{l}^{-{\nu^{*}\over
d}+{1\over 2e}}\Bigr{)}$. Therefore, the optimization problem in Proposition 2
can be replaced with $\mathop{\rm
arg~{}min}\limits_{N_{1},\cdots,N_{L}}\sum_{l=1}^{L}\Bigl{(}(h_{l}-h_{l-1})^{2\alpha}N_{l}^{-{\nu^{*}\over
d}+{1\over 2e}}\Bigr{)}$ under the restriction $\sum_{l=1}^{L}N_{l}t_{l}=T$
for positive $N_{l}$. By applying Lagrange multipliers, the optimal numbers of
points are given by

$N_{l}^{*}=\Bigl{(}\frac{t_{l}}{-r(h_{l}-h_{l-1})^{2\alpha}}\Bigr{)}^{\frac{1}{r-1}}\frac{T}{\sum_{j=1}^{L}t_{j}^{\frac{r}{r-1}}(-r(h_{l}-h_{l-1})^{2\alpha})^{\frac{1}{1-r}}}$

where $r=-\nu^{*}/d+1/2e$ $(<0)$.

|
---|---

(a)

(b)

Figure 2: The number of runs of at each level ( $L=5$). The blue, purple,
green, orange, red lines denote level 1, 2, 3, 4, and 5 respectively. $t_{l}$
and $h_{l}$ are in Table 1. $\nu_{l}=2.5$ for every $l$, $d=1$. The solid
lines denote the numerical solutions of optimization problem in Proposition 2
and the dashed line denote the closed from in Remark 1.

##  5 Numerical examples

In this section, we present the results of numerical examples in comparison
with existing multi-fidelity kriging models. Emulators by MLASCE yields more
accurate results while the existing methods show relatively poor performance
even with the sufficient amount of training data due to the inherent bias. In
our algorithm, $a_{l}$ in (17) is set to $t_{l}$ since the difference between
successive computer simulations gets smaller and smaller. Section 5.1
highlights the difference in performance of MLASCE and other methodology in a
simple example. Section 5.2 shows that MLASCE is more robust in the case where
different versions of a computer simulation have different degrees of
smoothness.

###  5.1 A toy model

A simple toy model is considered to illustrate the performance of MLASCE and
contrast it with the models in [19, 18] as the most sophisticated formulation
of multi-fidelity GPR and strategy of sequential design. Three levels of
computer simulations defined in $[0,\pi]$ are used (Figure 3(a)). This type of
nested functions are used as the numerical examples in [17]. We define the
simulations $f_{1}(x)=\mathrm{sin}x$,
$f_{2}(x)=f_{1}(x)+\xi\bigl{(}x,\frac{\pi}{3},0.4\bigr{)}$ and
$f_{3}(x)=f_{2}(x)-\frac{1}{2}\xi\bigl{(}x,\frac{\pi}{4},0.2\bigr{)}+\frac{1}{2}\xi\bigl{(}x,\frac{3\pi}{4},0.2\bigr{)}$
where
$\xi(x,a,\lambda)=\bigl{(}1+\frac{\sqrt{5}}{\lambda}|x-a|+\frac{5|x-a|^{2}}{3\lambda^{2}}\bigr{)}\mathrm{exp}\bigl{(}-\frac{\sqrt{5}}{\lambda}|x-a|\bigr{)}$.

|
---|---

(a)

(b)

Figure 3: The graphs of computer simulations (left) and emulators (right).

$f_{1}(x)$, $f_{2}(x)$ and $f_{3}(x)$ share the same trend $\mathrm{sin}x$
while $f_{2}(x)$ and $f_{3}$ have the idiosyncratic movements of
$\xi(x,a,\lambda)$ at $x=\pi/3,\pi/4,3\pi/4$ respectively. A similar
phenomenon can be found in a simulation of complex fluid dynamics as the
degree of discretization of the numerical simulator improves. The
computational costs for running these simulations are illustrated in Table 2.
These costs are artificially assumed and increase at an exponential rate of
the levels. We consider five cases of different total budget $340$, $380$,
$420$, $460$ and $500$. As discussed in the next part where the methods of
[19, 18] are compared with ours, their formulations require sufficiently large
amount of training data for every level hence the generous budget while MLASCE
is robust in small size of training data thanks to relatively few numbers of
hyperparameters.

Table 2: Computational costs of computer simulations $f_{1}(x)$ | $f_{2}(x)$ | $f_{3}(x)$ | $h_{1}(x)$ | $h_{2}(x)$ | $h_{3}(x)$
---|---|---|---|---|---
4 | 16 | 64 | 4 | 20 | 80

We use the incremental functions $h_{l},l=1,2,3$ in Implicit method defined as
$h_{1}(x)=f_{1}(x)$ and $h_{l}(x)=f_{l}(x)-f_{l-1}(x)\hskip 4.0pt(l=2,3)$.

####  5.1.1 Comparison with other existing methods

We prepare three Gaussian processes
$GP(0,K_{5/2,\lambda_{l},\sigma^{2}_{l}}(x,x^{\prime}))$ with a Matérn kernel
for $l=1,2,3$. MICE and fitting of GPs are implemented by python modules
“Multi-Output Gaussian Process Emulator” and “scikit-learn” and a nugget
parameter is chosen to be $10^{-8}$ for every level. The mean predictions
$m^{f_{3}}_{3}\left(x\right)=\sum_{l=1}^{3}m_{l,N_{l}}^{h_{l}}\left(x\right)$
of Implicit method with different budgets are illustrated in Figure 3(b) in
comparison with [19, 18] and the numbers of runs are also presented in Figure
4(b). As expected, the performance of the emulators improves as the budget
increases. Since MLASCE tends to run the function which yields larger
incremental value $\delta_{l}(x)$ more often than the other, the number of
training data points at level $3$ gets larger if $\delta_{3}$ is affordable.
Moreover, the numbers of training data points at level $1$ and $2$ stabilize
after sufficient numbers of data points are obtained while MLASCE focus on
low-fidelity simulator if the budget is small. The $L^{2}$ error is also
presented in Figure 4(a) and $L^{2}$ error is defined as the integral of
squared difference between the mean prediction and $f_{3}(x)$ in $[0,\pi]$.

|
---|---

(a)

(b)

Figure 4: The $L^{2}$ error (left) and numbers of runs of different emulators
(right).

We compare MLASCE with the model proposed by [19, 18] by looking at the
performances of these GP emulators under the same amount of computational
budget. Following the formulation in (5), we prepare three Gaussian processes
$GP(\beta_{0,l},K_{5/2,\lambda_{l},\sigma^{2}_{l}}(x,x^{\prime}))$ with a
Matérn kernel for $\delta_{1}(x),\delta_{2}(x)$ and $\delta_{3}(x)$ in (5)
where $\beta_{0,l}$ denotes a constant mean. We specify the adjustment effect
$\rho_{l}(x)$ as
$\rho_{l}(x)=\beta_{\rho_{l},0}+\beta_{\rho_{l},1}x+\beta_{\rho_{l},2}x^{2}$
for the method by [19] where $\beta_{\rho_{l},0}$, $\beta_{\rho_{l},1}$ and
$\beta_{\rho_{l},2}$ are constant coefficients. This specification is to
include nonlinear adjustment effect since the relationships among different
$f_{i}(x)$ seem nonlinear. On the other hand, the method of [18] assumes
constant for the adjustment effect $\rho_{l}(x)=\rho_{l}$. $f_{l}(x)$ is used
as a simulation instead of an increment $h_{l}(x)$. We assume noninformative
priors for all of the hyperparameters and maximum likelihood estimation is
employed.

As proposed in [19], a strategy of nested space filling design is presented in
which the designs of input at higher levels are incorporated in the ones at
lower levels while the method of [18] employs their own strategy of sequential
design. Their strategy of design and construction of the emulators are
implemented by R package ”MuFiCokriging”
(https://CRAN.R-project.org/package=MuFiCokriging). Although a rule for
deciding the numbers of design points is not presented in [19], the number of
data points at each level is determined to keep the good proportion of the
data size across every level (Figure 4(b)). More crucially, due to the large
amount of hyperparameters in their formulation, sufficient numbers of training
data points should be provided even for the most expensive simulators.
Otherwise fitting the hyperparameters fails and the emulator cannot be
constructed. This drawback also holds for [18] and it requires decent numbers
of initial design points at every level before implementing their strategy of
sequential design even though it originally intends to allocate limited
computational resources among cheap and expensive simulators. On the other
hand, in MLASCE, only the two hyperparameters $\lambda_{l},\sigma_{l}^{2}$ are
estimated at each level and the incremental function $\delta_{l}$ tends to
take smaller values hence fitted hyperparameters are relatively more stable.
Although the method of [19] tries to overcome the nonlinear correlations among
simulators of different fidelity by designing $\rho_{l}(x)$, this attempt
requires large amount of training data, which is practically impossible if the
simulator is expensive. On the other hand, constant $\rho_{l}(x)=\rho_{l}$ by
[18] ends up yielding much worse emulators. As discussed previously, the
budget is decided so that their method can construct an emulator and more data
is available for them since the cheaper simulator $f_{l}(x)$ is used instead
of $h_{l}(x)$.

The results are shown in Figure 3(b). Although the mean predictions improve as
the budget increases, MLASCE performs better than the method of [19, 18]
especially in the cases of the smaller budgets. This is more evident in terms
of $L^{2}$ error (Figure 4(a)). The key difference is whether the bias is
considered or not. As already mentioned, the existing co-kriging methods do
not take into account the error of the emulator and only consider the
posterior variance as an indicator of accuracy. This leads to relatively poor
performance as an emulator even though enough amount of training data is
provided. The sequential design strategy by [18] may not be helpful since it
places too much emphasis on the cheaper simulations, which is the original
intention of theirs, hence skews the emulation. This type of problem of the
existing co-kriging methods becomes more evident or even worse in the next
part, where the computer simulations with more levels of fidelity are dealt
with.

###  5.2 More levels and variations of smoothness

We focus on the functions with different smoothness and compare the
performance of emulators with kernels of different smoothness. As discussed,
it is unlikely to know in advance the exact degree of smoothness of the
simulation in practice and an estimation of the smoothness parameter $\nu$ in
a Matérn kernel is not reliable if the simulation is expensive to run. Thus,
as mentioned before, a feasible strategy is to choose relatively small $\nu$
so that the corresponding RKHS can incorporate wider class of functions. The
simulated function can be less smooth as the degree of fidelity in the
simulator refines hence flexible choice of $\nu$ may improve the quality of
the emulator. In the following examples, we see that the kernels with fixed
smoothness across every level lead to relatively inefficient surrogate.
Besides, we also assume more levels of computer simulations so that our method
is robust in such a case.

The five levels of functions are defined in $[0,\pi]$ and five different bump
effects are incorporated: $f_{1}(x)=\mathrm{sin}x$,
$f_{2}(x)=f_{1}(x)+\xi_{2}\bigl{(}x,\frac{\pi}{6}\bigr{)}+\xi_{2}\bigl{(}x,\frac{5\pi}{6}\bigr{)}$,
$f_{3}(x)=f_{2}(x)-\xi_{3}\bigl{(}x,\frac{\pi}{4}\bigr{)}-\xi_{3}\bigl{(}x,\frac{3\pi}{4}\bigr{)}$,
$f_{4}(x)=f_{3}(x)+\xi_{4}\bigl{(}x,\frac{\pi}{3}\bigr{)}+\xi_{4}\bigl{(}x,\frac{2\pi}{3}\bigr{)}$,
and
$f_{5}(x)=f_{4}(x)+\xi_{5}\bigl{(}x,\frac{\pi}{8}\bigr{)}-\xi_{5}\bigl{(}x,\frac{4\pi}{8}\bigr{)}+\xi_{5}\bigl{(}x,\frac{7\pi}{8}\bigr{)}$
where
$\xi_{2}\bigl{(}x,a\bigr{)}=\mathrm{exp}\Bigl{(}-\frac{1}{(\pi/8)^{2}-(x-a)^{2}}\Bigr{)}\times\boldsymbol{1}_{(|x-a|<\frac{\pi}{8})}$,
$\xi_{3}\bigl{(}x,a\bigr{)}=0.3\hskip
2.0pt\mathrm{e}^{-8(x-a)^{2}}\left(-|x-a|^{5}+1\right)$,
$\xi_{4}\bigl{(}x,a\bigr{)}=0.3\hskip
2.0pt\mathrm{e}^{-8(x-a)^{2}}\left(-|x-a|^{3}+1\right)$ and
$\xi_{5}\bigl{(}x,a\bigr{)}=0.15\hskip
2.0pt\mathrm{e}^{-12(x-a)^{2}}(x-a+1)^{2}\boldsymbol{1}_{(a-1<x\leq
a)}+0.3\hskip 2.0pt\mathrm{e}^{-12(x-a)^{2}}\bigl{(}1-{1\over
2}(x-a-1)^{2}\bigr{)}\boldsymbol{1}_{(a<x\leq a+1)}+0.3\hskip
2.0pt\mathrm{e}^{-12(x-a)^{2}}\boldsymbol{1}_{(a+1<x)}$.
$\boldsymbol{1}_{(\cdot)}$ denotes an indicator function.

As a function of $x$ and for a fixed $a$, $\xi_{2}$ is a function of class
$C^{\infty}$ (but not analytical) and $\xi_{3}$, $\xi_{4}$, $\xi_{5}$ is of
class $C^{4}$, $C^{2}$ and $C^{1}$. Thus, every $f_{i}$ has the same degree of
smoothness as that of $\xi_{i}$ while $f_{1}$ is analytical. An equivalent
Sobolev space $H^{s}(\mathbb{R}^{d}))$ is incorporated in
$C^{[s-d/2]}(\mathbb{R}^{d})$ for nonnegative and non-integer $s-d/2$ where
$[s-n/2]$ denotes the maximum integer which does not exceed $s-n/2$. Then, we
can see that the equivalent space $H^{3}(\mathbb{R})\subset C^{2}(\mathbb{R})$
for $\nu=5/2$ in the Matérn kernel and $H^{2}(\mathbb{R})\subset
C^{1}(\mathbb{R})$ for $\nu=3/2$ albeit indirect relationships. These
functions are illustrated in Figure 5(a). The bump term $\xi_{2}$ is very
smooth and has small effect on the simulation while the other bump terms have
more effects (0.3 for $\xi_{3}$ and $\xi_{4}$ and 0.15 for $\xi_{5}$ at
location $x=a$) and become less smooth as levels go up. The computational
costs are listed in Table 3 and definition of $h_{l}$ is same as in the
previous part.

Table 3: Computational costs of computer simulations. The costs increase at an exponential rate of the levels. | $l=1$ | $l=2$ | $l=3$ | $l=4$ | $l=5$
---|---|---|---|---|---
$f_{l}(x)$ | 0.5 | 2 | 8 | 32 | 128
$h_{l}(x)$ | 0.5 | 2.5 | 10 | 40 | 160

|
---|---

(a)

(b)

Figure 5: The graphs of computer simulations (left) and emulators (right).

To implement MLASCE in the same way as before, we prepare five Gaussian
processes $GP(0,K_{\nu_{l},\lambda_{l},\sigma^{2}_{l}}(x,x^{\prime}))$ with a
Matérn kernel for $l=1,\ldots,5$. We set $\nu_{1}=3.5$, $\nu_{2},\nu_{3}=2.5$
and $\nu_{4},\nu_{5}=1.5$ based on the smoothness of each computer simulation.
We obtain the mean prediction by implementing MICE and fitting the data to
these GPs. On the other hand, other five Gaussian processes
$GP(0,K_{2.5,\lambda_{l},\sigma^{2}_{l}}(x,x^{\prime}))$ with a Matérn kernel
for $l=1,\ldots,5$ are used to prepare another mean prediction. The smoothness
parameter $\nu_{l}$ is chosen to be fixed in this case. For comparison, three
other existing methods ([19] with Gaussian kernels and Matérn kernels with
$\nu=2.5$ and [18] with Matérn kernels $\nu=2.5$) are implemented. The
Gaussian kernel is too smooth compared with the computer simulations and
Matérn kernels with $\nu=2.5$ may be still overly smooth for $f_{5}$. We
assume five levels for this model and the specification of adjustment effect
$\rho_{l}(x)$ is the same as in the previous part. Maximum Likelihood method
is used in fitting and nugget parameters are $10^{-8}$. We consider the
different case of computational budget $1050$, $1150$, $1250$, $1350$ and
$1450$ thereby sufficient amount of initial design is prepared for [19, 18].

The mean predictions with the budget $1150$ are illustrated in Figure 5(b).
MLASCE with varied smoothness captures well the rough movements of $\xi_{5}$,
$\xi_{4}$ and $\xi_{3}$ though $\xi_{2}$ is evaluated as the less important
factor due to its small effect (The numbers of samples are plotted in Figure
6(a)). MLASCE with Matérn kernel $\nu_{l}=2.5$ still mimics the overall
movements but gives the slightly less efficient result. The autoregressive co-
kriging [19] has difficulty in capturing non-smooth movements although
relatively large amount of samples is prepared for every level (Figure 6(a)).
The emulator of [19] with Gaussian kernels produces strange kinky movements
around the points where $\xi_{i}(x)$ lies and their method with Matérn kernel
$\nu_{l}=2.5$ also has difficulty in approximating the true function. Due to
the inefficient specification of the multi-fidelity modelling, the strategy by
[18] does not provide efficient emulators if the budget is limited. Moreover,
even if the budget is rich enough, numerical difficulty may emerge due to too
many samples from cheap simulations. This is because their strategy of
sequential design emphasizes low-fidelity simulations too much as the
improvements from simulators are directly weighted with corresponding
computational costs. In this setting, the computational costs increase at an
exponential rate (the most expensive simulation is 256 times more expensive
than the cheapest) hence applying their strategy directly could be
problematic.

These results show that an emulator by MLASCE with flexible choice of small
smoothness parameters performs better than the one with fixed and overly
smooth kernels. More importantly, autoregressive co-kriging with Gaussian
kernels may result in poor performance if the target is non-smooth. Besides,
even if the kernel is not so smooth, it may still have difficulty in mimicking
non-smooth behaviors. The comparison of these emulators in terms of $L^{2}$
error is illustrated in Figure 6(b).

|
---|---

(a)

(b)

Figure 6: The numbers of runs of different emulators (left) and $L^{2}$ error
(right).

##  6 Conclusions

We presented in this paper a novel method of multilevel adaptive sequential
design of computer experiments (MLASCE) in the framework of multi-fidelity
Gaussian process emulators. This strategy allows us to allocate efficiently
limited computational resource over simulations of different levels of
fidelity and build the GP emulator accordingly. Furthermore, the allocation of
computational resource is shown to be the solution of a simple optimization
problem in a special case that turns out to match well the general optimal
allocation of resources. This allows ambitious computer experiments at
multiple fidelities to be vetted before launching the expensive suite of
simulations, thus reassuring modellers that their budget will be effective in
building the multilevel emulator. The strength of the suggested approach is
that it considerably improves the accuracy of the prediction by a GP emulator
compared to other existing multi-fidelity models. One of the most important
consequences of our proposed method is that the size of improvements from
using more sophisticated, hence more expensive, simulations is evaluated while
the convergence of the prediction is guaranteed.

First, we decompose the most refined simulation into the functions of
increments and scales of these functions are measured in the RKHS norms.
Matérn kernels are used in the covariance structures of GPs for each increment
and a principle for choosing the smoothness parameters is suggested. Second,
using the RKHS norms of the conditional means of increments, we introduce
MLASCE, an algorithm of sequential designs for simulations of different levels
of fidelity based on MICE. Third, we present the explicit form of optimal
numbers of runs in the special cases where the designs are located in tensor
grids.

MLASCE is applied to toy examples in comparisons with other existing multi-
fidelity GP emulators. From this illustration we can derive two conclusions.
First, when the computational budget is small (i.e. the number of high
fidelity runs is tiny), other existing multi-fidelity methods do not perform
as well as MLASCE, or sometimes cannot even construct an emulator. Our method
is more robust whenever the budget is small, because of the increased
flexibility in adaptively choosing the level and the input locations.
Secondly, the choice of smoothness parameters of the GP emulators of our
method is adjustable. The two cases of series of either fixed or different
smoothness parameters are considered. In the usual case where the smoothness
of a simulation is unknown, MLASCE with varied levels of smoothness parameters
across different levels shows more accurate results. Indeed, our method is
able to follow the usual growth in roughness of differences across two
consecutive levels whenever fidelities increase. Obviously, this is now an ad
hoc approach to the variation in smoothness across levels, so more research,
e.g in terms of eliciting prior knowledge about smoothness across levels would
improve our method further.

## References

  * [1] J. Beck and S. Guillas.  Sequential design with mutual information for computer experiments (mice): emulation of a tsunami model.  SIAM/ASA Journal on Uncertainty Quantification, 4(1):739–766, 2016\.
  * [2] A. Berlinet and C. Thomas-Agnan.  Reproducing kernel Hilbert spaces in probability and statistics.  Springer Science & Business Media, 2011.
  * [3] K. M. Bretthauer and B. Shetty.  The nonlinear resource allocation problem.  Operations research, 43(4):670–683, 1995.
  * [4] K. M. Bretthauer and B. Shetty.  The nonlinear knapsack problem–algorithms and applications.  European Journal of Operational Research, 138(3):459–472, 2002\.
  * [5] D. A. Cohn, Z. Ghahramani, and M. I. Jordan.  Active learning with statistical models.  In Advances in neural information processing systems, pages 705–712, 1995.
  * [6] P. S. Craig, M. Goldstein, J. C. Rougier, and A. H. Seheult.  Bayesian forecasting for complex systems using computer simulators.  Journal of the American Statistical Association, 96(454):717–729, 2001.
  * [7] P. S. Craig, M. Goldstein, A. Seheult, and J. Smith.  Constructing partial prior specifications for models of complex physical systems: [read before the royal statistical society at a meeting on’elicitation ‘on wednesday, april 16th, 1997, the president, professor afm smithin the chair].  Journal of the Royal Statistical Society: Series D (The Statistician), 47(1):37–53, 1998.
  * [8] J. A. Cumming and M. Goldstein.  Small sample bayesian designs for complex high-dimensional models based on information gained using fast approximations.  Technometrics, 51(4):377–388, 2009.
  * [9] C. D’Ambrosio and S. Martello.  Heuristic algorithms for the general nonlinear separable knapsack problem.  Computers & operations research, 38(2):505–513, 2011.
  * [10] M. B. Giles.  Multilevel monte carlo methods.  Acta Numerica, 24:259, 2015.
  * [11] D. Higdon, M. Kennedy, J. C. Cavendish, J. A. Cafeo, and R. D. Ryne.  Combining field data and computer simulations for calibration and prediction.  SIAM Journal on Scientific Computing, 26(2):448–466, 2004.
  * [12] T. Hsing and R. Eubank.  Theoretical foundations of functional data analysis, with an introduction to linear operators, volume 997.  John Wiley & Sons, 2015.
  * [13] M. E. Johnson, L. M. Moore, and D. Ylvisaker.  Minimax and maximin distance designs.  Journal of statistical planning and inference, 26(2):131–148, 1990\.
  * [14] M. C. Kennedy and A. O’Hagan.  Predicting the output from a complex computer code when fast approximations are available.  Biometrika, 87(1):1–13, 2000.
  * [15] M. C. Kennedy and A. O’Hagan.  Bayesian calibration of computer models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(3):425–464, 2001.
  * [16] A. Krause, A. Singh, and C. Guestrin.  Near-optimal sensor placements in Gaussian processes: Theory, efficient algorithms and empirical studies.  Journal of Machine Learning Research, 9(Feb):235–284, 2008.
  * [17] L. Le Gratiet.  Bayesian analysis of hierarchical multifidelity codes.  SIAM/ASA Journal on Uncertainty Quantification, 1(1):244–269, 2013\.
  * [18] L. Le Gratiet and C. Cannamela.  Cokriging-based sequential design strategies using fast cross-validation techniques for multi-fidelity computer codes.  Technometrics, 57(3):418–427, 2015.
  * [19] L. Le Gratiet and J. Garnier.  Recursive co-kriging model for design of computer experiments with multiple levels of fidelity.  International Journal for Uncertainty Quantification, 4(5), 2014\.
  * [20] D. J. MacKay.  Information-based objective functions for active data selection.  Neural computation, 4(4):590–604, 1992.
  * [21] R. E. Marsten and T. L. Morin.  A hybrid approach to discrete mathematical programming.  Mathematical programming, 14(1):21–40, 1978.
  * [22] T. L. Morin and R. E. Marsten.  An algorithm for nonlinear knapsack problems.  Management Science, 22(10):1147–1158, 1976.
  * [23] B. Peherstorfer, K. Willcox, and M. Gunzburger.  Survey of multifidelity methods in uncertainty propagation, inference, and optimization.  Siam Review, 60(3):550–591, 2018.
  * [24] P. Perdikaris, D. Venturi, J. O. Royset, and G. E. Karniadakis.  Multi-fidelity modelling via recursive co-kriging and gaussian–markov random fields.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2179):20150018, 2015.
  * [25] P. Z. Qian and C. J. Wu.  Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments.  Technometrics, 50(2):192–204, 2008.
  * [26] C. S. Reese, A. G. Wilson, M. Hamada, H. F. Martz, and K. J. Ryan.  Integrated analysis of computer and physical experiments.  Technometrics, 46(2):153–164, 2004.
  * [27] J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn.  Design and analysis of computer experiments.  Statistical science, pages 409–423, 1989.
  * [28] T. J. Santner, B. J. Williams, W. I. Notz, and B. J. Williams.  The design and analysis of computer experiments, volume 1.  Springer, 2003.
  * [29] T. W. Simpson, D. K. Lin, and W. Chen.  Sampling strategies for computer experiments: design and analysis.  International Journal of Reliability and applications, 2(3):209–240, 2001.
  * [30] M. L. Stein.  Interpolation of spatial data: some theory for kriging.  Springer Science & Business Media, 2012.
  * [31] A. Stuart and A. Teckentrup.  Posterior consistency for Gaussian process approximations of bayesian posterior distributions.  Mathematics of Computation, 87(310):721–753, 2018.
  * [32] A. L. Teckentrup.  Convergence of Gaussian process regression with estimated hyper-parameters and applications in bayesian inverse problems.  SIAM/ASA Journal on Uncertainty Quantification, 8(4):1310–1337, 2020.
  * [33] R. Tuo and W. Wang.  Kriging prediction with isotropic matérn correlations: robustness and experimental designs.  Journal of Machine Learning Research, 21(187):1–38, 2020.
  * [34] W. Wang, R. Tuo, and C. Jeff Wu.  On prediction properties of kriging: Uniform error bounds and robustness.  Journal of the American Statistical Association, 115(530):920–930, 2020.
  * [35] H. Wendland.  Scattered data approximation, volume 17.  Cambridge university press, 2004.
  * [36] H. Wendland and C. Rieger.  Approximate interpolation with applications to selecting smoothing parameters.  Numerische Mathematik, 101(4):729–748, 2005.
  * [37] C. K. Williams and C. E. Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * [38] Z.-m. Wu and R. Schaback.  Local error estimates for radial basis function interpolation of scattered data.  IMA journal of Numerical Analysis, 13(1):13–27, 1993.
  * [39] G. Wynne, F.-X. Briol, and M. Girolami.  Convergence guarantees for Gaussian process approximations under several observation models.  arXiv preprint arXiv:2001.10818, 2020.
  * [40] X. Yang, D. Barajas-Solano, G. Tartakovsky, and A. M. Tartakovsky.  Physics-informed cokriging: A Gaussian-process-regression-based multifidelity method for data-model convergence.  Journal of Computational Physics, 395:410–431, 2019.
  * [41] X. Yang, X. Zhu, and J. Li.  When bifidelity meets cokriging: An efficient physics-informed multifidelity method.  SIAM Journal on Scientific Computing, 42(1):A220–A249, 2020.

##  Appendix A Proof of Theorem 2

###### Proof.

Let $L_{N}^{\theta}$ be a closure of
$\\{\sum_{i=1}^{N}a_{i}K_{\nu,\theta}(\cdot,x_{i});a_{i}\in\mathbb{R}\\}$,
that is

$L^{\theta}_{N}=\overline{\\{\sum_{i=1}^{N}a_{i}K_{\nu,\theta}(\cdot,x_{i});a_{i}\in\mathbb{R}\\}}.$

$\\{K_{\nu,\theta}(\cdot,x_{i})\\}_{i=1}^{N}$ are nonorthogonal basis in
$L^{\theta}_{N}$ and $L^{\theta}_{N}$ is a closed subspace of
$\mathbb{H}_{\theta}$. We obtain the orthogonal basis functions
$\\{e_{i}\\}_{i=1}^{N}$ in $L^{\theta}_{N}$ which are specified by

$(e_{1},\ldots,e_{N})^{\top}{}=K_{\nu,\theta}(\cdot,X_{N})K_{\nu,\theta}^{-1/2}(X_{N},X_{N}).$

It is seen that $\\{e_{i}\\}_{i=1}^{N}$ are orthogonal by taking the inner
product

$\langle
K_{\nu,\theta}(\cdot,X_{N})K_{\nu,\theta}^{-1/2}(X_{N},X_{N}),K_{\nu,\theta}(\cdot,X_{N})K_{\nu,\theta}^{-1/2}(X_{N},X_{N})\rangle_{\mathbb{H}_{\theta}}=I_{N}$

where $I_{N}$ denotes the $N$ dimensional identity matrix. Since
$L^{\theta}_{N}$ is a closed subspace of $\mathbb{H}_{\theta}$, the
decomposition $\mathbb{H}_{\theta}=L^{\theta}_{N}\oplus{L^{\theta}_{N}}^{c}$
holds. We define the projection operator $P_{N}$ onto $L^{\theta}_{N}$ for
$y\in\mathbb{H}_{\theta}$ as

$\displaystyle P_{N}y=\sum_{i=1}^{N}\langle
y,e_{i}\rangle_{\mathbb{H}_{\theta}}e_{i}=K_{\nu,\theta}(\cdot,X_{N})K_{\nu,\theta}^{-1}(X_{N},X_{N})\langle
y,K_{\nu,\theta}(\cdot,X_{N})\rangle_{\mathbb{H}_{\theta}}$

As $N\to\infty$, $L_{N}^{\theta}$ itself coincides with $\mathbb{H}_{\theta}$
because $\mathbb{H}_{\theta}$ is separable and
$\text{Span}\\{K_{\nu,\theta}(\cdot,x_{i})|1\leq i\leq\infty\\}$ is dense in
$\mathbb{H}_{\theta}$. As a result, $P_{N}$ tends to the identity map
$I:\mathbb{H}_{\theta}\to\mathbb{H}_{\theta}$ in the operator norm since

$\displaystyle\sup_{\|y\|_{\mathbb{H}_{\theta}}=1}\|(I-P_{N})y\|_{\mathbb{H}_{\theta}}=\sup_{\|y\|_{\mathbb{H}_{\theta}}=1}\|\sum_{i=N+1}^{\infty}\langle
y,e_{i}\rangle_{\mathbb{H}_{\theta}}e_{i}\|_{\mathbb{H}_{\theta}}=\sup_{\|y\|_{\mathbb{H}_{\theta}}=1}\Bigl{(}\sum_{i=N+1}^{\infty}\langle
y,e_{i}\rangle^{2}_{\mathbb{H}_{\theta}}\Bigr{)}^{1/2}\to 0$

as $N\to\infty$ where $\\{e_{i}\\}_{i=1}^{\infty}$ denote the orthogonal basis
functions in $\mathbb{H}_{\theta}$. Applying Mercer’s theorem (Theorem 4.6.5
in [12]), $K_{\nu,\theta}(x,u)$ is represented as
$K_{\nu,\theta}(x,u)=\sum_{i=1}^{\infty}\lambda_{i}\phi_{i}(x)\phi_{i}(u)$
with eigenvalues $\\{\lambda_{i}\\}_{i=1}^{\infty}$ and eigenfunctions
$\\{\phi_{i}\\}_{i=1}^{\infty}$. Then,
$K_{\nu,\theta}(x,X_{N})K_{\nu,\theta}^{-1}(X_{N},X_{N})K_{\nu,\theta}(X_{N},u)$
becomes

$\displaystyle
K_{\nu,\theta}(x,X_{N})K_{\nu,\theta}^{-1}(X_{N},X_{N})K_{\nu,\theta}(X_{N},u)$
$\displaystyle=$ $\displaystyle
K_{\nu,\theta}(x,X_{N})K_{\nu,\theta}^{-1}(X_{N},X_{N})\Bigl{(}\sum_{i=1}^{\infty}\lambda_{i}\phi_{i}(X_{N})\phi_{i}(u)\Bigr{)}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{\infty}\lambda_{i}\bigl{(}P_{N}\phi_{i}(x)\bigr{)}\phi_{i}(u)$

where $\phi_{i}(X_{N})=(\phi_{i}(x_{1}),\ldots,\phi_{i}(x_{N}))^{\top}{}$. We
can establish $P_{N}\phi_{i}\to\phi_{i}$ pointwise for every $i$ since
$P_{N}\to I$ in the operator norm as $N\to\infty$. Thus,

$\displaystyle
K_{\nu,\theta}(x,u)-K_{\nu,\theta}(x,X)K_{\nu,\theta}^{-1}(X,X)K_{\nu,\theta}(X,u)=\sum_{i=1}^{\infty}\lambda_{i}\phi_{i}(x)\phi_{i}(u)-\sum_{i=1}^{\infty}\lambda_{i}\bigl{(}P_{N}\phi_{i}(x)\bigr{)}\phi_{i}(u)\to
0.$

Similarly,

$\displaystyle
K_{\nu,\theta}(x,X_{N})K_{\nu,\theta}^{-1}(X_{N},X_{N})y(X_{N})$
$\displaystyle=$ $\displaystyle
K_{\nu,\theta}(x,X_{N})K_{\nu,\theta}^{-1}(X_{N},X_{N})\langle
y,K_{\nu,\theta}(\cdot,X_{N})\rangle_{\mathbb{H}_{\theta}}$ $\displaystyle=$
$\displaystyle P_{N}y(x)\to y(x)$

∎

######  Remark 2.

Assuming the same conditions in Theorem 2, let
$\\{\hat{\theta}_{N}\\}_{N=1}^{\infty}$ be a convergent sequence in $\Theta$
with a fixed limit $\theta_{0}\in\Theta$ and $g_{N}(\theta)\coloneqq
m_{N,\theta}^{y}\left(\boldsymbol{x}\right)$. We further assume that
$\sup_{\theta\in\Theta}\sum_{i=1}^{2}|\frac{\partial
g_{N}}{\partial\theta_{i}}(\theta)|<M$ for every $N$ and some finite $M$
independent of $N$ and $\\{X_{N},y(X_{N})\\}$ where $\theta_{1}=\lambda$ and
$\theta_{2}=\sigma^{2}$. Based on these above assumptions and for any
$\boldsymbol{x}\in\mathbb{D}$, we can establish
$|y(\boldsymbol{x})-m_{N,\hat{\theta}_{N}}^{y}\left(\boldsymbol{x}\right)|\to
0$. To verify this, first we see

$\displaystyle|y(\boldsymbol{x})-m_{N,\hat{\theta}_{N}}^{y}\left(\boldsymbol{x}\right)|\leq|m_{N,\theta_{0}}^{y}\left(\boldsymbol{x}\right)-m_{N,\hat{\theta}_{N}}^{y}\left(\boldsymbol{x}\right)|+|y(\boldsymbol{x})-m_{N,{\theta}_{0}}^{y}\left(\boldsymbol{x}\right)|$

and the second term converges to $0$ by applying Theorem 2. As for the first
term,

$|m_{N,\hat{\theta}_{N}}^{y}\left(\boldsymbol{x}\right)-m_{N,\theta_{0}}^{y}\left(\boldsymbol{x}\right)|=|g_{N}(\hat{\theta}_{N})-g_{N}(\theta_{0})|<M\|\hat{\theta}_{N}-\theta_{0}\|_{2}\to
0$

since $M$ is independent of $N$ and $\\{X_{N},y(X_{N})\\}$ and
$\hat{\theta}_{N}\to\theta_{0}$.

If $\hat{\theta}_{N}$ converges to $\theta_{0}$ in probability, the
convergence of the first term is also in probability. If $\hat{\theta}_{N}$ is
determined by maximum likelihood method, assuming the convergence in
probability is plausible.

The condition on $g_{N}(\theta)$ is partially verified by taking note that

$\displaystyle\frac{\partial g_{N}}{\partial\theta_{i}}(\theta)=$
$\displaystyle\Bigl{(}$ $\displaystyle\frac{\partial
k_{\nu,\theta}}{\partial\theta_{i}}(x,X_{N})K_{\nu,\theta}(X_{N},X_{N})^{-1}$
$\displaystyle-$ $\displaystyle
K_{\nu,\theta}\left(x,X_{N}\right)K_{\nu,\theta}(X_{N},X_{N})^{-1}\frac{\partial
k_{\nu,\theta}}{\partial\theta_{i}}(X_{N},X_{N})K_{\nu,\theta}(X_{N},X_{N})^{-1}\Bigr{)}y(X_{N})$

is continuous in $\Theta$ because of the form of Matérn kernel and the
assumption that $K_{\nu,\theta}(X_{N},X_{N})^{-1}$ exists hence
$\sup_{\theta\in\Theta}\sum_{i=1}^{2}|\frac{\partial
g_{N}}{\partial\theta_{i}}(\theta)|<M_{N,X_{N},y(X_{N})}$. Moreover, if
$\sup_{\theta\in\Theta}\sum_{i=1}^{2}|\frac{\partial
g_{N}}{\partial\theta_{i}}(\theta)|<M$ does not hold, it implies that
$\sup_{\theta\in\Theta}|\frac{\partial g_{N}}{\partial\theta_{i}}(\theta)|$ is
diverging as the number of samples grows.
$\sup_{\theta\in\Theta}\sum_{i=1}^{2}|\frac{\partial
g_{N}}{\partial\theta_{i}}(\theta)|<M$ may be usually implicitly assumed for
the numerical stability in implementing the hyperparameter estimation
typically maximum likelihood method with the gradient.