# Fixed-budget optimal designs for multi-fidelity computer experiments

Gecheng Chen and Rui Tuo

###### Abstract

This work focuses on the design of experiments of multi-fidelity computer
experiments. We consider the autoregressive Gaussian process model proposed by
Kennedy and O’Hagan (2000) and the optimal nested design that maximizes the
prediction accuracy subject to a budget constraint. An approximate solution is
identified through the idea of multi-level approximation and recent error
bounds of Gaussian process regression. The proposed (approximately) optimal
designs admit a simple analytical form. We prove that, to achieve the same
prediction accuracy, the proposed optimal multi-fidelity design requires much
lower computational cost than any single-fidelity design in the asymptotic
sense. Numerical studies confirm this theoretical assertion.

###### keywords:

Gaussian process models, autoregressive models, multi-level approximation,
uncertainty quantification

††journal: JSPI

organization=Wm Michael Barnes ’64 Department of Industrial & Systems
Engineering,addressline=Texas A&M University, city=College Station,
postcode=77845, state=Texas, country=U.S.

This work introduces a modified autoregressive model to formulate computer
code with infinite fidelity levels.

This work proposes a fixed precision optimal design for the modified
autoregressive models.

This work analyzes the upper bound of the total cost for the proposed MLGP
design and compares it with the single fidelity designs.

##  1 Introduction

With the advance of mathematical modeling and numerical computation, computer
simulations have been widely used to mimic complex physical systems. Compared
with physical experiments, computer models are often much less costly and more
flexible to run. Conceptually, each computer model is based on a mathematical
model in terms of, for instance, a set of partial differential equations. The
aim of a computer code run is to find a numerical solution to this
mathematical model. The design and analysis of computer experiments, an area
drawing increasing attention recently, aim to use statistical approaches to
collect and analyze computer simulation responses and facilitate decision-
making. We refer to the books by Santner et al. (2003); Gramacy (2020) for an
introduction to this field. One of the most important tasks in this area is to
build surrogate models of computer simulations. A surrogate model is a
statistical model constructed using the computer outputs over a selected set
of input sites. A point evaluation of a surrogate model is much less costly
than running the corresponding computer simulation code. Thus we can use the
surrogate model to explore and analyze the underlying computer response
function, which expedites the decision-making.

Many computer simulation codes allow users to specify an accuracy level of the
numerical responses, in terms of the degree of discretization, or the number
of iterations, etc. Generally, there is a trade-off between the accuracy and
the computational cost: a high-accuracy computer run is more costly than a
low-accuracy run (Qian and Wu, 2008; Xiong et al., 2013). Many statisticians
have considered the problem of integrating the computer responses from
multiple accuracy (also known as fidelity) levels to make predictions. The
first model was suggested by Kennedy and O’Hagan (2000), who proposed a
Gaussian-process-based autoregressive model to integrate the output data from
computer experiments with multiple fidelity levels. This model has been widely
used in real-world applications, e.g., simulation oil pressure (Kennedy and
O’Hagan, 2000), designing of linear cellular alloy (Xiong et al., 2013) and
prediction of band gaps (Pilania et al., 2017). Alternative models multi-
fidelity computer experiments include Bayesian hierarchical models (Qian and
Wu, 2008), deep Gaussian process model Cutajar et al. (2019), and
nonstationary Gaussian process models (Tuo et al., 2014). In this article, we
focus on the autoregressive models given the wide acceptance of these models.

Despite the widespread use of multi-fidelity computer experiments, there
remain some unanswered questions. First, it is unclear why one would prefer a
multi-fidelity computer experiment over a simple single-fidelity experiment. A
conventional explanation says that the low-fidelity data can help explore the
underlying response function quickly, and the high-fidelity data are then used
to refine the model. This assertion is usually confirmed by numerical studies,
but there is a lack of quantitative analysis of this relationship.

The second question is regarding the choice of the set of input sites for the
computer simulation, known as a design of experiment. Studying the design of
computer experiments is considered important because each computer run can be
costly, and a good design can help collect more information under the same
computational cost (Santner et al., 2003; Fang et al., 2005). The importance
of design of experiments is even more dramatic for multi-fidelity computer
experiments, because a poor design can even have an adverse effect on the
performance.

A number of designs for multi-fidelity computer experiments are proposed in
the literature. Qian (2009) proposed the nested Latin hypercube designs
(NLHDs). The design points of an NLHD on one fidelity level form a Latin
hypercube, and those from different fidelity levels are nested. He and Qian
(2011) proposed the nested orthogonal array-based Latin hypercube designs to
achieve stratification in both bivariate and univariate margins. Rennen et al.
(2010) constructed nested maximin Latin hypercube designs using a numerical
search algorithm. See also Qian (2009, 2012); Hwang et al. (2016); Qian et al.
(2006); Tuo et al. (2014); Xiong et al. (2013); He and Qian (2011); Ehara and
Guillas (2021); Sung et al. (2022) for more related works. Despite the rich
literature, none of these works addressed the aforementioned questions. Some
existing works provided theoretical justifications in terms of numerical
integration, e.g., He and Qian (2016), but it is still not clear why these
designs are suitable for the autoregressive models.

There is another class of approaches from the area of applied mathematics for
a related problem, known as multi-level Monte Carlo (MLMC) (Giles, 2008, 2015;
Heinrich, 2001). These methods are algorithms for computing expectations that
arise in stochastic simulations with different levels of accuracy. MLMC
reduces the computational cost of the traditional Monte Carlo method by taking
most samples with low accuracy, and only a few samples taken at high accuracy.
Multifidelity Monte Carlo (MFMC) (Ng and Willcox, 2014), based on the control
variable techniques for Monte Carlo simulation, resembles the MLMC method for
more general engineering models. More related works can be found in Haji-Ali
et al. (2016); Peherstorfer et al. (2016a, b); Geraci et al. (2017);
Teckentrup et al. (2015). Even though these multi-level based designs are
equipped with rigorous mathematical foundations, they are not based on
Gaussian process models so that they cannot be used for posterior-based
analysis and uncertainty quantification.

This work aims at underpinning the multi-fidelity computer experiments by
analyzing their predictive error theoretically. Inspired by the idea of multi-
level methods, we formulate a model framework for multi-fidelity computer
experiments. We then consider the fixed precision optimal design defined as
the design with the minimum total simulation cost subject to an accuracy
guarantee. Our main theoretical contribution is to obtain the order of
magnitude of the simulation cost as the required predictive error of the
autoregressive model goes to zero. The analysis is based on the error bounds
for Gaussian process regression provided by Tuo and Wang (2020). We prove that
a well-designed multi-fidelity experiment will require much lower simulation
cost in the order of magnitude to reach the same accuracy level, compared with
the best single-fidelity experiment. This result confirms the superiority of
multi-fidelity computer experiments. Based on the theory, we propose a novel
design scheme referred to as the multi-level Gaussian process (MLGP) designs.
Compared with the existing optimal design methods, the proposed method is easy
to implement because it does not require a numerical search. Numerical studies
confirm the superiority of the proposed method over existing alternatives.

The remaining part of this work is organized as follows. Section 2 reviews the
autoregressive models. In Section 3, we introduce the formulation of the
modified autoregressive model and in Section 4, we introduce the formulation
of the modified autoregressive model and propose the MLGP design method. In
Section 5, we give the detailed steps for implementing the MLGP method in
applications. In Section 6, we conduct some numerical studies, which show the
nice performance of the proposed MLGP. The conclusions and some discussion are
provided in Section 7. Technical proofs are presented in the Appendix.

##  2 Review on Gaussian process regression and autoregressive models

In this section, we review the Gaussian process regression and the
autoregressive models for multi-fidelity computer experiments.

###  2.1  Gaussian process regression

Let $Z$ be a Gaussian process on $\mathbb{R}^{d}$ with mean zero and
covariance function $r(x_{1},x_{2})$. We call $Z$ stationary (Santner et al.,
2003) if $r$ can be expressed as a function of $x_{1}-x_{2}$, i.e.,
$r(x_{1},x_{2})=\sigma^{2}\Phi(x_{1}-x_{2}),$ where $\sigma^{2}$ is the
variance, and $\Phi$ is a correlation function with $\Phi(0)=1$; otherwise, we
call it nonstationary. Given scattered evaluations
$(x_{1},Z(x_{1})),\ldots,(x_{n},Z(x_{n}))$, the Gaussian process regression
method (also known as kriging) estimates $Z(x)$ for any $x\in\mathbb{R}^{d}$
using the conditional expectation

$\displaystyle\hat{Z}(x):=E(Z(x)|Z(x_{1}),\ldots,Z(x_{n}))=a^{T}(x)R^{-1}Y,$
(1)

where $a(x):=(r(x-x_{1}),\ldots,r(x-x_{n}))^{T},R=(r(x_{j}-x_{k}))_{jk}$ and
$Y=(Z(x_{1}),\ldots,Z(x_{n}))^{T}$. In the area of computer experiments, (1)
is used to reconstruct the response function $Z(x)$, where
$\\{x_{1},\ldots,x_{n}\\}$ denotes the set of designs points,
$\\{Z(x_{1}),\ldots,Z(x_{n})\\}$ the corresponding response, and $x$ an
untried point.

###  2.2 Autoregressive model

Kennedy and O’Hagan (2000) proposed a model to fit multi-fidelity computer
experiments using Gaussian process regression. The main objective of this
model is to incorporate the computer outputs with different fidelity levels
and predict the response value of the computer code at an untried point of the
highest fidelity level.

Suppose we have $K$ levels of computer codes, with increasing accuracies,
denoted as $y_{1},\dots,y_{K}$. Kennedy and O’Hagan (2000) considered the
following autoregressive model to link each pair of consecutive $y_{i}$’s:

$\displaystyle y_{i}(x)=\rho_{i-1}y_{i-1}(x)+\delta_{i}(x),\quad for\quad
i=2,\dots,K,$ (2)

where $\rho_{i}$’s are known or unknown autoregressive coefficients, and
$y_{1}$ and $\delta_{t}$’s as mutually independent Gaussian processes for all
t.

This autoregressive model contains some hyper-parameters. The parameter
estimation is usually done using maximum likelihood (Forrester et al., 2007)
or Bayesian methods (Kennedy and O’Hagan, 2000). Both approaches involve
certain iterative algorithms and thus are expensive in computation. Here we
skip the details of the hyper-parameters estimation because they are
irrelevant to our main topic. We only consider the prediction for the computer
outputs at the highest fidelity given the hyper-parameters. The autoregressive
model has been widely used in real-world applications for multi-fidelity
computer experiments, e.g., simulation oil pressure (Kennedy and O’Hagan,
2000), designing of linear cellular alloy (Xiong et al., 2013) and prediction
of bandgaps (Pilania et al., 2017).

##  3 Autoregressive models for computer code with infinite fidelity levels

In this work, we focus on computer experiments with infinite fidelity levels.
There are many examples of this kind. For example, in finite element analysis,
the accuracy is governed by the mesh density, which, at least conceptually,
has infinite levels. Other examples include the step length in a finite
difference algorithm, the number of iterations in an iterative algorithm, etc.
Such a parameter is referred to as a “tuning parameter” in Tuo et al. (2014).

Let us use the finite element method (FEM) as an example to introduce our
intuition. FEM generally uses piecewise linear functions over prespecified
local patches to approximate general functions, such as the partial
differential equation (PDE) solutions. The number of local patches is governed
by the chosen mesh size. A typical multi-fidelity computer experiment uses a
sequence of consecutively refined meshes. Figure 1 shows the finite elements
corresponding to the first three mesh sizes in the set of refining meshes
$\\{1/2,1/4,1/8,...,1/2^{n},...\\}$ in a unit square.

Figure 1: Consecutively refined meshes. The number of the triangle local
patches quadruples as the mesh size is halved.

Here comes the question: what is a suitable configuration of autoregressive
models for the above FEM example which can potentially have infinite fidelity
levels? Clearly, we need to make some modifications to (2). The most obvious
one follows from the fact that there are infinite fidelity levels in the above
example. That is, conceptually, we should have computer response functions
$y_{i}$ for $i=1,\ldots,K,\ldots$. Next, we presume that the FEM solution
converges to the exact solution of the PDE as the mesh size goes to zero,
which is usually the case. Then the underlying quantity of interest is the
function $y_{\infty}:=\lim_{i\rightarrow\infty}y_{i}$. To guarantee the
existence of this limit, more constraints must be imposed in the
autoregressive model (2). First, we must have
$\lim_{i\rightarrow\infty}\rho_{i}=1$ as the existence of $y_{\infty}$
requires $\lim_{i\rightarrow\infty}(y_{i}-y_{i-1})=0$. In this work, we assume
$\rho_{i}=1$ for all $i$ for simplicity. Then we have the expression

$y_{\infty}=y_{1}+\sum_{i=2}^{\infty}\delta_{i}.$

Note that
$\text{Var}(y_{\infty})=\text{Var}(y_{1})+\sum_{i=2}^{\infty}\text{Var}(\delta_{i})$.
To ensure that $y_{\infty}$ has a finite variance, $\delta_{i}$ must have
decreasing variances. Here we model $\delta_{i}\sim
GP(0,\lambda^{2i}\sigma^{2}\Phi_{i})$ for some correlation function
$\Phi_{i}$, variance $\sigma^{2}$ and $0<\lambda<1$. The exponential decay in
the variance can be justified because the mesh size in Figure 1 also decays
exponentially as the refinement continues. Since we have not seen any data in
the experimental design stage, for simplicity, we assume that all the
$\delta_{i}$’s share the same correlation function $\Phi$.

We now define the modified version of the autoregressive model as,

$\displaystyle y_{i}(x)=y_{i-1}(x)+\delta_{i}(x),\quad for\quad
i=0,\dots,K,\dots,$ (3)

where $y_{-1}=0$, and $\delta_{i}\sim GP(0,\lambda^{2i}\sigma^{2}\Phi)$ are
mutually independent.

##  4 Multiple fidelity designs based on MLGP

It is hopeful that a multi-fidelity experiment can perform better than a
single-fidelity experiment, in the sense that the computational cost is lower
to achieve a fixed predictive accuracy. Otherwise, there is no value in using
a multi-fidelity design. In this section, the goal is to understand the best
prediction performance of a multi-fidelity experimental design.

First, we introduce an assumption for the computational cost of each
experimental trial. At the same time, the number of finite elements increases
exponentially as the mesh refines iteratively so that the resulting cost for
each untried point also increases exponentially. In view of this, we assume
the cost for each sample in the $i$-th level as $C_{b}a^{i}\leq C_{i}\leq
C_{a}a^{i},$ where $C_{a}$ and $C_{b}$ are constants, $a>1$ is the cost ratio.
Additionally, the difference in the approximation values given by successive
mesh sizes decreases as the mesh refines.

###  4.1 A fixed precision optimal design for multi-fidelity models

In this section, we introduce the problem setting for constructing optimal
designs for the modified autoregressive model with fixed precision. As
discussed in the above section, we use $y_{\infty}$ to denote the underlying
function. (3) allows us to reduce the multi-fidelity problem to a sequence of
single-fidelity problems. If we can reconstruct $\delta_{i}$ (i.e.,
$y_{i}-y_{i-1}$) for each $i$, we are able to reconstruct $y_{\infty}$ using
the summation of all these constructed functions. To reconstruct $\delta_{i}$,
a simple idea is to evaluate the function $\delta_{i}$ over a set of selected
points, denoted as $\mathcal{X}_{i}$ correspondingly, via the Gaussian process
regression. We assume the underlying function in the $K$-th layer, i.e.,
$y_{K}=\sum_{i=0}^{K}\delta_{i}$, is sufficiently accurate. One can
reconstruct the $y_{K}$ by $\hat{y}_{K}$ defined as (4) to approximate
$y_{\infty}$.

$\displaystyle\hat{y}_{K}(x)=E\left[\sum_{i=0}^{K}\delta_{i}(x)\Big{|}y_{K}(\mathcal{X}_{K}),...,y_{0}(\mathcal{X}_{0})\right]$
(4)

Denote the desired accuracy level, in terms of the lack of fit, as $\epsilon$
($\epsilon<1$). Then we can formulate the general objective function for
constructing optimal designs for the modified autoregressive model as

$\displaystyle\min_{n_{i},K}\quad$ $\displaystyle\sum_{i=0}^{K}n_{i}C_{i}$ (5)
s.t. $\displaystyle
E\left\|y_{\infty}-\hat{y}_{K}\right\|^{2}_{L^{2}}\leq\epsilon^{2},$ (6)

where $n_{i}\in\mathbb{N}$ denotes the sample size of $\mathcal{X}_{i}$, and
$\left\|y_{\infty}-\hat{y}_{K}\right\|^{2}_{L^{2}}$ is the expected squared
predictive error of $\hat{y}_{K}$. This objective function can be interpreted
as: we want to minimize the cost caused by the designs $\\{n_{i},K\\}$ which
satisfies the accuracy requirement.

###  4.2 MLGP method

The main difficulties for solving (5) include: 1. $y_{\infty}$ is unknown; 2.
$\hat{y}_{K}$ shown in (4) does not have a closed form under an arbitrary
design; 3. the integer optimization problem is NP-hard. In this section, we
make some approximations and restrictions to make (5) solvable. In detail, to
make the expected squared error in (6) tractable, we only consider nested
designs Kennedy and O’Hagan (2001); Qian (2009) to simplify $\hat{y}_{K}$.
Additionally, we replace the left-hand side of (6) with its upper bound to get
an approximated solution. We also relax $n_{i}$’s and $K$ to be real numbers
at first to search for the optimal solution, and then round down them to the
final optimal solution.

The widely used nested designs can refine the low-accuracy experiments with
the high-accuracy experiments to obtain a better model because the response is
also available in low-accuracy experiments at the same design points (Qian,
2009). Denote $\mathcal{X}_{0},\mathcal{X}_{1},\ldots,\mathcal{X}_{K}$ as a
sequence of nested sets of points, then the correlation between them can be
formulated as
$\mathcal{X}_{K}\subset\cdots\subset\mathcal{X}_{1}\subset\mathcal{X}_{0}$. A
good property of the nested design given in Theorem 1 helps simplify the
approximation (4).

######  Theorem 1.

Suppose mutually independent $\delta_{i}\sim GP(0,\sigma^{2}\Phi_{i})$,
$i=0,1,\dots,K$, where $\sigma^{2}$ is the variance, $\Phi$ is the correlation
function and $0<\lambda^{2}<1$. Let
$\mathcal{X}_{K}\subset\cdots\subset\mathcal{X}_{1}\subset\mathcal{X}_{0}$ be
a sequence of nested sets of points. Then for each $x$, the estimator
$\hat{y}_{K}(x)$ defined in (4) admits the expression

$\displaystyle\hat{y}_{K}(x)=\sum_{i=0}^{K}E\left[(y_{i}-y_{i-1})(x)\Big{|}(y_{i}-y_{i-1})(\mathcal{X}_{i})\right].$
(7)

Theorem 1 shows that if a nested design is used, Gaussian process regression
for an autoregressive model can be done by solving a set of single-fidelity
Gaussian process regression problems. Specifically, (7) is a sum of $K+1$ GP
predictors. Theorem 1 gives us a good incentive to use nested designs, and
then the multi-level interpolation approximation (4) can be rewritten as (7).

The next step is to identify the optimal nested designs for the (5). Because
$\hat{y}_{K}$ depends only on $y_{0},\ldots,y_{K}$, we can see that
$y_{\infty}-y_{K}$ is independent of $y_{K}-\hat{y}_{K}$ and therefore,

$\displaystyle
E\left\|y_{\infty}-\hat{y}_{K}\right\|^{2}_{L^{2}}=E\left\|y_{\infty}-y_{K}\right\|^{2}_{L^{2}}+E\left\|y_{K}-\hat{y}_{K}\right\|^{2}_{L^{2}}.$
(8)

To achieve the desired accuracy, a simple idea is to bound each of the right-
hand side terms of (8) by $\epsilon^{2}/2$, i.e., we convert the constrain (6)
into

$\displaystyle E\left\|y_{\infty}-y_{K}\right\|^{2}_{L^{2}}<\epsilon^{2}/2,$
(9)

and

$\displaystyle E\left\|y_{K}-\hat{y}_{K}\right\|^{2}_{L^{2}}<\epsilon^{2}/2.$
(10)

where (9) reflects the error caused by using $y_{K}$ to approximate
$y_{\infty}$ and (10) reflects the error caused by the Gaussian process
interpolation. Next, we will handle these two parts separately.

First, we focus on (9). Given the modified autoregressive model with shrinking
variance, it is easy to compute the upper bound of the error of the left-hand
side of (9).

######  Theorem 2.

Given mutually independent $\delta_{i}\sim GP(0,\lambda^{2i}\sigma^{2}\Phi)$,
$i=0,1,\dots,K,\dots$, with $0<\lambda^{2}<1$ and $\sigma^{2}>0$, we have

$\displaystyle E\left\|y_{\infty}-y_{K}\right\|_{L^{2}}^{2}\leq
C_{l}^{2}\lambda^{2K}\sigma^{2},$ (11)

where $C_{l}=\sqrt{C_{s}}/(1-\lambda)$ is a constant, and $C_{s}$ is the
volume of the region $\mathcal{X}$.

Based on Theorem 2, it is sufficient to bound the upper bound of the left-hand
side of (11) by $\epsilon^{2}/2$ to make (9) hold, that is,

$\displaystyle E\left\|y_{\infty}-y_{K}\right\|^{2}_{L^{2}}\leq
C_{l}^{2}\lambda^{2K}\sigma^{2}\leq\epsilon^{2}/2.$

Then we can choose

$\displaystyle
K=\left\lceil\log_{\lambda^{2}}\frac{\epsilon^{2}}{2C_{l}^{2}\sigma^{2}}\right\rceil,$
(12)

as the total number of levels $K$, where $\left\lceil a\right\rceil$ denotes
the largest integer which is smaller or equal to $a$.

Next, we turn to (10), which is associated with the error from the Gaussian
process modeling. Here we make a similar assumption as the Tuo and Wang (2020)
to our nested designs. Before introducing the assumption, we make some
definitions. For a set of design points
$\mathcal{X}=\\{x_{1},x_{2},...,x_{n}\\}\subset\Omega$, define the fill
distance as

$h_{\mathcal{X},\Omega}=\sup_{x\in\Omega}\inf_{x_{j}\in\mathcal{X}}||x-x_{j}||,$

and the separation radius as

$q_{\mathcal{X}}=\min_{1\leq j\neq k\leq n}||x_{j}-x_{k}||/2.$

######  Assumption 1.

$\mathcal{X}_{K}\subset...\subset\mathcal{X}_{1}\subset\mathcal{X}_{0}$ is a
sequence of nested sets of points, denote the corresponding fill distance and
separation radius as $h_{\mathcal{X}_{i},\Omega}$ and $q_{\mathcal{X}_{i}}$.
There exists a constant $A$ such that
$h_{\mathcal{X}_{i},\Omega}/q_{\mathcal{X}_{i}}\leq A$ holds for all $i$.

Given scattered evaluations $(x_{1},Z(x_{1})),\ldots,(x_{n},Z(x_{n}))$ where
the design is quasi-uniform, if we use the Matérn kernel to reconstruct the
underlying function $Z(x)$, then there exist two positive constants $p$ and
$q$ such that Tuo and Wang (2020):

$E\left\|Z(x)-\hat{Z}(x)\right\|^{2}_{L^{2}}<pn^{-2\nu/d},$ (13)

and

$\displaystyle E\left\|Z(x)-\hat{Z}(x)\right\|^{2}_{L^{2}}>qn^{-2\nu/d},$ (14)

where $\hat{Z}(x)$ is the reconstructed function, $\nu$ is the smoothness
parameter of the Matérn kernel and $d$ is the dimension of the design.

Based on the error bound of the original Gaussian process regression, we
propose the upper bound of the interpolation error in the modified
autoregressive model in Theorem 4.

######  Theorem 3.

Assume a infinite series $\delta_{i}\sim GP(0,\lambda^{2i}\sigma^{2}\Phi)$ and
they are mutually independent. Denote $y_{K}$ as the summation of the first
$(K+1)$ terms, and $\hat{y}_{K}$ defined in (7) as the emulation of $y_{K}$
based on the set of points $\mathcal{X}_{i}$ and Matérn kernel, then under
Assumption 1 we have

$\displaystyle
E\left\|y_{K}-\hat{y}_{K}\right\|^{2}_{L^{2}}\leq\sum_{i=0}^{K}p\lambda^{2i}\sigma^{2}n_{i}^{-2\nu/d},$
(15)

where $n_{i}$ is the number of points of $X_{i}$, $d$ is the number of
dimension, $p$ is a constant, $\nu$ is the smoothness parameter of $\Phi$.

Based on Theorem 4, it is sufficient to bound the upper bound of the left-hand
side of (16) by $\epsilon^{2}/2$ to make (10) hold. Then we can approximate
the solution of the objective function (5). Based on Tuo and Wang (2020), we
state the upper bound of the interpolation error in the modified
autoregressive model in Theorem 4.

######  Theorem 4.

Assume a infinite series $\delta_{i}\sim GP(0,\lambda^{2i}\sigma^{2}\Phi)$ and
they are mutually independent. Denote $y_{K}$ as the summation of the first
$(K+1)$ terms, and $\hat{y}_{K}$ defined in (7) as the emulation of $y_{K}$
based on the set of points $\mathcal{X}_{i}$ and Matérn kernel, then under
Assumption 1 we have

$\displaystyle
E\left\|y_{K}-\hat{y}_{K}\right\|^{2}_{L^{2}}\leq\sum_{i=0}^{K}p\lambda^{2i}\sigma^{2}n_{i}^{-2\nu/d},$
(16)

where $n_{i}$ is the number of points of $X_{i}$, $d$ is the number of
dimension, $p$ is a constant, $\nu$ is the smoothness parameter of $\Phi$.

Based on Theorem 4, it is sufficient to bound the upper bound of the left-hand
side of (16) by $\epsilon^{2}/2$ to make (10) hold. Then we can approximate
the solution of objective function (5) with (6) by solving a relaxed version,
i.e.,

$\begin{split}\min_{n_{i}}\quad&\sum_{i=0}^{K}n_{i}C_{i}\\\
\mbox{s.t.}\quad&\sum_{i=0}^{K}p\lambda^{2i}\sigma^{2}n_{i}^{-2\nu/d}\leq\epsilon^{2}/2\end{split}$
(17)

where $K=\lceil 2\log_{\lambda^{2}}\frac{\epsilon}{2C_{l}\sigma}\rceil$. One
can see that (17) is an integer programming problem which is NP-hard. In view
of this, we consider a relaxed version of it. First we assume $n_{i}$’s can
choose any real number, then (17) becomes a convex optimization problem whose
optimal solution is

$\displaystyle
n_{i}=\left(\frac{\epsilon^{2}}{2p\sigma^{2}S}\right)^{-\frac{d}{2\nu}}\left(\frac{a^{i}}{\lambda^{2i}}\right)^{-\frac{d}{d+2\nu}},\quad
i=0,...,K,$ (18)

where
$S=\sum_{i=0}^{K}(a^{i})^{\frac{\nu}{d+2\nu}}(\lambda^{2i})^{\frac{d}{d+2\nu}}$,
$p$ is a constant. The proof of (18) is given in the Appendix.

Note that $\left\lceil n_{i}\right\rceil$’s also locate in the feasible region
of (17). We regard the $\left\lceil n_{i}\right\rceil$’s shown in (19) as a
nearly optimal solution to (17).

$\displaystyle
n_{i}=\left\lceil\left(\frac{\epsilon^{2}}{2p\sigma^{2}S}\right)^{-\frac{d}{2\nu}}\left(\frac{a^{i}}{\lambda^{2i}}\right)^{-\frac{d}{d+2\nu}}\right\rceil,\quad
i=0,...,K,$ (19)

In summary, we consider the designs generated based on (19) and (12) as the
optimal designs for the modified autoregressive model. We will discuss the
implementation details in Section 5.

###  4.3 Comparison with single-fidelity experiments

In this section, we will show that the MLGP needs less cost than any single-
level method to achieve a given accuracy provided that the accuracy level is
small enough. Let us denote $a=\exp{(\alpha)}$ and
$\lambda^{2}=\exp{(-\beta)}$ for convenience ($\alpha,\beta>0$). In this
section, we use the symbol “$\lesssim$” to represent that the left-hand side
is less than the right-hand side multiplied by a constant depending only on
$\Phi$. First, we give an upper bound for the total cost for the MLGP method
that reaches a fixed accuracy level $\epsilon$.

######  Theorem 5.

Assume $y_{i}(x)=y_{i-1}(x)+\delta_{i}(x)$, for $i=0,\dots,K,\dots$, where
$\delta_{i}{\sim}GP(0,\lambda^{2i}\sigma^{2}\Phi)$ and they are mutually
independent, $\lambda^{2}<1$ and $y_{-1}=0$. If we use the design (12) and
(19) to reach a fixed accuracy level $\epsilon$, then the upper bound for the
total cost is summarized as

$\displaystyle
C_{\epsilon}\lesssim\left\\{\begin{matrix}\epsilon^{-d/\nu},&2\alpha\nu<d\beta\\\
\epsilon^{-d/\nu}|\ln\epsilon^{2}|^{\frac{d+2\nu}{2\nu}},&2\alpha\nu=d\beta\\\
\epsilon^{-2\alpha/\beta},&2\alpha\nu>d\beta\end{matrix}\right.$ (20)

Next, we discuss the lower bound for the total cost if we implement the
single-fidelity method to reach the accuracy level $\epsilon$.

######  Theorem 6.

Assume $y_{i}(x)=y_{i-1}(x)+\delta_{i}(x)$, for $i=0,\dots,K,\dots$, where
$\delta_{i}{\sim}GP(0,\lambda^{2i}\sigma^{2}\Phi)$ and they are mutually
independent, $\lambda^{2}<1$ and $y_{-1}=0$. If we only spread design samples
in one layer to reach a fixed accuracy level $\epsilon$, then the lower bound
for the total cost is

$\displaystyle
C_{\epsilon}^{SL}\gtrsim\epsilon^{-(\frac{d}{\nu}+\frac{2\alpha}{\beta})}.$
(21)

##  5 Fixed budget designs

For many real applications, a more realistic assumption is that the total
budget, rather than the accuracy level, is fixed. In this section, we will
introduce how to implement the MLGP design method given a fix budget.

First, assume we know the values for the hyper-parameters $\lambda^{2}$,
correlation function, and $\nu$. In practice, $\lambda^{2}$ can be decided
based on expert knowledge. For example, we can choose the decreasing ratio of
the mesh sizes as $\lambda^{2}$ in the FEM case. The correlation function and
the corresponding parameter(s) need to be decided based on the smoothness of
the underlying function, see Wilson and Adams (2013); Schulz et al. (2018) for
more information.

In order to obtain the designs based on (12) and (19) under a certain budget
$B$, we need to have the exact values for the parameters $p$, $C_{l}$,
$\sigma$ and the desired accuracy level $\epsilon$. These parameters are often
hard to decide in real-world applications. However, the number of available
layers is often determined in real-world applications. Based on this, we just
assume $K$ is determined. Please note that this assumption does not mean we
need to use all the layers in our designs because $n_{i}$ can be very small
for a certain layer if the MLGP method thinks that the corresponding layer is
not needed.

The next thing is to calculate $n_{i}$. Denote
$\tau:=p\sigma^{2}/\epsilon^{2}$ and $\eta:=p/C_{l}^{2}$, then we rewrite (19)
as

$\displaystyle n_{i}=\left\lceil(2\tau
S)^{\frac{d}{2\nu}}(\frac{a^{i}}{\lambda^{2i}})^{-\frac{d}{d+2\nu}}\right\rceil,\quad
i=0,1,...,K$ (22)

where $\tau=\frac{\eta}{2\lambda^{2K}}$.

In order to avoid the estimation of those constants ($p$, $C_{l}$, $\sigma$
and $\epsilon$), one can tune $\eta$ to get the corresponding optimal design
under the fixed budget $B$. By choosing a value for $\eta$, we can get a value
for $\tau$ which leads to corresponding $n_{i},i=0,...,K$ according to (22).
The resulting design that triggers the highest cost and meets the budget
requirement can be chosen as the optimal design $n_{i}^{*}$’s. In
applications, the user can try different values for $\eta$ in a pre-specified
set (e.g., the integers from $1$ to $100$) to search for the optimal design
$n_{i}^{*}$’s.

Note that the design $n_{i}^{*}$’s does not use all the budget in most cases
because we need to choose an integer for $n_{i}^{*}$’s when implementing (22).
Denote the unallocated budget as

$\displaystyle B_{K}=B-\sum_{i=0}^{K}n_{i}^{*}C_{i}.$ (23)

The final step is to re-allocate $B_{K}$ to modify the design
$n^{*}_{i},i=0,...,K$. We implement a greedy algorithm to solve the following
objective function (24) to reallocate the budget to further reduce the error
in the left-hand side of the constraint of (17).

$\displaystyle\min_{n_{i}^{**}}\quad$
$\displaystyle\sum_{i=0}^{K}\lambda^{2i}\sigma^{2}((n_{i}^{**})^{-2\nu/d}-(n_{i}^{*})^{-2\nu/d})$
(24) s.t. $\displaystyle\sum_{i=0}^{K}(n_{i}^{**}-n_{i}^{*})C_{i}\leq B_{K}$
(25)

Then we get the final optimal designs $n^{**}_{i},i=0,...,K$.

Denote a $K$ level design as a vector $N=[n_{1},...,n_{K}]$, and a function
$G(N^{**})=\sum_{i=0}^{K}\lambda^{2i}\sigma^{2}((n_{i}^{**})^{-2\nu/d}-(n_{i}^{*})^{-2\nu/d})$.
Also denote a vector with the $i$-th element equal to $1$ and other elements
equal to $0$ (i.e., $1^{i}=[0,...,1,...0]$) as one sample in the $i$-th layer.
Now the detailed steps of the re-allocating starting from the design $N^{*}$
are as follows:

  * 1.

Search $i$ exhaustively in the range $[0,K]$ to minimize $G(N^{*}+1^{i})$ on
the condition that $C_{i}\leq B_{K}$;

  * 2.

Add one sample to the $i$-th layer of $N^{*}$, that is, $N^{*}=N^{*}+1^{i}$;

  * 3.

Calculate the unallocated budget $B_{K}=B_{K}-C_{i}$;

  * 4.

Repeat the above steps until $B_{K}=0$ or $B_{K}$ cannot afford any additional
sample.

######  Remark 1.

Before reaching the optimal designs, we need to determine the values for
several parameters, including $a$, $C_{a}$, $C_{b}$, $B$, $K$, $\lambda$,
$\nu$ and $\eta$. Some of the parameters are given or estimated by specific
applications, such as $a$, $C_{b}$, $C_{a}$, $B$ and $K$. Some of the
parameters need to be chosen by the users based on experience or other
knowledge, i.e., $\lambda$ and $\nu$. We may misspecify them because we do not
know the true values for these two parameters. However, empirical studies show
that the designs generated by the proposed MLGP can still outperform existing
alternatives even with misspecified values. See Section 6.2 for more details.
As for the tuning of $\eta$, screening only over integer values should give
decent tuning results. Note that we still have a re-allocation step after this
search that will further modify the design, so we still have a good chance to
get a nice design even if we do not reach the optimal $\eta$ exactly. If need,
when the optimal $\eta$ lies between two consecutive integers, we can consider
a finer discretization over the search space to obtain a better $\eta$.
Besides, some other searching methods are available to find the optimal $\eta$
more efficiently, such as the binary search algorithm which is applied in our
experiments.

Given the detailed steps for calculating the $n_{i}^{**}$’s for each layer, we
still need to decide how to spread these samples in each layer. Please recall
that the error bound in Theorem 4 only holds when Assumption 1 holds. But it
is very hard to maintain this assumption in applications rigorously because it
involves a lot of optimization steps. we consider approximately space-filling
sequences, such as Halton sequences (Halton, 1964) which is a type of low-
discrepancy sequences. Halton sequences enjoy the following two advantages.
First, it can give designs of any sample sizes. Second, the sequences have
analytical expressions which are easy to compute.

###  5.1 Examples

In this section, we visualize the designs generated by the MLGP under
different budgets or different $\lambda^{2}$ in $[0,1]^{2}$. We choose $a=4$
and $C_{a}=C_{b}=1$ for the experiments in this section.

In the first experiment, we use $B=192$, $K=2$ and Matérn correlation function
with $\nu=1.25$ to generate designs under four different $\lambda^{2}$ using
the MLGP. Tabel 5.1 and Figure 2 show the design structures when $\lambda^{2}$
changes from $1/8$ to $3/4$. Note that the column ”Design structure” shows the
number of samples from the lowest to the highest layer. For example, the
design $20,7,3$ means a nested design with $20$ samples (empty circles in
Figure 2) in the $0$-th layer, $7$ samples (filled circles in Figure 2) in the
$1$st layer and $3$ samples (filled triangles in Figure 2) in the $2$nd layer.
The subplot a, b, c and d are corresponding to $\lambda^{2}=1/8$,
$\lambda^{2}=1/4$, $\lambda^{2}=1/2$ and $\lambda^{2}=3/4$, respectively. It
can be seen that the MLGP allocates more budget on the higher layers with the
$\lambda^{2}$ increases. This phenomenon exactly matches the changing of the
true function. This is because the variance of $\delta_{i}$’s in (7) in the
high layers take less proportion in the underlying function $y_{K}$ if
$\lambda^{2}$ is small. In other words, the information of the higher layers
plays a less important role in the overall accuracy. As a result, if we want
to get a low error as much as possible under a fixed budget, we should pay
more budget on the samples in the lower layers and vice versa. This exclusive
flexibility of the MLGP can help the users a lot to construct suitable designs
when facing different $\lambda^{2}$ while the multi-level and single-level
methods can only build designs with fixed proportions.

In the second experiment, we use $K=2$, $\lambda^{2}=1/2$, Matérn correlation
function with $\nu=1.25$ to generate designs under budgets from $96$ to $240$
using the MLGP. Table 5.1 and Figure 3 show the design structures generated by
the MLGP under different $B$.

In the third experiment, we test our method under large $\nu$’s for the Matérn
correlation function. We use $K=3$, $\lambda=1/2$ and budget $B=4760$ to
generate designs under different $\nu$’s from $5$ to $50$. Table 5.1 shows the
designs in this experiment, we can see that our method still works well under
large $\nu$’s.

Table 1: Designs under different $\lambda^{2}$ $\lambda^{2}$ | Design structure
---|---
$1/8$ | 48, 16, 5
$1/4$ | 32, 16, 6
$1/2$ | 24, 14, 7
$3/4$ | 16, 12, 8
Table 2: Designs under different budgets $B$ $B$ | Design structure
---|---
$96$ | 20, 7, 3
$144$ | 24, 10, 5
$192$ | 24, 14, 7
$240$ | 28, 17, 9
Table 3: Designs under different smoothness parameters $\nu$ $\nu$ | Design structure
---|---
$5$ | 24,16,16,7
$10$ | 16, 9, 17, 7
$25$ | 16, 9, 9, 8
$50$ | 16, 9, 9, 8

Figure 2: Designs generated by MLGP under different $\lambda^{2}$ Figure 3:
Designs generated by MLGP under different budgets $B$

##  6 Numerical study

In this section, some numerical studies will be conducted to show the
superiority of the proposed MLGP. In Section 6.1, we will first compare the
performance of the proposed MLGP, classical multi-fidelity designs, and
single-level designs in 1-dimensional cases when we know all hyper-parameters
($\lambda^{2}$, correlation function, and $\nu$). Then we will test the
performance of these three methods if one of the hyper-parameters is
misspecified in Section 6.2. In Section 6.4, we will compare the proposed MLGP
with Nested Latin Hypercube designs and single-level designs in a
2-dimensional case. Emulators given by the proposed MLGP yield more accurate
results than other methods in all cases.

###  6.1 Numerical study in one-dimensional space

In the section, we prepare a four-level function $f(x)$ (i.e., $K=3$) defined
in $[0,15]$ as the simulation function. Each level is assumed as a mutually
independent Gaussian process with zero mean and Matérn correlation function.
We assume $\lambda^{2}=1/3$, the correlation matrix is Matérn correlation
function with $\nu=1.25$. In this experiment, we choose the cost ratio $a=8$
and the cost for each sample in the $0$-th layer to be $1$, i.e.,
$C_{a}=C_{b}=1$.

In this experiment, we assume all the true values hyper-parameters are known.
We calculate three groups of designs based on the proposed MLGP method, the
classical multi-fidelity design method, and the single-level design method
under different budgets. After calculating the number of samples in each
layer, we generate the corresponding (nested) design samples based on the
Halton sequences for all three methods. Then we train a Gaussian process
regression (GPR) emulation based on each design to test the prediction
accuracy on a testing set with $200$ samples. Note that we do not use maximum
likelihood estimation (MLE) to estimate the variance and the length parameter.
The accuracy of GPR emulations is measured by the root mean squared error
(RMSE) defined as $RMSE=\sqrt{N^{-1}\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}},$
where $N=200$, $y_{i}$, and $\hat{y}_{i}$ are the true label and the
predictive label for the $i$-th testing sample, respectively.

In order to make the results more convincing and show the robustness of the
proposed MLGP algorithm, we repeat this experiment $30$ times starting from
different realizations of $f(x)$ and compare the mean RMSE of these methods.

Figure 4: RMSEs of the MLGP, classical multi-fidelity method and single-level
method for the 1-dimensional case under different budgets.

Figure 4 shows the average RMSEs of the three methods under different budgets.
The blue line shows the performance of the single-level method, and it almost
keeps the same even with different budgets. This is because in all cases, the
best accuracy of the single-level method is achieved by only sampling at the
$2$-th layer. The red line and the black line represent the performance of the
multi-level method and the MLGP, respectively. One can easily see that the
designs given by the MLGP lead to the emulations with the lowest RMSEs
compared with other methods, which proves the superiority of the proposed
MLGP.

###  6.2 Numberical study for the situations with misspecified hyper-
parameters

We may not get the exact values for all the hyper-parameters in some cases. To
show the robustness of the MLGP, we also test the performance of the proposed
MLGP when we incorporate a misspecified hyper-parameter. First, we use the
same parameters as the first experiment to generate the four-level function
$f(x)$, i.e., $K=3$, $\lambda^{2}=1/3$, Matérn correlation function with
$\nu=1.25$. Also, we choose the cost ratio $a=8$ and $C_{a}=C_{b}=1$. We
generate designs based on the proposed MLGP using all the same parameters but
a misspecified $\lambda^{2}=1/2$. Then we test the GPR emulation trained by
this design and compare the accuracy of the testing set with other methods.
Next, we use the hyper-parameters $K=3$, $\lambda^{2}=1/3$, Matérn correlation
function with $\nu=1.25$ to generate the four-level function $f(x)$ and
generate designs based on the proposed MLGP using all the same parameters but
a misspecified $\nu=1.5$. GPR emulation is also used to compare the
performance of the MLGP and the other two methods.

Table 4 shows the results for these two experiments. One can see that the MLGP
still gets the lowest error even though the misspecified hyper-parameters are
incorporated.

Table 4: Mean RMSEs for the proposed MLGP method, classical multi-fidelity design method and single-level design method with misspecified $\lambda^{2}$ or $\nu$ Situation | Method | Average RMSE | Designs structure
---|---|---|---
| MLGP | 0.2831 | 88, 32, 21, 6
misspecified $\lambda^{2}$ | Multi-fidelity | 0.3441 | 56, 28, 14, 7
| Single-level | 0.3789 | 0, 0, 74, 0
| MLGP | 0.1652 | 88, 40, 20, 6
misspecified $\nu$ | Multi-fidelity | 0.2127 | 56, 28, 14, 7
| Single level 2 | 0.2064 | 0, 0, 74, 0

###  6.3 Case study in two-dimensional space

In this section, we will compare the performance of the proposed MLGP, the
nested Latin hypercube design (NLHD), and the single-level design in a two-
dimensional case. We prepare a three-level function $f(x)$ (i.e., $K=2$)
defined in $[0,1]^{2}$ as the simulation function. Each level is assumed as a
mutually independent Gaussian process (GP) with zero mean and Matérn
correlation function with $\nu=1.25$. We assume $\lambda^{2}=1/2$. In this
experiment, we choose the cost ratio $a=4$ and the cost for each sample in the
$0$-th layer to be $1$, i.e., $C_{a}=C_{b}=1$.

In this experiment, we assume all the true parameters are known, i.e.,
$\lambda^{2}=1/2$, the correlation function is Matérn correlation function
with $\nu=1.25$. For the proposed MLGP and single-level methods, we first
calculate the design structure and then generate the corresponding (nested)
samples using the Halton sequences. As for the NLHD, we use a design given in
Table 7. Then we train a GPR emulation based on each design to test the
prediction accuracy on a testing set with $500$ samples based on the root mean
squared error (RMSE). In this experiment, $30$ different realizations of
$f(x)$ are generated to test the average RMSE of the GPR emulations. Table 5
shows the average RMSEs and the design structures of these three methods. We
can see that the MLGP yields the lowest RMSE, which demonstrates the
effectiveness of the MLGP in two-dimensional cases.

Table 5: Results for the proposed MLGP method, Nested Latin hypercube design (NLHD) and single level design method Method | Average RMSE | Designs structure
---|---|---
MLGP | 0.1510 | 20, 7, 3
Multi-fidelity | 0.1735 | 32, 8, 2
Single-level | 0.1658 | 0, 0, 6
Table 6: Nested Latin hypercube designs in two-dimensional space Level 3 |  (0.166670,0.863640), (0.712120,0.318180),
---|---
Level 2 |  (0.166670,0.863640), (0.712120,0.318180), (0.590910,0.712120), (0.984850,0.590910), (0.863640,0.984850), (0.045455,0.166670), (0.439390,0.045455), (0.287880,0.469700),
Level 3 |  (0.166670,0.863640), (0.712120,0.318180), (0.590910,0.712120), (0.984850,0.590910), (0.863640,0.984850), (0.045455,0.166670), (0.439390,0.045455), (0.287880,0.469700), (0.954550,0.409090), (0.833330,0.681820), (0.196970,0.136360), (0.348480,0.924240), (0.681820,0.954550), (0.651520,0.530300), (0.136360,0.378790), (0.530300,0.439390), (0.893940,0.833330), (0.621210,0.075758), (0.318180,0.772730), (0.257580,0.287880), (0.772730,0.106060), (0.075758,0.742420), (0.378790,0.196970), (0.106060,0.560610), (0.560610,0.227270), (0.803030,0.500000), (0.469700,0.621210), (0.227270,0.651520), (0.924240,0.257580), (0.409090,0.348480), (0.500000,0.893940), (0.742420,0.803030),

###  6.4 Case study in higher dimensional space

In this section, we will compare the performance of the proposed MLGP and the
single-level design in four-dimensional and eight-dimensional cases. We
prepare a three-level function $f(x)$ (i.e., $K=2$) defined in $[0,1]^{2}$ as
the simulation function. Each level is assumed as an i.i,d Gaussian process
(GP) with zero mean and Matérn correlation function with $\nu=1.25$. We assume
$\lambda^{2}=1/2$. In this experiment, we choose the cost ratio $a=4$ and the
cost for each sample in the $0$-th layer to be $1$, i.e., $C_{a}=C_{b}=1$.

In this experiment, we assume the available budget $B=144$, $\lambda^{2}=1/2$
and the correlation function is Matérn correlation function with $\nu=1.25$.
For the proposed MLGP and single-level methods, we first calculate the design
structure and then generate the corresponding (nested) samples using the
Halton sequences. Then we train a GPR emulation based on each design to test
the prediction accuracy on a testing set with $500$ samples based on the root
mean squared error (RMSE). In this experiment, $30$ different realizations of
$f(x)$ are generated to test the average RMSE of the GPR emulations. Table 7
shows the average RMSEs and the design structures of these two methods. We can
see that the MLGP yields the lowest RMSE, which demonstrates the effectiveness
of the MLGP in higher dimensional cases.

Table 7: Results for the proposed MLGP method and the single-level design method in higher dimensional cases | $d=4$ | $d=8$
---|---|---
| Avg. RMSE | Structure | Avg. RMSE | Structure
MLGP | 0.3702 | 24, 10, 5 | 1.074 | 12, 9, 6
Single-level | 0.4621 | 0, 0, 9 | 1.181 | 0, 0, 9

##  7 Conclusion

In this work, we propose new experimental designs for multi-fidelity computer
experiments, called the multi-level Gaussian process (MLGP) designs, and prove
that they are theoretically superior to any single-fidelity designs. In our
theoretical analysis, we assume that a quasi-uniform nested design is used. It
is not clear how to obtain such a nested design efficiently. Chen and Xiong
(2017) proposed a numerical algorithm to search for good nested designs. In
our work, we use Halton sequences, which are computationally efficient and
nearly quasi-uniform. The proposed designs are readily applicable for a wide
range of multi-fidelity computer experiments such as the heat exchanger
problem described in (Qian et al., 2006).

## Acknowledgement

The authors are grateful for the AE and reviewers for very helpful comments.
This research is supported by NSF DMS-1914636, DMS-2312173, and CNS-2328395.

## References

  * Santner et al. (2003) T. Santner, B. Williams, W. Notz, The Design and Analysis of Computer Experiments, Springer Verlag, 2003\.
  * Gramacy (2020) R. B. Gramacy, Surrogates: Gaussian process modeling, design, and optimization for the applied sciences, Chapman and Hall/CRC, 2020.
  * Qian and Wu (2008) P. Z. Qian, C. J. Wu,  Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments,  Technometrics 50 (2008) 192–204.
  * Xiong et al. (2013) S. Xiong, P. Z. Qian, C. J. Wu,  Sequential design and analysis of high-accuracy and low-accuracy computer codes,  Technometrics 55 (2013) 37–46.
  * Kennedy and O’Hagan (2000) M. C. Kennedy, A. O’Hagan,  Predicting the output from a complex computer code when fast approximations are available,  Biometrika 87 (2000) 1–13.
  * Pilania et al. (2017) G. Pilania, J. E. Gubernatis, T. Lookman,  Multi-fidelity machine learning models for accurate bandgap predictions of solids,  Computational Materials Science 129 (2017) 156–163.
  * Cutajar et al. (2019) K. Cutajar, M. Pullin, A. Damianou, N. Lawrence, J. González,  Deep gaussian processes for multi-fidelity modeling,  arXiv preprint arXiv:1903.07320 (2019).
  * Tuo et al. (2014) R. Tuo, C. J. Wu, D. Yu,  Surrogate modeling of computer experiments with different mesh densities,  Technometrics 56 (2014) 372–380.
  * Fang et al. (2005) K. Fang, R. Li, A. Sudjianto, Design and Modeling for Computer Experiments, volume 6, Chapman & Hall/CRC, 2005.
  * Qian (2009) P. Z. Qian,  Nested latin hypercube designs,  Biometrika 96 (2009) 957–970.
  * He and Qian (2011) X. He, P. Z. Qian,  Nested orthogonal array-based latin hypercube designs,  Biometrika 98 (2011) 721–731.
  * Rennen et al. (2010) G. Rennen, B. Husslage, E. R. Van Dam, D. Den Hertog,  Nested maximin latin hypercube designs,  Structural and Multidisciplinary Optimization 41 (2010) 371–395.
  * Qian (2012) P. Z. Qian,  Sliced latin hypercube designs,  Journal of the American Statistical Association 107 (2012) 393–399.
  * Hwang et al. (2016) Y. Hwang, X. He, P. Z. Qian,  Sliced orthogonal array-based latin hypercube designs,  Technometrics 58 (2016) 50–61.
  * Qian et al. (2006) Z. Qian, C. C. Seepersad, V. R. Joseph, J. K. Allen, C. Jeff Wu,  Building surrogate models based on detailed and approximate simulations (2006).
  * Ehara and Guillas (2021) A. Ehara, S. Guillas,  An adaptive strategy for sequential designs of multilevel computer experiments,  arXiv preprint arXiv:2104.02037 (2021).
  * Sung et al. (2022) C.-L. Sung, Y. Ji, T. Tang, S. Mak,  Stacking designs: designing multi-fidelity computer experiments with confidence,  arXiv preprint arXiv:2211.00268 (2022).
  * He and Qian (2016) X. He, P. Z. Qian,  A central limit theorem for nested or sliced latin hypercube designs,  Statistica Sinica (2016) 1117–1128.
  * Giles (2008) M. B. Giles,  Multilevel monte carlo path simulation,  Operations research 56 (2008) 607–617.
  * Giles (2015) M. B. Giles,  Multilevel monte carlo methods,  Acta numerica 24 (2015) 259–328.
  * Heinrich (2001) S. Heinrich,  Multilevel monte carlo methods,  in: International Conference on Large-Scale Scientific Computing, Springer, 2001, pp. 58–67.
  * Ng and Willcox (2014) L. W. Ng, K. E. Willcox,  Multifidelity approaches for optimization under uncertainty,  International Journal for numerical methods in Engineering 100 (2014) 746–772.
  * Haji-Ali et al. (2016) A.-L. Haji-Ali, F. Nobile, R. Tempone,  Multi-index monte carlo: when sparsity meets sampling,  Numerische Mathematik 132 (2016) 767–806.
  * Peherstorfer et al. (2016a) B. Peherstorfer, T. Cui, Y. Marzouk, K. Willcox,  Multifidelity importance sampling,  Computer Methods in Applied Mechanics and Engineering 300 (2016a) 490–509.
  * Peherstorfer et al. (2016b) B. Peherstorfer, K. Willcox, M. Gunzburger,  Optimal model management for multifidelity monte carlo estimation,  SIAM Journal on Scientific Computing 38 (2016b) A3163–A3194.
  * Geraci et al. (2017) G. Geraci, M. S. Eldred, G. Iaccarino,  A multifidelity multilevel monte carlo method for uncertainty propagation in aerospace applications,  in: 19th AIAA non-deterministic approaches conference, 2017, p. 1951.
  * Teckentrup et al. (2015) A. L. Teckentrup, P. Jantsch, C. G. Webster, M. Gunzburger,  A multilevel stochastic collocation method for partial differential equations with random input data,  SIAM/ASA Journal on Uncertainty Quantification 3 (2015) 1046–1074.
  * Tuo and Wang (2020) R. Tuo, W. Wang,  Kriging prediction with isotropic Matérn correlations: Robustness and experimental designs.,  J. Mach. Learn. Res. 21 (2020) 187–1.
  * Forrester et al. (2007) A. I. Forrester, A. Sóbester, A. J. Keane,  Multi-fidelity optimization via surrogate modelling,  Proceedings of the royal society a: mathematical, physical and engineering sciences 463 (2007) 3251–3269.
  * Kennedy and O’Hagan (2001) M. Kennedy, A. O’Hagan,  Bayesian calibration of computer models,  Journal of the Royal Statistical Society: Series B 63 (2001) 425–464.
  * Wilson and Adams (2013) A. Wilson, R. Adams,  Gaussian process kernels for pattern discovery and extrapolation,  in: International conference on machine learning, PMLR, 2013, pp. 1067–1075.
  * Schulz et al. (2018) E. Schulz, M. Speekenbrink, A. Krause,  A tutorial on gaussian process regression: Modelling, exploring, and exploiting functions,  Journal of Mathematical Psychology 85 (2018) 1–16.
  * Halton (1964) J. H. Halton,  Algorithm 247: Radical-inverse quasi-random point sequence,  Communications of the ACM 7 (1964) 701–702.
  * Chen and Xiong (2017) D. Chen, S. Xiong,  Flexible nested latin hypercube designs for computer experiments,  Journal of Quality Technology 49 (2017) 337–353.

##  Appendix A Proofs for Theorem 1-3

In this section, we prove Theorem 1-3.

###  A.1 Proof of Theorem 1

###### Proof.

By definition,

$\displaystyle\hat{y}_{K}(x)$
$\displaystyle=\mathbb{E}\left[\sum_{i=0}^{K}\delta_{i}(x)\Big{|}y_{K}(\mathcal{X}_{K}),...,y_{0}(\mathcal{X}_{0})\right]$
$\displaystyle=\sum_{i=0}^{K}\mathbb{E}\left[\delta_{i}(x)\Big{|}y_{K}(\mathcal{X}_{K}),...,y_{0}(\mathcal{X}_{0})\right]$

Because $\delta_{i}$’s are mutually independent, for each $i$ we can get

$\displaystyle\mathbb{E}\left[\delta_{i}(x)\Big{|}y_{K}(\mathcal{X}_{K}),\ldots,y_{0}(\mathcal{X}_{0})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\delta_{i}(x)\Big{|}y_{K}(\mathcal{X}_{K})-y_{K-1}(\mathcal{X}_{k}),\ldots,y_{1}(\mathcal{X}_{1})-y_{0}(\mathcal{X}_{0}),y_{0}(\mathcal{X}_{0})\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\delta_{i}(x)\Big{|}y_{i}(\mathcal{X}_{i})-y_{i-1}(\mathcal{X}_{i})\right].$

Then we can rewrite the $\hat{y}_{K}(x)$ as

$\hat{y}_{K}(x)=\sum_{i=0}^{K}\mathbb{E}\left[\delta_{i}(x)\Big{|}y_{i}(\mathcal{X}_{i})-y_{i-1}(\mathcal{X}_{i})\right].$

∎

###  A.2 Proof of Theorem 2

###### Proof.

By definition,

$\displaystyle\mathbb{E}\left\|y_{K}-y_{\infty}\right\|^{2}_{L_{2}}$
$\displaystyle=\mathbb{E}\left\|\sum_{i=K+1}^{\infty}\delta_{i}\right\|^{2}_{L_{2}}$
$\displaystyle=\sum_{i=K+1}^{\infty}\mathbb{E}\left\|\delta_{i}\right\|^{2}_{L_{2}}$
$\displaystyle=\sum_{i=K+1}^{\infty}\mathbb{E}\left(\int_{\mathfrak{X}}\delta_{i}^{2}dx\right)$
$\displaystyle=\sum_{i=K+1}^{\infty}\int_{\mathfrak{X}}\mathbb{E}(\delta_{i}^{2})dx$
$\displaystyle=\sum_{i=K+1}^{\infty}C_{s}\lambda^{2i}\sigma^{2}$
$\displaystyle=C_{l}^{2}\lambda^{2K}\sigma^{2},$

where $\mathfrak{X}$ denotes the input domain of $x$, $C_{s}$ is the area of
$\mathfrak{X}$ and constant $C_{l}=\sqrt{C_{s}}/(1-\lambda)$. ∎

###  A.3 Proof of Theorem 4

###### Proof.

Based on the Theorem 1, we have

$\displaystyle\mathbb{E}\left\|y_{K}-\hat{y}_{K}\right\|^{2}_{L_{2}}$
$\displaystyle=\mathbb{E}\left\|\sum_{i=0}^{k}(I_{i}(\delta_{i})-\delta_{i})\right\|^{2}_{L_{2}}$
$\displaystyle=\int_{\mathfrak{X}}\sum_{i=0}^{k}(I_{i}(\delta_{i})-\delta_{i})^{2}dx$
$\displaystyle=\int_{\mathfrak{X}}\sum_{i=0}^{k}\mathbb{E}((I_{i}(\delta_{i})-\delta_{i})^{2})dx$
$\displaystyle\leq\sum_{i=0}^{k}p\sigma^{2}\lambda^{2i}n_{i}^{-\nu/d},$

where $n_{i}$ is the number of samples in the $i$-th step, $p$ is a constant
only based on the correlation function, $d$ is the dimension of the input
space and $\nu$ is a constant. ∎

##  Appendix B Solve the optimization problem

Consider the optimization problem

$\displaystyle\min_{n_{i}}\quad$ $\displaystyle\sum_{i=0}^{K}n_{i}C_{i}$ (26)
s.t.
$\displaystyle\sum_{i=0}^{K}\lambda^{2i}\sigma^{2}n_{i}^{-\nu/d}\leq\epsilon^{2}/2$
(27)

We form the Lagrange function of (26) as

$\displaystyle
L(n_{0},n_{1},...,n_{K},\eta)=\sum_{i=0}^{K}n_{i}C_{a}a^{i}+\eta(\sum_{i=0}^{K}p\lambda^{2i}\sigma^{2}n_{i}^{-\nu/d}-\epsilon^{2}/2).$

To find a local optimum, we need to have $\nabla L=0$, leading to $K+2$
conditions

$\displaystyle\frac{\partial L}{\partial
n_{i}}=C_{a}a^{i}-\frac{p\eta\nu}{d}\lambda^{2i}\sigma^{2}n_{i}^{-(\nu/d+1)}=0,\quad
i=0,1,...,K$ (28) $\displaystyle\frac{\partial
L}{\partial\eta}=\sum_{i=0}^{K}p\lambda^{2i}\sigma^{2}n_{i}^{-\nu/d}-\epsilon^{2}/2=0.$
(29)

Solving the $(K+1)$ equations in (28) leads to

$\displaystyle
n_{i}=(\frac{pC_{a}d}{\nu\sigma^{2}}\frac{a^{i}}{\lambda^{2i}}\eta^{-1})^{-\frac{d}{d+\nu}},\quad
i=0,1,...,K.$ (30)

Substitute (30) into (29). Then we have

$\displaystyle\eta=\frac{C_{a}d}{\nu}(\sigma^{2})^{\frac{d}{\nu}}(\frac{p\epsilon^{2}}{2S})^{-\frac{d+\nu}{\nu}},$
(31)

where

$\displaystyle
S=\sum_{i=0}^{K}(a^{i})^{\frac{\nu}{d+\nu}}(\lambda^{2i})^{\frac{d}{d+\nu}}.$

Substitute (31) into (30). We get

$\displaystyle
n_{i}=(\frac{\epsilon^{2}}{2p\sigma^{2}S})^{-\frac{d}{\nu}}(\frac{a^{i}}{\lambda^{2i}})^{-\frac{d}{d+\nu}},\quad
i=0,1,...,K.$

Since $n_{i}$ is an integer, we choose

$\displaystyle n_{i}^{*}=\lceil n_{i}\rceil$

to get an approximate solution to the original problem.

##  Appendix C Analysis of cost for the MLGP

In this section, we make the analysis of the cost based on the $K$ and
$n_{i}^{*}$ given in the above section. The total cost is given by

$\displaystyle
C_{\epsilon}=\sum_{i=0}^{K}n_{i}^{*}C_{i}\leq\sum_{i=0}^{K}n_{i}^{*}C_{a}a^{i}.$

Since $\lceil n_{i}\rceil\leq n_{i}+1$, we have

$\displaystyle C_{\epsilon}$ $\displaystyle\leq$
$\displaystyle\sum_{i=0}^{K}Ca(n_{i}+1)a^{i}$ (32) $\displaystyle\eqsim$
$\displaystyle\sum_{i=0}^{K}(S\sigma^{2}(\epsilon^{2}/2)^{-1})^{\frac{d}{\nu}}\sum_{i=0}^{K}(a^{\frac{\nu}{d+\nu}}\lambda^{\frac{2d}{d+\nu}})^{i}+\sum_{i=0}^{K}a^{i}$
$\displaystyle=$
$\displaystyle(2\sigma^{2}\epsilon^{-2})^{\frac{d}{v}}S^{\frac{d+\nu}{\nu}}+\sum_{i=0}^{K}a^{i}.$

We denote the first part on the right-hand side of (32) as part (I) and the
second part as (II), and consider them separately.

If $(a^{i})^{\frac{\nu}{d+\nu}}(\lambda^{2i})^{\frac{d}{d+\nu}}=1$ (i.e.,
$\alpha\nu=d\beta$), then

$\displaystyle(I)$ $\displaystyle=$
$\displaystyle(2\sigma^{2}\epsilon^{-2})^{\frac{d}{v}}(K+1)^{\frac{d+\nu}{\nu}}$
$\displaystyle\eqsim$
$\displaystyle\epsilon^{-2d/\nu}(2\log_{\lambda^{2}}\frac{\epsilon}{2C_{l}\sigma})^{\frac{d+\nu}{\nu}}$
$\displaystyle\eqsim$
$\displaystyle\epsilon^{-2d/\nu}|\ln\epsilon^{2}|^{\frac{d+\nu}{\nu}};$
$\displaystyle(II)$ $\displaystyle=$
$\displaystyle\frac{a^{K}(1-a^{-(K+1)})}{1-a^{-1}}$ $\displaystyle\eqsim$
$\displaystyle a^{K}\eqsim(\epsilon^{2})(\ln a/\ln\lambda^{2})$
$\displaystyle=$ $\displaystyle\epsilon^{-2d/\nu}.$

Since
$\epsilon^{-2d/\nu}|\ln\epsilon^{2}|^{\frac{d+\nu}{\nu}}>\epsilon^{-2d/\nu}$,
$C_{\epsilon}\lesssim\epsilon^{-2d/\nu}|\ln\epsilon^{2}|^{\frac{d+\nu}{\nu}}$.

If $(a^{i})^{\frac{\nu}{d+\nu}}(\lambda^{2i})^{\frac{d}{d+\nu}}<1$ (i.e.,
$\alpha\nu<d\beta$), $S$ will converge to a limitation which is independent of
$K$. As a result, we have

$\displaystyle(I)\lesssim\epsilon^{-2d/\nu};$
$\displaystyle(II)\lesssim\epsilon^{-2\alpha/\beta}.$

Since $d/\nu>\alpha/\beta$, $C_{\epsilon}\lesssim\epsilon^{-2d/\nu}$.

If $(a^{i})^{\frac{\nu}{d+\nu}}(\lambda^{2i})^{\frac{d}{d+\nu}}>1$ (i.e.,
$\alpha\nu>d\beta$), we have

$\displaystyle(I)$ $\displaystyle\lesssim$
$\displaystyle\epsilon^{-2d/\nu}(\exp{\frac{\alpha\nu-d\beta}{d+\nu}})^{\frac{K(d+\nu)}{\nu}}$
$\displaystyle\lesssim$
$\displaystyle\epsilon^{-2d/\nu}\exp{(-\frac{(\ln\epsilon^{2}(\alpha\nu-d\beta))}{\nu\beta})}$
$\displaystyle=$ $\displaystyle\epsilon^{-2\alpha/\beta};$ $\displaystyle(II)$
$\displaystyle\eqsim$ $\displaystyle\epsilon^{-2\alpha/\beta}.$

As a result, $C_{\epsilon}\lesssim\epsilon^{-2\alpha/\beta}$.

In summary,

$\displaystyle
C_{\epsilon}\lesssim\left\\{\begin{matrix}\epsilon^{-2d/\nu},&\alpha\nu<d\beta\\\
\epsilon^{-2d/\nu}|\ln\epsilon^{2}|^{\frac{d+\nu}{\nu}},&\alpha\nu=d\beta\\\
\epsilon^{-2\alpha/\beta},&\alpha\nu>d\beta\end{matrix}\right.$

##  Appendix D Analysis of cost for single-level

The error caused by single-level approximation can be bounded as

$\displaystyle\mathbb{E}(\left\|y_{\infty}-\hat{y}_{K}^{SL}\right\|^{2})=\mathbb{E}(\left\|y_{\infty}-y_{K}^{SL}\right\|^{2}+\left\|y_{K}^{SL}-\hat{y}_{K}^{SL}\right\|^{2}),$

where $\hat{y}_{K}^{SL}=I(\delta_{K})$. Then the objective function for this
single-level approximation can be summarized as

$\displaystyle{2}\min_{K,n_{K}}\quad$ $\displaystyle n_{K}C_{K}$ (33) s.t.
$\displaystyle\mathbb{E}\left\|y_{\infty}-\hat{y}_{K}^{SL}\right\|^{2}_{L^{2}}\leq\epsilon$
(34)

Based on Theorems 1 and 2, we have

$\displaystyle\mathbb{E}\left\|y_{\infty}-\hat{y}_{K}^{SL}\right\|^{2}\gtrsim\lambda^{2K}\sigma^{2}+\sigma^{2}n_{K}^{-\nu/d},$

where $n_{K}$ is the number of samples. It should be sufficient to bound the
upper bound of the error by $\epsilon$, that is,

$\displaystyle\lambda^{2K}+n_{K}^{-\nu/d}\lesssim\epsilon^{2}.$ (35)

From the (35), we can get the lower bound for $n_{K}$, that is,

$\displaystyle n_{K}\gtrsim(\epsilon^{2}-\lambda^{2K})^{-d/\nu}.$

Then we have the lower bound for the cost of single-level approximation:

$\displaystyle C_{\epsilon}^{SL}=n_{K}C_{k}\gtrsim
a^{K}(\epsilon^{2}-\lambda^{2K})^{-d/\nu}.$ (36)

We denote $a=\exp{(\alpha)}$ and $\lambda^{2}=\exp{(-\beta)}$ for convenience,
$\alpha,\beta>0$. Then wen can rewrite (36) as

$\displaystyle\ln C_{\epsilon}^{SL}\gtrsim\alpha
K-\frac{d}{\nu}\ln{(\epsilon^{2}-e^{-\beta K})}.$ (37)

Now we want to find the lower bound of $C_{\epsilon}^{SK}$, that is, we need
to find the lower bound for the right-hand side of (37). Consider the
following function with respect to $K$:

$\displaystyle T(K)=\alpha K-\frac{d}{\nu}\ln{(\epsilon^{2}-e^{-\beta K})}.$

Set the derivative of $T$ with respect to $K$ equal to $0$:

$\displaystyle\frac{dT}{dK}$ $\displaystyle=$
$\displaystyle\alpha+\frac{d\beta}{\nu}(1-\frac{\epsilon^{2}}{\epsilon^{2}-e^{-\beta
K}})=0.$

We can get the solution
$K_{0}=-\frac{1}{\beta}\ln\frac{\epsilon^{2}\alpha\nu}{\alpha\nu+\beta d}$.
When $K<K_{0}$, $T(K)$ is decreasing and when $K>K_{0}$, $T(K)$ is increasing.
Substitute $K=0$ into $T(K)$, we can get the lower bound of the cost for the
single-level:

$\displaystyle
C_{\epsilon}^{SL}\gtrsim\epsilon^{-2(\frac{d}{\nu}+\frac{\alpha}{\beta})}.$