[table]capposition=top

# Multi-Fidelity Cost-Aware Bayesian Optimization

Zahra Zanjani Foumani  Department of Mechanical and Aerospace Engineering,
University of California, Irvine  Mehdi Shishehbor  Department of Mechanical
and Aerospace Engineering, University of California, Irvine  Amin Yousefpour
Department of Mechanical and Aerospace Engineering, University of California,
Irvine  Ramin Bostanabad Corresponding Author<EMAIL_ADDRESS>
Gitlab repository: https://gitlab.com/S3anaz/multi-fidelity-cost-aware-
bayesian-optimization/-/tree/main/ Department of Mechanical and Aerospace
Engineering, University of California, Irvine

Abstract
Bayesian optimization (BO) is increasingly employed in critical applications
such as materials design and drug discovery. An increasingly popular strategy
in BO is to forgo the sole reliance on high-fidelity data and instead use an
ensemble of information sources which provide inexpensive low-fidelity data.
The overall premise of this strategy is to reduce the overall sampling costs
by querying inexpensive low-fidelity sources whose data are correlated with
high-fidelity samples. Here, we propose a multi-fidelity cost-aware BO
framework that dramatically outperforms the state-of-the-art technologies in
terms of efficiency, consistency, and robustness. We demonstrate the
advantages of our framework on analytic and engineering problems and argue
that these benefits stem from our two main contributions: $(1)$ we develop a
novel acquisition function for multi-fidelity cost-aware BO that safeguards
the convergence against the biases of low-fidelity data, and $(2)$ we tailor a
newly developed emulator for multi-fidelity BO which enables us to not only
simultaneously learn from an ensemble of multi-fidelity datasets, but also
identify the severely biased low-fidelity sources that should be excluded from
BO.

Keywords: Bayesian optimization; multi-fidelity modeling; emulation; resource
allocation; and Gaussian process.

##  1 Introduction

Bayesian optimization (BO) is an iterative and sample-efficient global
optimization technique that has been successfully applied to a wide range of
applications including materials discovery [1, 2, 3, 4], design of chemical
systems such as catalysts [5], hyperparameter tuning in machine learning (ML)
models [6], robot motion control [7], and updating internet-scale software
systems [8]. While BO is very effective, the total optimization cost can still
be high if only an expensive source (e.g., experiments or costly simulations)
is sampled during the optimization. To reduce the overall data collection
costs in such scenarios an increasingly popular strategy is to formulate
multi-fidelity (MF) methods that use multiple data sources which typically
have different levels of accuracy and cost. Assuming low-fidelity (LF) sources
are cheaper to query, the overall premise of these methods is to reduce the
total sampling costs by leveraging the correlations between low- and high-
fidelity (HF) data. In this paper, we propose a multi-fidelity cost-aware
(MFCA) BO framework that optimizes an expensive objective function using an
ensemble of data sources with arbitrary levels of accuracy and cost. We
provide a new perspective on probabilistic learning from multiple sources
which endows our framework with five major advantages over existing MF BO
techniques: $(a)$ safeguarding the convergence against the biases of the LF
sources even if they are extremely inexpensive to query (i.e., if the majority
of the samples are LF), $(b)$ learning the relative fidelity of the sources
rather than requiring a priori determination of such relations by the user,
$(c)$ dispensing with the assumptions that aim to relate the fidelity and cost
of a data source, and $(d)$ improving numerical stability and efficiency.

As schematically demonstrated in Figure 1 and detailed in Section 2.2, BO has
two main ingredients that interact sequentially to optimize a black-box and
expensive-to-evaluate objective function. These two ingredients include an
emulator (i.e., a probabilistic surrogate) and an acquisition function (AF).
The optimization process starts by fitting the emulator to a small initial
dataset that is typically obtained via design of experiments. The emulator is
next used in the AF to determine the candidate input(s) whose corresponding
output(s) must be obtained by querying the expensive function. Given the new
sample(s), the training dataset is updated and the entire fitting-searching-
sampling process is repeated until a convergence criterion is met (e.g.,
resources are exhausted).

Figure 1: Bi-fidelity BO: The left box of the figure shows two different
scenarios (without LF data and with LF Data) for a function with two numerical
inputs ($x_{1}$ , $x_{2}$) and one categorical input ($t$) that has two levels
($A$ and $Z$). In each scenario, $x_{2}$ is fixed to simplify the
visualizations and the only difference between the upper and lower plots is
the level of $t$. In the left panel only HF data is used in BO while in the
right panel two data sources are available where one of them has a lower
fidelity and sampling cost. In the multi-fidelity case, the correlations
between the data sources are used to decrease the overall sampling costs by
reducing the reliance on the HF source.

While many emulators such as Bayesian neural networks are available, Gaussian
processes (GPs) are typically used in BO since they very efficiently learn
from small data, are easy and fast to train, provide prediction uncertainties,
and have interpretable parameters [9, 10, 11, 12]. Use of GPs in BO has
increased even more because of the recent works that enable them to handle
categorical variables [13, 14, 15], high-dimensional inputs [16, 17], large
datasets [18, 19, 20, 21], and non-stationary noise [22, 23]. As for the AF,
there are many choices [24, 25] such as expected improvement (EI), probability
of improvement (PI), and knowledge gradient (KG). The primary difference among
these AFs is that they select the candidate input(s) while taking different
approaches for balancing exploitation (i.e., sampling based on the best
predictions of the emulator) and exploration (i.e., sampling to reduce
prediction uncertainty). This selection involves integration which can
sometimes be analytically computed (e.g., when GP and EI are chosen as the
emulator and AF, respectively).

To optimally use an ensemble of information sources in BO, two conditions must
be met: the emulator should leverage the cross-source correlations (which are
hidden in the datasets) to more accurately surrogate all the data sources
(esp. the HF one) and the AF should appropriately calculate the value or
utility of a to-be-sampled data point based on its source and evaluation cost.
Satisfying these two conditions in many realistic applications is nontrivial
for the following reasons:

  * •

The global optima (input and output) of LF sources differ from those of the HF
source, see Figure 4 for a one dimensional illustration.

  * •

Some LF sources have major biases (which are a priori unknown to the analyst)
and must be excluded from the search process from the very beginning.
Including such LF sources increases the overall sampling costs and may also
result in convergence to an incorrect solution.

  * •

If highly cheap LF sources are available, a naively designed AF chooses to
sample from them very frequently since the information value of a candidate
point is inversely scaled by the cost of its source. This heavily biased
sampling can force BO to converge to the optima of those sources rather than
the HF source.

As reviewed in Section 2.2, existing multi-fidelity BO technologies partially
address these challenges by ad hoc tuning and making simplifying assumptions.
These assumptions often include presuming a direct relationship between
fidelity level and sampling cost, assuming LF sources are always useful, or
manually adjusting the sampling costs (e.g., converting the $1000/50/1$ cost
ratio between three sources to $1/0.5/0.1$). These manual changes are quite
laborious and result in either convergence to an incorrect solution or higher
overall costs compared to single-fidelity (SF) BO which solely leverages the
HF source.

In this paper, we provide new perspectives for learning from multi-fidelity
sources in the context of BO. In particular, we argue that $(a)$ the emulator
must fuse the multi-source data in a nonlinearly learnt manifold to maximally
leverage cross-source correlations and also indicate trustworthy LF sources
that do not deteriorate BO’s performance, and $(b)$ the AF should use the
available information on the LF sources (HF source) solely for exploration
(exploitation). As demonstrated in Section 4, these contributions endow our
framework with significant performance improvements over existing
technologies.

The rest of the paper is organized as follows. We review the relevant
technical background in Section 2 and introduce our approach in Section 3. We
test the performance of our MFCA BO framework against the state-of-the-art on
a set of analytic and two engineering problems in Section 4 and finally
conclude the paper in section 5. We also provide a nomenclature and list of
symbols in Appendix E. The Gitlab repository associated with this project
hosts supplementary materials.

##  2 Technical Background

In this Section, we first provide some background on latent map Gaussian
process (LMGP) which is one of the key components of our MFCA BO framework.
Then, we elaborate on the two main ingredients of BO in Section 2.2 where we
also review some of the existing methods for handling MF data.

###  2.1 Latent Map Gaussian Processes (LMGPs)

For metamodeling via GPs, one assumes that the training data comes from a
multivariate normal distribution with parametric mean and covariance functions
and then uses closed-form conditional distribution formulas for prediction.
Below, we first detail the process for estimating these parameters when the
input space contains categorical and/or numerical variables. Then, we provide
the prediction formulas.

Assume the training data is a realization from a GP and that the following
relation holds:

$\begin{split}y(\boldsymbol{x})=\boldsymbol{m}(\boldsymbol{x})\boldsymbol{\beta}+\xi(\boldsymbol{x})\end{split}$
(1)

where $\boldsymbol{x}=[x_{1},x_{2},\ldots,x_{dx}]^{T}$ is the input vector,
$y(\boldsymbol{x})$ is the output/response,
$\boldsymbol{m}(\boldsymbol{x})=\left[m_{1}(\boldsymbol{x}),\ldots,m_{d\beta}(\boldsymbol{x})\right]$
are a set of parametric basis functions,
$\boldsymbol{\beta}=[\beta_{1},\beta_{2},...,\beta_{d\beta}]^{T}$ are unknown
coefficients, and $\xi(\boldsymbol{x})$ is a zero-mean GP whose covariance
function or kernel is:

$\begin{split}\operatorname{cov}\left(\xi(\boldsymbol{x}),\xi\left(\boldsymbol{x}^{\prime}\right)\right)=c\left(\boldsymbol{x},\boldsymbol{x}^{\prime}\right)=\sigma^{2}r\left(\boldsymbol{x},\boldsymbol{x}^{\prime}\right)\end{split}$
(2)

where $\sigma^{2}$ is the variance of the process and $r(.,.)$ is a parametric
correlation function. A common choice for $r(.,.)$ is:

$\begin{split}r(\boldsymbol{x},\boldsymbol{x}^{\prime})=\exp\\{-\sum_{i=1}^{dx}10^{\omega_{i}}(x_{i}-x_{i}^{\prime})^{2}\\}=\exp\\{-(\boldsymbol{x}-\boldsymbol{x}^{\prime})^{T}10^{\boldsymbol{\Omega}}(\boldsymbol{x}-\boldsymbol{x}^{\prime})\\}\end{split}$
(3)

where
$\boldsymbol{\omega}=[\omega_{1},\ldots,\omega_{dx}]^{T},-\infty<\omega_{i}<\infty$
are the scale parameters and
$\boldsymbol{\Omega}=\operatorname{diag}(\boldsymbol{\omega})$. GP modeling
largely depends on the choice of the correlation function which measures the
weighted distance between any two inputs, see Eq. 3. As recently motivated in
[26], to directly use GPs for MF modeling we must extend them such that they
can handle categorical inputs. This extension primarily relies on
reformulating $r(.,.)$ and can be accomplished in different ways. In this
paper, we employ LMGPs [26] since $(1)$ they have been shown to outperform
other approaches, and $(2)$ they provide a visualizable and interpretable
manifold which can be used to detect discrepancies among data sources (this
manifold helps us to exclude highly biased sources from BO).

Let us denote categorical inputs via
$\boldsymbol{t}=[t_{1},\ldots,t_{dt}]^{T}$ where variable $t_{i}$ has $l_{i}$
distinct levels. For example, $t_{1}=\\{Male,Female\\}$ and
$t_{2}=\\{Persian,American,Spanish\\}$ are two categorical inputs where
$l_{1}=2$ and $l_{2}=3$. To handle mixed inputs,
$\boldsymbol{u}=[x_{1},...,x_{dx},t_{1},...,t_{dt}]^{T}$, LMGP learns a
unified parametric function that maps each combination of categorical
variables to a point in a manifold or latent space111Multiple mapping
functions can also be used and we leverage this in Section 3.1 where two
mapping functions are learned to build two manifolds.. This mapping function
can be incorporated into any standard correlation function such as the
Gaussian which is reformulated as follows for mixed inputs:

$\begin{split}r(\boldsymbol{u},\boldsymbol{u}^{\prime})=\exp\\{-(\boldsymbol{x}-\boldsymbol{x}^{\prime})^{T}\boldsymbol{\Omega}(\boldsymbol{x}-\boldsymbol{x}^{\prime})-\|\boldsymbol{z}(\boldsymbol{t})-\boldsymbol{z}(\boldsymbol{t}^{\prime})\|_{2}^{2}\\}\end{split}$
(4)

or equivalently,

$\begin{split}r(\boldsymbol{u},\boldsymbol{u}^{\prime})=\exp\\{-\sum_{i=1}^{dx}10^{\omega_{i}}(x_{i}-x_{i}^{\prime})^{2}\\}\times\exp\\{-\sum_{i=1}^{dz}(z_{i}(\boldsymbol{t})-z_{i}(\boldsymbol{t}^{\prime}))^{2}\\}\end{split}$
(5)

where $\|\cdot\|_{2}$ denotes the Euclidean 2-norm and
$\boldsymbol{z}(\boldsymbol{t})=\left[z_{1}(\boldsymbol{t}),\ldots,z_{dz}(\boldsymbol{t})\right]_{1\times
dz}$ is the location in the learned manifold corresponding to the specific
combination of the categorical variables denoted by $\boldsymbol{t}$. To find
these latent points, LMGP first assigns a unique prior representation (a
unique vector) to each combination of categorical variables. Then, it learns a
linear transformation222More complex transformations based on, e.g., deep
neural networks can also be used. that maps these unique vectors into a
compact manifold (aka latent space) with dimensionality $dz$:

$\begin{split}\boldsymbol{z(t)=\zeta(t)A}\end{split}$ (6)

where $\boldsymbol{t}$ denotes a specific combination of the categorical
variables, $\boldsymbol{z(t)}$ is the $1\times dz$ posterior latent
representation of $\boldsymbol{t}$, $\boldsymbol{\zeta}(\boldsymbol{t})$ is
the unique prior vector representation of $\boldsymbol{t}$, and
$\boldsymbol{A}$ is a rectangular matrix that maps $\boldsymbol{\zeta(t)}$ to
$\boldsymbol{z(t)}$. There are various methods for constructing the prior
vectors $\boldsymbol{\zeta}$ and we refer the reader to [26] for more details.
In this paper, we use grouped one-hot encoding which makes
$\boldsymbol{\zeta}(\boldsymbol{t})$ and $\boldsymbol{A}$ to be of sizes
$1\times\sum_{i=1}^{dt}l_{i}$ and $\sum_{i=1}^{dt}l_{i}\times dz$,
respectively. For instance, in the above example the grouped one-hot encoded
version of the combination $\boldsymbol{t}=[Female,American]^{T}$ is
$\boldsymbol{\zeta(\boldsymbol{t})}=[0,1,0,1,0]$ where the first two numbers
encode the levels of $t_{1}$ while the last three numbers indicate the levels
of $t_{2}$.

To emulate via LMGP, the hyperparameters ($\boldsymbol{\beta,A,\omega}\text{,
and }\sigma^{2}$) must be estimated via the training data. To find these
estimates, we utilize maximum likelihood estimation (MLE) due to its
computational efficiency over fully Bayesian techniques. The MLE estimates the
hyperparameters such that they maximize the likelihood of the $n$ training
data being generated by $y(\boldsymbol{x})$, that is:

$\begin{split}[\widehat{\boldsymbol{\beta}},\hat{\sigma},\widehat{\boldsymbol{\omega}},\widehat{\boldsymbol{A}}]=\underset{\boldsymbol{\beta},\sigma^{2},\omega,\boldsymbol{A}}{\operatorname{argmax}}\left|2\pi\sigma^{2}\boldsymbol{R}\right|^{-\frac{1}{2}}\times\exp\left\\{\frac{-1}{2}(\boldsymbol{y}-\boldsymbol{M}\boldsymbol{\beta})^{T}\left(\sigma^{2}\boldsymbol{R}\right)^{-1}(\boldsymbol{y}-\boldsymbol{M}\boldsymbol{\beta})\right\\}\end{split}$
(7)

or equivalently:

$\begin{split}[\widehat{\boldsymbol{\beta}},\hat{\sigma},\widehat{\boldsymbol{\omega}},\widehat{\boldsymbol{A}}]=\underset{\boldsymbol{\beta},\sigma^{2},\omega,\boldsymbol{A}}{\operatorname{argmin}}\frac{n}{2}\log\left(\sigma^{2}\right)+\frac{1}{2}\log(|\boldsymbol{R}|)+\frac{1}{2\sigma^{2}}(\boldsymbol{y}-\boldsymbol{M}\boldsymbol{\beta})^{T}\boldsymbol{R}^{-1}(\boldsymbol{y}-\boldsymbol{M}\boldsymbol{\beta}),\end{split}$
(8)

where $\log(\cdot)$ is the natural logarithm, $|\cdot|$ denotes the
determinant operator, $\boldsymbol{y}=[y^{1},\ldots,y^{n}]^{T}$ is the
$n\times 1$ vector of outputs in the training data, $\boldsymbol{R}$ is the
$n\times n$ correlation matrix with the $(i,j)^{th}\text{ element
}R_{ij}=r(\boldsymbol{x}^{i},\boldsymbol{x}^{j})\text{ for }i,j=1,\ldots,n$ ,
and $\boldsymbol{M}$ is the $n\times d\beta$ matrix with the $(i,j)^{th}$
element $M_{ij}=m_{j}(\boldsymbol{x}^{i})\text{ for }i=1,\ldots,n\text{ and
}j=1,\ldots,d\beta$.

Using the method of profiling [27], Eq. 8 can be simplified to:

$\begin{split}[\widehat{\boldsymbol{\omega}},\widehat{\boldsymbol{A}}]=\underset{\boldsymbol{\omega},\boldsymbol{A}}{\operatorname{argmin}}\hskip
8.53581ptn\log\left(\hat{\sigma}^{2}\right)+\log(|\boldsymbol{R}|)=\underset{\boldsymbol{\omega},\boldsymbol{A}}{\operatorname{argmin}}\hskip
8.53581ptL,\end{split}$ (9)

where
$\hat{\sigma}^{2}=\frac{1}{n}(\boldsymbol{y}-\boldsymbol{M}\hat{\boldsymbol{\beta}})^{T}\boldsymbol{R}^{-1}(\boldsymbol{y}-\boldsymbol{M}\hat{\boldsymbol{\beta}})$
and
$\hat{\boldsymbol{\beta}}=[\boldsymbol{M}^{T}\boldsymbol{R}^{-1}\boldsymbol{M}]^{-1}[\boldsymbol{M}^{T}\boldsymbol{R}^{-1}\boldsymbol{y}]$.
Eq. 9 can be efficiently solved via a gradient-based optimization technique
[28, 12].

Once the hyperparameters are estimated, the conditional distribution formulas
are used to predict the response distribution at the arbitrary point
$\boldsymbol{x}^{*}$. The mean and variance of this normal distribution are:

$\begin{split}\mathbb{E}[y(\boldsymbol{x}^{*})]=\mu(\boldsymbol{x}^{*})=\boldsymbol{m}(\boldsymbol{x}^{*})\widehat{\boldsymbol{\beta}}+\boldsymbol{r}^{T}(\boldsymbol{x}^{*})\boldsymbol{R}^{-1}(\boldsymbol{y}-\boldsymbol{M}\widehat{\boldsymbol{\beta}})\end{split}$
(10)

$\begin{split}\operatorname{cov}(y(\boldsymbol{x}^{*}),y(\boldsymbol{x}^{*}))=\sigma^{2}(\boldsymbol{x}^{*})=\hat{\sigma}^{2}(1-\boldsymbol{r}^{T}(\boldsymbol{x}^{*})\boldsymbol{R}^{-1}\boldsymbol{r}(\boldsymbol{x}^{*})+\boldsymbol{g}(\boldsymbol{x}^{*})(\boldsymbol{M}^{T}\boldsymbol{R}^{-1}\boldsymbol{M})^{-1}\boldsymbol{g}(\boldsymbol{x}^{*})^{T})\end{split}$
(11)

where $\mathbb{E}$ denotes expectation,
$m(\boldsymbol{x}^{*})=[m_{1}(\boldsymbol{x}^{*}),\ldots,m_{d\beta}(\boldsymbol{x}^{*})]$,
$\boldsymbol{r}(\boldsymbol{x}^{*})$ is an $n\times 1$ vector with the
$i^{th}$ element $r(\boldsymbol{x}^{i},\boldsymbol{x}^{*})$, and
$\boldsymbol{g}(\boldsymbol{x}^{*})=\boldsymbol{m}(\boldsymbol{x}^{*})-\boldsymbol{M}^{T}\boldsymbol{R}^{-1}\boldsymbol{r}(\boldsymbol{x}^{*})$.

###  2.2 Bayesian Optimization (BO)

BO is increasingly used to optimize expensive-to-evaluate (and typically
black-box) functions. As opposed to gradient-based or heuristic optimization
techniques that only rely on function evaluations (or predictions of a
surrogate), BO leverages the probabilistic predictions of a sequentially
updated emulator. In this paper, we focus on single-response functions over
unconstrained and bounded domains but note that BO can also handle multi-
response (or multi-task) [29], composite [30], or constrained problems [31].

As summarized in Algorithm 1, BO has an iterative nature where an emulator is
first fitted to some initial training data. This emulator is then queried via
the AF which calculates the expected utility of any point in the input space,
i.e., $\mathbb{E}[I(\boldsymbol{x})]$ where $I(\boldsymbol{x})$ is the user-
defined utility function. These queries are used in the auxiliary optimization
problem333Gradient-based optimization techniques are almost always used at
this stage. which aims to find the candidate point with the maximum expected
utility in the input space. Once this point is found, the expensive function
is queried and the resulting (input, output) pair is used to update the
training data. Given the augmented dataset, the above process is repeated
until a convergence metric is met.

Except for some special cases, solving the auxiliary optimization problem in
Algorithm 1 is highly costly since each evaluation of its objective function
(i.e., $\mathbb{E}[I(\boldsymbol{x})]$) amounts to integration which cannot be
calculated analytically. Fortunately, the special cases perform quite well in
most practical applications and hence are used frequently (we also employ them
in our framework).

Algorithm 1 is strictly sequential in that the dataset is augmented with a
single sample at each iteration. To leverage parallel computing, one can
augment the dataset with a pool of samples which jointly maximize the expected
utility [32]. Additionally, Algorithm 1 is myopic in that the AF does not
consider the effect of the to-be-evaluated sample on the emulator in the
future iterations. This myopic nature is addressed in look-ahead AFs such as
KG (detailed below) which typically provide improved performance at the
expense of significantly increasing the cost of solving the auxiliary
optimization problem.

Given: Initial data
$\mathcal{D}^{k}=\\{(\boldsymbol{x}^{i},y^{i})\\}_{i=1}^{k}$, expensive black-
box function $f(\boldsymbol{x})$

Define: Utility function $I(\boldsymbol{x})$, stop conditions

while  _stop conditions not met_ do

  1. 1.

Train the GP emulator using $\mathcal{D}^{k}$

  2. 2.

Define the acquisition function
$\alpha(\boldsymbol{x})=\mathbb{E}[I(\boldsymbol{x})\mid\mathcal{D}^{k}]$

  3. 3.

Solve the auxiliary optimization problem:
$\boldsymbol{x}^{k+1}=\underset{{\boldsymbol{x}\in\mathbb{X}}}{\operatorname*{argmax}}\hskip
5.69054pt\alpha(\boldsymbol{x})$

  4. 4.

Query $f(\boldsymbol{x})$ at $\boldsymbol{x}^{k+1}$ to obtain $y^{k+1}$

  5. 5.

Update data:
$\mathcal{D}^{k+1}\leftarrow\mathcal{D}^{k}\cup(\boldsymbol{x}^{k+1},y^{k+1})$

  6. 6.

Update counter: $k\leftarrow k+1$

end while

Output: Updated data
$\mathcal{D}^{k}=\left\\{\left(\boldsymbol{x}^{i},y^{i}\right)\right\\}_{i=1}^{k}$,
GP emulator

Algorithm 1 Strictly Sequential and Myopic Bayesian Optimization for
Maximization

####  2.2.1 Single-Fidelity Acquisition Functions

Many different AFs have been proposed for diverse applications and in this
Section we review three of the most widely used ones: EI, PI, and KG. While
the first two AFs are myopic444EI does have a look-ahead version [33]., KG is
look-ahead. Any AF calculates the expected value of a user-defined utility
function conditioned on the available data $D$, that is:

$\begin{split}\alpha(\boldsymbol{x})=\mathbb{E}[I(\boldsymbol{x})\mid\mathcal{D}]\end{split}$
(12)

Our proposed AF for MFCA BO (see Section 3.2) leverages EI as well as PI and
hence is also myopic. In Section 4 we compare our AF to EI, PI, and KG and
indicate that it consistently outperforms them.

PI is an AF that favors exploitation [34], i.e., it rewards samples that
improve $y^{*}$ which is the best function value seen so far. For instance,
when maximizing the expensive black-box function $f(\boldsymbol{x})$, this AF
uses the following utility function:

$\begin{split}I_{PI}(\boldsymbol{x})=\begin{cases}1&y(\boldsymbol{x})>y^{*}\\\
0&y(\boldsymbol{x})\leq y^{*}\end{cases}\end{split}$ (13)

where $y(\boldsymbol{x})$ is the emulator-based prediction at
$\boldsymbol{x}$. Based on Eq. 13, if $y(\boldsymbol{x})$ is less than
$y^{*}$, the point $\boldsymbol{x}$ has zero utility. Assuming a GP is used
for emulation, $y(\boldsymbol{x})$ follows a normal distribution whose mean
and variance are given in Eq. 10 and Eq. 11, respectively. Using the
reparameterization trick (see Section A) we can show that the resulting AF
based on $I_{PI}(\boldsymbol{x})$ is [35, 36]:

$\begin{split}\alpha_{PI}(\boldsymbol{x})=\Phi\left(\frac{\mu(\boldsymbol{x})-y^{*}}{\sigma(\boldsymbol{x})}\right)\end{split}$
(14)

where $\mu(\boldsymbol{x})$ and $\sigma(\boldsymbol{x})$ are defined in Eq. 10
and Eq. 11, and $\Phi(z)$ is the cumulative density function (CDF) of the
standard normal random variable $z$. Eq. 14 clearly indicates that
$\alpha_{PI}(\boldsymbol{x})$ favors exploitation because $\Phi(z)$ is
maximized at locations where the predictions are close to $y^{*}$ and have
small uncertainty.

In contrast to PI which discards the magnitude of improvement (regardless of
the magnitude of $y(\boldsymbol{x})$, $I_{PI}(\boldsymbol{x})$ is either $0$
or $1$), EI rewards large improvements over $y^{*}$ by adopting the following
utility function:

$\begin{split}I_{EI}(\boldsymbol{x})=\max(y(\boldsymbol{x})-y^{*},0)\end{split}$
(15)

The corresponding AF can now be obtained by substituting Eq. 15 in Eq. 12 and
using the reparameterization trick (see Section A for the details):

$\begin{split}\alpha_{EI}(\boldsymbol{x})=(\mu(\boldsymbol{x})-y^{*})\Phi(\frac{\mu(\boldsymbol{x})-y^{*}}{\sigma(\boldsymbol{x})})+\sigma(\boldsymbol{x})\phi(\frac{\mu(\boldsymbol{x})-y^{*}}{\sigma(\boldsymbol{x})})\end{split}$
(16)

where $\mu(\boldsymbol{x})$ and $\sigma(\boldsymbol{x})$ are given in Eq. 10
and Eq. 11, respectively, and $\phi(z)$ is the probability density function
(PDF) of $z$. Eq. 16 clearly demonstrates that $\alpha_{EI}(\boldsymbol{x})$
strikes a balance between exploration and exploitation when it is used as the
objective function of the auxiliary optimization problem in Algorithm 1: while
the second term on the right-hand side directly deals with uncertainty and
hence encourages exploration, the first term favors exploitation [37, 38].

Another widely used AF is KG which, unlike PI and EI, is look-ahead because it
chooses $\boldsymbol{x}^{k+1}$ (see Algorithm 1) based on the effect of the
yet-to-be-seen observation (i.e., $y^{k+1}$ which follows a normal
distribution) on the optimum value predicted by the emulator. Following the
terminology and setup of Algorithm 1, this AF quantifies the expected utility
of $\boldsymbol{x}$ at iteration $k+1$ as:

$\begin{split}\alpha_{KG}(\boldsymbol{x})=\mathbb{E}_{p(y\mid\boldsymbol{x},\mathcal{D}^{k})}[\max{\mu^{k+1}}(\boldsymbol{x})]-\max{\mu^{k}}(\boldsymbol{x})\end{split}$
(17)

where $\max{\mu^{k}}(\boldsymbol{x})$ denotes the maximum mean prediction of
the GP trained on $\mathcal{D}^{k}$. The expectation operation in Eq. 17
appears due to the fact that $y^{k+1}$ is not observed yet and
$\alpha_{KG}(\boldsymbol{x})$ is relying on the predictive distribution
provided by the GP that is trained on $\mathcal{D}^{k}$. This expectation
cannot be calculated analytically and hence a Monte Carlo estimate is used in
practice:

$\begin{split}\alpha_{KG}(\boldsymbol{x})\approx\frac{1}{M}\sum_{m=1}^{M}\max{{\mu^{k+1}}^{m}(\boldsymbol{x})}-\max{\mu^{k}}(\boldsymbol{x})\end{split}$
(18)

where $\max{{\mu^{k+1}}^{m}(\boldsymbol{x})}$ is calculated by first drawing a
sample at $\boldsymbol{x}$ from the GP that is trained on $\mathcal{D}^{k}$
and then retraining the GP on $\mathcal{D}^{k}\cup(\boldsymbol{x},y^{m})$
where $y^{m}$ is response of the drawn sample. In practice, a small value must
be chosen for $M$ since maximizing $\alpha_{KG}(\boldsymbol{x})$ over the
input space at each iteration of BO is very expensive. We refer the readers to
[39, 40] for more information on KG and its implementation.

####  2.2.2 Existing Multi-fidelity BO Techniques

As stated in Section 1, the overall computational efficiency of BO can be
increased by leveraging inexpensive LF datasets. MF BO has been successfully
used in many applications such as hyperparameter tuning [41, 42, 43, 44],
finding Pareto fronts in multi-objective optimizations [45, 46, 47], and
solving non-linear state-space models [48, 49]. For MFBO, both the emulator
and the AF must accommodate the multi-source and unbalanced555LF sources
contribute more samples to the training data since they are typically much
cheaper to query from. nature of the data.

Co-Kriging which is an extension of Kriging (or GP) is a popular emulator that
handles MF data by reformulating the covariance function in Eq. 2 as follows
(assuming there are three data sources denoted by $A$, $B$, and $C$):

$\begin{split}\operatorname{cov}([y_{A}(\boldsymbol{x}),y_{B}(\boldsymbol{x}),y_{C}(\boldsymbol{x})]^{T},[y_{A}(\boldsymbol{x}^{\prime}),y_{B}(\boldsymbol{x}^{\prime}),y_{C}(\boldsymbol{x}^{\prime})]^{T})=\Sigma\otimes
r\left(\boldsymbol{x},\boldsymbol{x}^{\prime}\right)\end{split}$ (19)

where $\otimes$ denotes the Kronecker product and $\Sigma$ is a symmetric
positive-definite matrix of size $3\times 3$. This reformulation assumes that
all the responses (regardless of the source) follow a multi-variate normal
distribution and that the matrix $\Sigma_{ij}$ captures the overall
correlation between sources $i$ and $j$666GPs can handle multi-response
datasets in a similar manner, see [50].. While this method can fuse any number
of data sources, it fails to accurately capture cross-source correlations
since the matrix $\Sigma$ has insufficient learning capacity.

Another well-known MF emulation method is that of Kennedy and O’Hagan who fuse
bi-fidelity datasets by learning a discrepancy function that aims to explain
the differences between HF and LF sources. While this bi-fidelity emulator has
proved useful in a wide range of applications [51, 52], it has some major
drawbacks such as inability to jointly learn from more than two sources,
numerical issues, and assuming a priori additive relation between the
discrepancy function and the two data sources.

The bi-fidelity approach of Kennedy and O’Hagan can be viewed as a special
case of hierarchical MF modeling where it is assumed that the relative
accuracy between all the data sources is known. Space mapping techniques
belong to this category, but they are rarely used for sequential sampling, BO,
or MF modeling (see [53] for a bi-fidelity example). These techniques are
typically employed in solving partial differential equations, particularly to
accelerate the convergence of an HF simulation (e.g., based on fine
discretization) by initializing it via the results of an LF simulation.

Upon reformulating the covariance function in Eq. 2, GPs can also be used for
hierarchical MF modeling. For instance, the single-task MF GP of the popular
BoTorch package adopts an additive covariance function that relies on
introducing two user-defined quantitative features [24, 20]. The first
feature, denoted by $x_{a}$, is restricted to the $[0,1]$ range and assigns a
fidelity value to a source based on the user’s belief (larger values
correspond to higher fidelities). The second feature, denoted by $x_{b}$, is
the iteration fidelity parameter and benefits MF BO specifically in the
context of hyperparameter tuning of large machine learning models. The
covariance function directly uses these two additional features as follows:

$\begin{split}\operatorname{cov}(\boldsymbol{x},\boldsymbol{x}^{\prime})=c_{0}(\boldsymbol{x},\boldsymbol{x}^{\prime})+e_{1}(x_{a},x^{\prime}_{a})c_{1}(\boldsymbol{x},\boldsymbol{x}^{\prime})+e_{2}([x_{a},x_{b}]^{T},[x^{\prime}_{a},x^{\prime}_{b}]^{T})c_{2}(\boldsymbol{x},\boldsymbol{x}^{\prime})+e_{3}(x_{b},x^{\prime}_{b})c_{3}(\boldsymbol{x},\boldsymbol{x}^{\prime})\end{split}$
(20)

where $c_{i}(\boldsymbol{x},\boldsymbol{x}^{\prime})$ are Matern kernels777The
parameters of these kernels are endowed with Gamma priors in BoTorch. that
characterize the spatial correlations across the numerical inputs and
$e_{i}(\cdot)$ are user-defined functions that model the cross-source
correlations. $e_{1}(x_{a},x^{\prime}_{a})$ and $e_{3}(x_{b},x^{\prime}_{b})$
are bias kernels that aim to take the discrepancies among the sources into
account while $e_{2}([x_{a},x_{b}]^{T},[x^{\prime}_{a},x^{\prime}_{b}]^{T})$
models the fidelities’ interaction kernel (see Appendix B for more details).

Despite being useful [54], Eq. 20 has some limitations. For instance, it
requires a priori knowledge about the exact hierarchy of fidelities and how
they should be encoded as a numerical feature (i.e., $x_{a}$). Additionally,
the manually-defined functions $e_{i}(\cdot)$ are insufficiently flexible to
learn complex cross-source relations and they also do not provide any learned
metric that quantifies which LF sources are useful for MF BO.

Compared to the few emulators described above, the diversity of the AFs in MF
BO is more since they are often tailored to the application, see [55, 24, 25,
56]. Many of these developments leverage existing AFs that are used in SF BO
such as EI, PI, upper confidence bound [57], Thompson sampling [58], or GP-
predictive entropy search [59]. One specific example is most likely EI (MLEI)
[60] which is tailored to direct policy search problems where the EI in Eq. 16
is first scaled by multiple context-specific priors and then the resulting AFs
are optimized to determine the next candidate point. As another example, Wu et
al [44] develop trace-aware KG to accelerate the hyperparameter tuning process
of machine learning models whose training relies on minimizing the loss
function (defined as the expected prediction error on the validation data). MF
BO is useful in this process since the evaluation accuracy of the loss
function can be controlled by parameters such as the number of iterations and
training/validation data points. Correspondingly, trace-aware KG adjusts these
parameters to use LF but inexpensive evaluations of the loss function (and it
trace) during training. We highlight that EI is also widely used in
hyperparameter tuning problems [61, 62, 63].

##  3 Proposed Approach

As was previously stated, the emulator and AF are the two fundamental
components of any BO framework. In this Section, we first discuss the
rationale for using LMGP as the emulator of our MF BO framework in Section 3.1
and then introduce our novel cost-aware AF in Section 3.2. We elaborate on the
convergence conditions and provide an algorithmic summary of our framework in
Section 3.3.

###  3.1 Multi-fidelity Emulation via LMGP

As schematically illustrated in Figure 2, MF emulation via LMGPs is quite
straightforward: Assuming there are $ds$ data sources, we augment the inputs
with the additional categorical variable
$s=\\{^{\prime}1^{\prime},\cdots,^{\prime}ds^{\prime}\\}$ whose $j^{th}$
element corresponds to data source $j$ for $j=1,\cdots,ds$. After this
augmentation, the inputs and outputs of all the datasets are concatenated as
(following the notation of Figure 2):

$\begin{split}\boldsymbol{U}=\left[\begin{array}[]{cc}\boldsymbol{U}_{1}&{}^{\prime}\mathbf{1}^{\prime}_{n_{1}\times
1}\\\ \boldsymbol{U}_{2}&{}^{\prime}\mathbf{2}^{\prime}_{n_{2}\times 1}\\\
\vdots&\vdots\\\
\boldsymbol{U}_{ds}&{}^{\prime}\mathbf{ds}^{\prime}_{n_{ds}\times
1}\end{array}\right]\text{ and
}\boldsymbol{y}=\left[\begin{array}[]{c}\boldsymbol{y}_{1}\\\
\boldsymbol{y}_{2}\\\ \vdots\\\
\boldsymbol{y}_{ds}\end{array}\right]\end{split}$ (21)

where the subscripts $1,2,...,ds$ correspond to the data sources, $n_{j}$ is
the number of samples obtained from $s(j)$ (i.e., source $j$),
$\boldsymbol{U}_{j}$ and $\boldsymbol{y}_{j}$ are, respectively, the
$n_{j}\times(dx+dt)$ feature matrix and the ${n_{j}\times 1}$ vector of
responses obtained from $s(j)$, and ${}^{\prime}\boldsymbol{j}^{\prime}$ is a
categorical vector of size ${n_{j}\times 1}$ whose elements are all set to
${}^{\prime}j^{\prime}$. Once the $\\{\boldsymbol{U},\boldsymbol{y}\\}$
dataset is built, it is directly fed into LMGP to build an MF emulator.

Figure 2: Multi-fidelity emulation with LMGP: The training data is built by
first augmenting the inputs with the categorical feature $s$ that denotes the
data source of a sample and then concatenating all the inputs and outputs, see
Eq. 21. LMGP can use one or more manifolds to encode the categorical variables
into a quantitative space. For MF emulation, we recommend using two manifolds
to simplify the visualization of cross-source relations: one manifold for $s$
and the other for the rest of the categorical variables, i.e.,
$\boldsymbol{t}$.

We argue that, compared to the existing techniques (see Section 2.2.2), LMGPs
provide a more flexible and accurate mechanism to build MF emulators, see
Figure 3 for a comparison study on an analytic example. This superiority is
because LMGP learns the relations between the sources (which are hidden in the
combined datasets) in a manifold. This manifold is learned nonlinearly and
hence has a much higher representation power than methods that rely on linear
operations, e.g., the matrix $\Sigma$ in co-Kriging that linearly scales the
correlations, see Eq. 19. Similarly, LMGP has major advantages over single-
task GPs (STGPs) reviewed in Section 2.2.2 because $(1)$ it does not assume
any hierarchy across the data sources, $(2)$ the cross-source relations are
encoded via learned latent variables which have significantly higher
representation power than a single user-defined scalar variable (see Eq. 20)
that directly affects the covariance function of the underlying GP and
requires knowledge of the relative source fidelities.

Figure 3: Emulator comparison: We compare the prediction accuracy of LMGP
against single-task multi-fidelity GP (STGP) and Co-Kriging (Co-K) on an
analytic problem with three sources (see Borehole in Table 1 where HF, LF3,
and LF4 are used). While LMGP and STGP use $(5,5,50)$ initial data for (HF,
LF3, LF4) sources, Co-K uses $(50,50,50)$ to achieve comparable performance.
Each emulator is trained $10$ times by randomizing the initial data. It is
evident that LMGP consistently outperform other methods in emulating all the
sources. Prediction accuracy is measured by calculating the mean squared error
(MSE, on $log$ scale) on unseen test data.

As shown in Figure 3, LMGPs build MF emulators that more accurately learn all
the sources (rather than just the HF source). While we can alter the
formulations in Section 2.1 such that LMGPs prioritize learning the HF source,
we do not believe this is a good general decision in the context of MF BO. The
reasoning behind our belief is that the quality of LF predictions (obtained
via the emulator) greatly affects the exploration nature of BO and, hence, its
convergence behavior. The empirical results in Section 4 strongly support this
reasoning.

Another major advantage of LMGPs is that their learned manifold provides an
intuitive and visualizable global metric for comparing the relative
discrepancies/similarities among the data sources, see Figure 6 for four
examples with different number of data sources. This manifold is particularly
useful in detecting anomalous sources whose data adversely affects MF BO. We
demonstrate this in Section 4 with some examples where we use the initial MF
data (i.e., before starting the BO iterations) to correctly predict whether SF
BO outperforms MF BO.

Given the importance of identifying relative discrepancies among data sources,
we slightly modify the correlation function of LMGP to learn two manifolds
where the first one encodes the original categorical variables (collectively
denoted by $\boldsymbol{t}$ in Section 2.1) while the second one encodes the
data source identifier (denoted by $s$ in Figure 2). The new correlation
function is (compare to Eq. 5):

$\displaystyle r(\begin{bmatrix}\boldsymbol{x}\\\ \boldsymbol{t}\\\
s\end{bmatrix},\begin{bmatrix}\boldsymbol{x^{\prime}}\\\
\boldsymbol{t^{\prime}}\\\
s^{\prime}\end{bmatrix})=\exp\\{-\sum_{i=1}^{dx}10^{\omega_{i}}(x_{i}-x_{i}^{\prime})^{2}\
-\sum_{i=1}^{dz}(z_{i}(\boldsymbol{t})-z_{i}(\boldsymbol{t}^{\prime}))^{2}-\sum_{i=1}^{dh}(h_{i}(s)-h_{i}(s^{\prime}))^{2}\\}$
(22)

where $\boldsymbol{h}(s)=[h_{1},\cdots,h_{dh}]^{T}$ is the latent
representation of data source $s$. Similar to Eq. 6, this latent
representation is obtained by post-multiplying a prior vector by the
parametric matrix $\boldsymbol{A}_{h}$:

$\begin{split}\boldsymbol{h}(s)=\boldsymbol{\zeta}(s)\boldsymbol{A}_{h}\end{split}$
(23)

In our studies, we always design the prior by one-hot-encoding the categorical
variable $s=\\{^{\prime}1^{\prime},\cdots,^{\prime}ds^{\prime}\\}$ that
identifies the data source and estimate all of LMGP’s hyperparameters via MLE:

$\begin{split}[\widehat{\boldsymbol{\omega}},\widehat{\boldsymbol{A}},\widehat{\boldsymbol{A}}_{h}]=\underset{\boldsymbol{\omega},\boldsymbol{A},\boldsymbol{A}_{h}}{\operatorname{argmin}}\hskip
8.53581ptn\log\left(\hat{\sigma}^{2}\right)+\log(|\boldsymbol{R}|)\end{split}$
(24)

where the elements of $\boldsymbol{R}$ are now built using Eq. 22.

We now use Eq. 22 to explain the relation between the latent fidelity
representations, i.e., $\boldsymbol{h}(s)$, and the relative fidelity of the
data sources. At the same inputs, the correlation between the estimated
outputs of sources $s$ and $s^{\prime}$ is:

$\displaystyle 0\leq
corr(y_{s}(\boldsymbol{x},\boldsymbol{t}),y_{s^{\prime}}(\boldsymbol{x},\boldsymbol{t}))=r(\begin{bmatrix}\boldsymbol{x}\\\
\boldsymbol{t}\\\ s\end{bmatrix},\begin{bmatrix}\boldsymbol{x}\\\
\boldsymbol{t}\\\
s^{\prime}\end{bmatrix})=\exp\\{0-0-\sum_{i=1}^{d_{h}}(h_{i}(s)-h_{i}(s^{\prime}))^{2}\\}\leq
1$ (25)

So sources with similar fidelities which provide highly correlated responses,
must have similar latent representations, i.e.,
$\boldsymbol{h}_{i}(s)\simeq\boldsymbol{h}_{i}(s^{\prime})$. This relation is
illustrated in Figure 6 where sources with similar fidelities are encoded by
close-by points in the manifold.

###  3.2 Multi-source Cost-aware Acquisition Function

The choice of AF affects the performance of BO quite significantly. This
choice is especially important in MF BO because, in addition to balancing
exploration and exploitation, the AF has to consider the biases of LF data and
source-dependent sampling costs. To demonstrate these challenges, consider the
analytic example in Figure 4(a) where, while the two functions are correlated,
the LF source’s global optimum (the location and the corresponding $y$ value)
is quite different than that of the HF source. Since LF sources are typically
much cheaper than the HF source, a naive AF (that merely scales the expected
utility based on the cost) forces MF BO to converge to the global optimum of
the LF source, see Figure 4(b) for a one-dimensional example.

(a) Biased LF source

(b) Effect of AF on convergence

Figure 4: Double-well potential: The LF source has a systematic bias because
its optima do not match with those of the HF source. However, since the two
curves have similar trends, an effective AF can leverage this correlation to
reduce the overall data acquisition costs. The LF and HF sources have a cost-
per-sample of, respectively, $1$ and $1000$ in this example.

Contrary to existing approaches, we argue that the key to addressing the above
challenges is to quantify the information value of LF and HF data based on
different metrics which are then compared against each other to determine the
candidate input and the corresponding source. In particular, we propose to use
the LF sources exclusively for exploration to leverage their correlations with
the HF source while preventing them from dominating the convergence behavior
of MF BO. Additionally, we propose to exclusively employ the HF source for
exploitation to maximally use its trustworthy samples888These samples may be
corrupted via $\epsilon\sim\mathcal{N}(0,\sigma^{2})$ where $\sigma^{2}$ is
the (unknown) noise variance. during optimization.

To develop the AF for the $j^{th}$ LF source with $j=[1,\cdots,ds]$ and $j\neq
l$ where $l$ denotes the HF source, we follow Section 2.2.1 and define the
improvement (for a maximization problem) as $y_{j}(\boldsymbol{x})–y_{j}^{*}$
where $y_{j}(\boldsymbol{x})$ denotes the LMGP-based prediction at
$\boldsymbol{x}$ for source $j$ and $y_{j}^{*}$ is the best function value in
the obtained dataset from source $j$. We use the reparameterization trick to
rewrite this improvement as:

$\begin{split}\frac{y_{j}(\boldsymbol{x})-\mu_{j}(\boldsymbol{x})}{\sigma_{j}(\boldsymbol{x})}=z\sim\mathcal{N}(0,1)\Rightarrow
y_{j}(\boldsymbol{x})-y_{j}^{*}=(\mu_{j}(\boldsymbol{x})-y_{j}^{*})+\sigma_{j}(\boldsymbol{x})z\end{split}$
(26)

where $\mu_{j}(\boldsymbol{x})$ and $\sigma_{j}(\boldsymbol{x})$ are defined
in Eq. 10 and Eq. 11, respectively. In Eq. 26 the
($\mu_{j}(\boldsymbol{x})-y_{j}^{*}$) and $\sigma_{j}(\boldsymbol{x})z$ terms
control the exploitation and exploration aspects of the improvement,
respectively. We now define our utility function that focuses on exploration
by dropping the first term on the far right-hand-side of Eq. 26:

$\begin{split}I_{LF}(\boldsymbol{x};j)=\begin{cases}\sigma_{j}(\boldsymbol{x})z&y_{j}(\boldsymbol{x})>y_{j}^{*}\\\
0&y_{j}(\boldsymbol{x})\leq y_{j}^{*}\end{cases}\end{split}$ (27)

which is used for the $j^{th}$ LF source in our framework. We obtain the
corresponding AF by substituting $I_{LF}(\boldsymbol{x};j)$ in Eq. 12:

$\begin{split}\alpha_{LF}(\boldsymbol{x};j)=\int_{-\infty}^{\infty}I_{LF}(\boldsymbol{x};j)\phi(z)\mathrm{d}z=\int_{-\infty}^{\infty}\sigma_{j}(\boldsymbol{x})z\phi(z)\mathrm{d}z\end{split}$
(28)

The integral is zero for $y_{j}(\boldsymbol{x})<y_{j}^{*}$ so we find the
corresponding switch point in terms of $z$:

$\begin{split}y_{j}(\boldsymbol{x})=y_{j}^{*}\Rightarrow\mu_{j}(\boldsymbol{x})+\sigma_{j}(\boldsymbol{x})z=y_{j}^{*}\Rightarrow
z_{0}=\frac{y_{j}^{*}-\mu_{j}(\boldsymbol{x})}{\sigma_{j}(\boldsymbol{x})}\end{split}$
(29)

Inserting Eq. 29 in Eq. 28 yields:

$\begin{split}\begin{aligned}
&\alpha_{LF}(\boldsymbol{x};j)=\int_{z_{0}}^{\infty}\sigma_{j}(\boldsymbol{x})z\phi(z)\mathrm{d}z=\int_{z_{0}}^{\infty}\frac{\sigma_{j}(\boldsymbol{x})z}{\sqrt{2\pi}}e^{-\frac{z^{2}}{2}}\mathrm{d}z\\\
&=\frac{\sigma_{j}(\boldsymbol{x})}{\sqrt{2\pi}}\int_{z_{0}}^{\infty}ze^{-\frac{z^{2}}{2}}\mathrm{d}z=\frac{\sigma_{j}(\boldsymbol{x})}{\sqrt{2\pi}}\int_{z_{0}}^{\infty}(e^{-\frac{z^{2}}{2}})^{\prime}\mathrm{d}z\\\
&=-\frac{\sigma_{j}(\boldsymbol{x})}{\sqrt{2\pi}}[e^{-\frac{z^{2}}{2}}]_{z_{0}}^{\infty}=\sigma_{j}(\boldsymbol{x})\phi(z_{0})=\sigma_{j}(\boldsymbol{x})\phi(\frac{y_{j}^{*}-\mu_{j}(\boldsymbol{x})}{\sigma_{j}(\boldsymbol{x})})\end{aligned}\end{split}$
(30)

Comparison between this AF and Eq. 16 illustrates that our proposed AF for LF
sources is the same as the exploration part of EI.

For the HF source, we propose to use PI as the AF because it focuses on
exploitation and is computationally efficient (the efficiency is due to the
analytic form of PI which dispenses with expensive numerical integration).
Assuming source $l$ provides the HF data, this AF is given by:

$\begin{split}\alpha_{HF}(\boldsymbol{x};l)=\Phi(\frac{\mu_{l}(\boldsymbol{x})-y_{l}^{*}}{\sigma_{l}(\boldsymbol{x})})\end{split}$
(31)

We use $\alpha_{HF}(\boldsymbol{x};l)$ and $\alpha_{LF}(\boldsymbol{x};j)$ as
defined above in each iteration of BO to solve $ds$ auxiliary optimization
problems (assuming there are $ds$ data sources) to find the candidate location
with the highest expected utility from each source. We then find the final
candidate point and the corresponding source by comparing the cost-normalized
version of these $ds$ values. Hence, assuming source $l$ provides HF data
among the $ds$ sources, our AF that considers source-dependent costs and
fidelities is:

$\alpha_{MFCA}(\boldsymbol{x};j)=\begin{cases}\nicefrac{{\alpha_{LF}(\boldsymbol{x};j)}}{{O(j)}}&j=[1,\cdots,ds]\And
j\neq l\\\ \nicefrac{{\alpha_{HF}(\boldsymbol{x};l)}}{{O(l)}}&j=l\end{cases}$
(32)

and

$[\boldsymbol{x}^{k+1},j^{k+1}]=\underset{\boldsymbol{x},j}{\operatorname*{argmax}}\hskip
5.69054pt\alpha_{MFCA}(\boldsymbol{x};j)$ (33)

where $O(j)$ is the cost associated with taking a single sample from source
$j$ and $\boldsymbol{x}^{k+1}$ is the point that source
$j^{k+1}\in\\{1,\cdots,ds\\}$ must be evaluated in the current BO iteration.
We now point out several important aspects of Eq. 32.

The AFs in Eq. 30 and Eq. 31 quantify the value of a sample by comparing it to
the best available sample for the corresponding source (and not the best
sample across all the sources). The advantage of this source-wise comparison
in each of the AFs is that it encourages sampling from sources that provide
larger values (which is desirable for a maximization problem). However, this
formulation enables LF sources whose optima are larger than the HF source to
dominate the optimization process where not only more samples are taken from
these LF sources (which may also cause numerical issues), but also the
converged solution does not belong to the profile of the HF source. This issue
is exacerbated once the AFs are divided by the data collection costs (see Eq.
32) since LF samples are (typically) much cheaper than HF data.

The abovementioned issues are addressed with three mechanisms in our MFCA
approach. Firstly, we always report the points sampled from the HF source as
the final optimization history (this choice ensures that the final solution
indeed belongs to the profile of the HF source but it does not guarantee
global optimality999BO does not guarantee global optimality anyways.).
Secondly, we always use the fidelity manifold of the LMGP that is trained on
the initial data to detect the LF sources that should not be used in BO due to
their severe discrepancy (we demonstrate the benefit of this exclusion with
examples in Section 4). Thirdly, we have designed the core of our AFs for the
HF and LF sources based on, respectively, the CDF and the scaled PDF of the
standard normal variable $z$. As detailed below and empirically shown in
Section 4, the intricate relation between these two functions during the
optimization reduces the effect of LF sources’ biases on the convergence.

As illustrated in Figure 5, $\phi(z)$ and $\Phi(z)$ have comparable values up
to $\mathbb{E}[z]=0$ but the ratio $\Phi(z)/\phi(z)$ increasingly grows as $z$
realizes larger values. This trend indicates that if an HF candidate point
sufficiently improves $y_{l}^{*}$, then Eq. 33 queries the HF source at that
point to obtain a new sample for the next BO iteration. The frequency of this
query during the optimization process is controlled by the data collection
cost and $\sigma_{j}(\boldsymbol{x})$, see Eq. 30 and Eq. 31, respectively. If
HF samples are highly costly, the auxiliary optimization in Eq. 33 reduces the
sampling frequency. However, unlike existing approaches such as BoTorch, this
reduction does not translate into ”never sampling from the HF source” (even if
the cost ratios are as large as $1000$, see Section 4) because
$\Phi(z)/\phi(z)$ can be quite large. Regarding the
$\sigma_{j}(\boldsymbol{x})$ term in Eq. 30, we note that it encourages
exploring the regions where LMGP provides highly uncertain predictions for an
LF source. This scenario happens when an LF source is rarely sampled and there
are insufficient correlations between that source and other sources.

In summary, our proposed AF in Eq. 32, while involving intricate interactions
between the fidelities and costs, has a simple form which is analytic (and
hence computationally efficient) and interpretable. As we illustrate in
Section 4 this AF, combined with LMGP, dramatically improves the performance
of our MFCA BO framework.

Figure 5: Standard normal variable: The PDF and CDF of $z\sim\mathcal{N}(0,1)$
have comparable magnitudes up to the mean of $z$ (i.e., $0$) but increasingly
differ after $0$.

###  3.3 Convergence Metric

Similar to AFs, convergence criteria of MF BO techniques are traditionally
tailored to the application since many factors (e.g., budget constraints,
numerical issues, or convergence history) affect the results. We believe the
emulator and AF of our framework alleviate many of the convergence issues
associated with MF BO and hence use two simple convergence metrics: overall
costs and maximum number of iterations without improvement. The former is a
rather generic metric but it can result in considerably high number of
iterations if one of the LF sources is extremely inexpensive to query. For
this reason, we recommend using additional metrics (such as the second one
above) that track convergence.

Algorithm 2 summarizes our framework for MFCA BO. Compared to Algorithm 1, the
major differences are in the choice of the emulator and AF which now can
handle multi-source data that have different levels of fidelity and cost. In
addition, a pre-processing step is added which leverages the fidelity manifold
of LMGP to detect the LF sources that must be excluded from the BO process due
to their large discrepancies with respect to the HF source.

Given: Initial multi-fidelity data
$\mathcal{D}^{k}=\\{(\boldsymbol{x}^{i},y^{i}\\}_{i=1}^{k}$, black-box
functions $f(\boldsymbol{x};j)$ and their corresponding sampling costs $O(j)$
where $j=[1,...,ds]$

Goal: Optimizing high-fidelity function (source $l\in[1,...,ds]$)

Define: Utility functions (see Eq. 13 and Eq. 27) and stop conditions

Step $0$: Train an LMGP and exclude highly biased low-fidelity sources based
on its fidelity manifold

while  _stop conditions not met_ do

  1. 1.

Train an LMGP using $\mathcal{D}^{k}$

  2. 2.

Define the multi-fidelity cost-aware acquisition function (see Equations 30,
31, and 32):

$\alpha_{MFCA}(\boldsymbol{x};j)=\begin{cases}\alpha_{LF}(\boldsymbol{x};j)/O(j)&j\in\\{1,\cdots,ds\\}\quad\&\quad
j\neq l\\\ \alpha_{HF}(\boldsymbol{x};l)/O(l)&j=l\end{cases}$

* 3.

Solve the auxiliary optimization problem:

$[\boldsymbol{x}^{k+1},j^{k+1}]=\underset{\boldsymbol{x},j}{\operatorname{argmax}}\hskip
5.69054pt\alpha_{MFCA}(\boldsymbol{x};j)$

* 4.

Query $f(\boldsymbol{x};j)$ at point $\boldsymbol{x}^{k+1}$ from source
$j^{k+1}$ to obtain $y^{k+1}$ 5. Update data:
$\mathcal{D}^{k+1}\leftarrow\mathcal{D}^{k}\cup(\boldsymbol{x}^{k+1},y^{k+1})$
6. Update counter: $k\leftarrow k+1$

end while

Output: Updated data
$\mathcal{D}^{k}=\\{(\boldsymbol{x}^{i},y^{i})\\}_{i=1}^{k}$, LMGP

Algorithm 2 Multi-fidelity Cost-aware Bayesian Optimization for Maximization

##  4 Results and Discussions

We compare the performance of the following four methods on four analytic and
two real-world examples (Sections 4.1 and 4.2, respectively):

  * •

LMGPCA: Our proposed MFCA BO approach.

  * •

LMGPEI: Single-fidelity BO whose emulator and AF are LMGP and EI,
respectively.

  * •

LMGPPI: Single-fidelity BO whose emulator and AF are LMGP and PI,
respectively.

  * •

BoTorch: Multi-fidelity BO with BoTorch where STGP and KG are used as emulator
and AF.

The first three methods are myopic and hence computationally much faster than
BoTorch which is lookahead (see Section 2.2.2 for details on BoTorch).
Assuming the computational costs of BO (which are mainly associated with
emulation and solving the auxiliary optimization problem) are negligible
compared to querying any of the data sources, we compare the above methods in
terms of cost and accuracy which are defined as, respectively, the overall
data collection costs and the ability to find the global optimum of the HF
source.

Our rationales for comparing LMGPCA against the above three methods are to
demonstrate: $(1)$ the advantages of employing LF sources in BO, $(2)$ that
our designed myopic AF (see Eq. 32) improves both the sampling cost and
accuracy compared to even lookahead methods such as BoTorch, and $(3)$ the
importance of excluding highly biased LF sources from BO. We also note that
the emulators in BoTorch cannot handle categorical inputs and hence we compare
LMGPCA to LMGPEI and LMGPPI in Section 4.2 where the two examples have
categorical variables.

For all the methods, we terminate the optimization process if either of the
following conditions are met: $(1)$ the overall sampling cost reaches a pre-
determined maximum budget, or $(2)$ the best HF sample (i.e., $y_{l}^{*}$ in
Eq. 31) does not improve over $50$ consecutive iterations. These conditions
are quite simple and straightforward; allowing us to focus on the effects of
AF and the emulator on the performance.

In Section 4.1, we set the maximum budget to $40000$ units for LMGPCA, LMGPEI,
and LMGPPI while for BoTorch we choose $50000$ since it is, as demonstrated
below, highly inefficient and requires more samples to provide reasonable
accuracy. Additionally, in all of our examples an HF sample costs $1000$ so
LMGPEI and LMGPPI are terminated based on the maximum budget condition. In
Section 4.1, we set the maximum budget to $1000$ and $1800$ for the two
examples since their data collection cost are much lower than the examples in
Section 4.1.

###  4.1 Analytic Examples

As detailed in Table 1 in Appendix C, we consider four analytic examples
(Double-well Potential, Rosenbrock, Borehole, and Wing) whose input
dimensionality ranges from $1$ to $10$. All examples are single-response and
the number of data sources varies between $2$ and $5$. The source-dependent
sampling costs and number of initial data points are also detailed in Table 1.
To compare the robustness and effectiveness of the four BO methods described
above, we use relatively small initial datasets (especially from the HF
source) and consider various cost ratios (the maximum cost ratio between two
sources equals $1000$). For each example, we quantify the effect of random
initial data by repeating the optimization process $20$ times for each method
(all initial data are generated via Sobol sequence).

(a) Double-well Potential

(b) Rosenbrock

(c) Borehole

(d) Wing

Figure 6: Fidelity manifolds of analytic examples: We train LMGPs on the
initial multi-source data to (inversely) learn the relation between the
sources. The encoding described in Section 3.1 is used and hence each source
is represented with a point in the fidelity manifold of LMGP. For each
example, the plot corresponds to a randomly selected case among the $20$
repetitions (only one plot is shown due to consistency across the
repetitions). Based on the distances in the fidelity manifold (see also Eq.
22), we conclude that the second and third sources are highly biased in
Borehole and hence must be excluded from the MF BO process.

Table 1 also enumerates the relative accuracy of the LF sources of each
example by calculating the relative root mean squared error (RRMSE) between
them and the corresponding HF source based on $1000$ samples (these RRMSEs are
not used in BO in any way). In the case of Borehole we observe that, unlike
Wing, the source ID, true fidelity level (based on the RRMSEs), and sampling
cost are not related. For instance, the first LF source is the least accurate
and most expensive among all the LF sources in Borehole.

Per Step $0$ in Algorithm 2, we always use the initial data to train an LMGP
to identify the useful LF sources. Based on Figure 6, we expect the LF sources
to be beneficial in all the examples except for Borehole since the latent
points of the first and second sources are distant from the latent position of
the HF source (we test this expectation in Appendix D). Hence, hereafter we
exclude these highly biased LF sources, i.e., both LMGPCA and BoTorch use
three sources (HF, LF3, and LF4 in Table 1) to optimize Borehole. We note that
BoTorch does not provide any mechanism to detect highly biased LF sources, but
we also leverage this LMGP-based insight to boost the performance of BoTorch
and, in turn, better highlight the effectiveness of our proposed AF.

(a) Double-well Potential

(b) Rosenbrock

(c) Borehole

(d) Wing

Figure 7: Convergence histories: The plots illustrate the best HF sample
(i.e., $y_{l}^{*}$ in Eq. 31) found by each method as a function of sampling
costs accumulated during the BO iterations (the cost of initial data is
included). LMGPCA consistently outperforms other methods in all the examples,
especially in high-dimensional cases (i.e., Borehole and Wing). The solid
thick curves indicate the average behavior across the $20$ repetitions (the
variations associated with BoTorch in Double-well Potential are extremely
small and hence not visible). LMGPCA and BoTorch use three sources in
Borehole. The second termination condition (i.e., maximum of $50$ BO
iterations without improvement in $y_{l}^{*}$) is disabled for BoTorch in Wing
to illustrate its convergence trajectory.

Figure 7 summarizes the convergence histories by tracking the best HF estimate
found by each method (i.e., $y_{l}^{*}$ in Eq. 31) as a function of
accumulated sampling costs. It is evident that LMGPCA consistently outperforms
all other methods across the four examples. In particular, LMGPCA demonstrates
the advantage of leveraging inexpensive LF sources in BO by accelerating the
convergence without sacrificing the accuracy (compare LMGPCA to LMGPEI and
LMGPPIfor any of the examples in Figure 7). Additionally, unlike BoTorch, the
performance of LMGPCA is robust to the input dimensionality and sampling
costs. For instance, BoTorch estimates the optimum as $y_{l}^{*}=10.02$ in
Figure 7(c) while the ground truth is $3.98$. The reason behind this
inaccuracy is that BoTorch fails to find an HF sample whose information value
is large enough to overcome its high sampling cost and, as a result, cheap LF
sources are largely queried. However, these queries do not improve $y_{l}^{*}$
and hence the second strop condition (maximum number of repetitions without
improvement of $y_{l}^{*}$) terminates BoTorch after $50$ iterations. Removing
this stop condition, while significantly increasing the number of iterations,
does not improve the accuracy of BoTorch. To demonstrate the effect of this
removal, we only consider the maximum budget constraint for BoTorch in Wing
and observe that the ground truth is again not found, see Figure 7(d). These
issues are resolved in LMGPCA by the intricate interplay between $\Phi(z)$ and
$\phi(z)$ as explained in Section 3.2.

(a) Double-well Potential

(b) Rosenbrock

(c) Borehole

(d) Wing

Figure 8: Accumulated costs before improvements plateau: The box-plots
illustrate the accumulated costs up to and including the iteration at which
the best HF sample is first obtained (i.e., these box-plots do not consider
the two termination criteria). On the right axis, the converged solution
(averaged across the $20$ repetitions) and ground truth are demonstrated via,
respectively, the blue marker and the horizontal dashed line. In all four
examples, LMGPCA finds the optimum faster than other methods. Comparison
between the axes indicates that small accumulated costs do not imply superior
performance since the converged solution might be a local optimum, as is the
case for BoTorch in 8(c).

To exclude the effect of the termination criteria from the results, we provide
the accumulated cost up to and including the iteration at which each method
finds its best estimate, see Figure 8. In terms of finding the true optimum
(compare the blue dots to the horizontal dashed line), LMGPCAoutperforms all
other methods and is followed by LMGPEI and then LMGPPI(the high accuracy of
the SF methods in finding the ground truth is expected since they only sample
from the HF source and neither of the two termination criteria are stringent).
However, BoTorch either terminates at an incorrect solution or is even
costlier than SF methods. This poor performance is expected since, even though
BoTorch is lookahead, its AF cannot handle the high cost-ratios across the
sources and its emulator does not effectively learn the nonlinear relations
between the HF and LF sources.

(a) Double-well Potential

(b) Rosenbrock

(c) Borehole

(d) Wing

Figure 9: Number of iterations at convergence: As expected, BoTorch and LMGPCA
need more iterations to converge compared to the SF methods. However, the
difference in the case of LMGPCA is very small as our method is quite
efficient in leveraging the LF sources. It is noted that, since only one
sample is obtained per iteration, these plots are also representative of the
total number of samples collected via BO.

To provide more insight into the mechanics of MF BO methods, we also report
the number of iteration at convergence, see Figure 9. As expected, BoTorch and
LMGPCA require more iterations as they aim to leverage cheap LF sources to
reduce the overall costs. LMGPCA needs fewer iterations for convergence than
BoTorch except in Figure 9(b). This behavior is a result of the termination
criteria: the ground truth of Rosenbrock is $-456.3$ while LMGPCA rapidly
converges to a very close solution ($-456$) and then takes highly cheap LF
samples to improve the best sample. These improvements are quite small ($0.01$
per iteration) and hence many iterations are needed for convergence.

(a) Double-well Potential

(b) Rosenbrock

(c) Borehole

(d) Wing

Figure 10: Number of samples taken from each source: As opposed to BoTorch,
LMGPCA optimally and automatically adjust the sampling frequency from each
source. For instance, BoTorch does not sample from two sources in Borehole
since they are much more expensive than the second LF source. However, LMGPCA
not only samples from all sources, but also adjusts the sampling frequency
from the LF sources based on their relative accuracy, initial data, and cost
(note that LF2 is ten times cheaper to query than LF1 in Borehole).

We highlight that, as long as the associated costs are low, a large number of
iterations is not reflective of bad performance since the goal of MF BO is to
reduce the overall sampling costs and not necessarily the total number of
samples. To demonstrate this, we provide the per-source sampling frequencies
for LMGPCA and BoTorch in Figure 10 which demonstrates that LMGPCA
automatically adjusts its sampling mechanism based on the initial data and the
relative accuracy of the LF sources (compared to the corresponding HF source)
and their costs. For instance, unlike BoTorch which avoids querying the HF
source in Borehole, LMGPCA leverages all sources where the number of samples
taken from each source depends on its cost and (in the case of LF sources)
relative accuracy. In particular, we observe that LMGPCA samples almost
equally from the LF sources even though the second source is $10$ times
cheaper (note that these LF sources correspond to the third and fourth sources
in Table 1). This behavior may seem undesirable at the first glance especially
since the LF sources have the same relative accuracy (see RRMSEs in Table 1)
but a closer look indicates that it is primarily caused by the number of
initial samples: since there are $10$ times more data points from the second
LF source, the emulator of LMGPCA correctly provides large prediction
uncertainties which, in turn, results in a large expected utility for the
first LF source, see Eq. 30. These discussions also hold for Wing, see Figure
10(d) where, unlike BoTorch, LMGPCA samples from the LF sources based on their
relative accuracy (which is learnt internally by its emulator) as well as
cost.

Finally, we demonstrate the performance of LMGPCA in balancing exploration and
exploitation. To this end, source-wise sampling orders made by LMGPCA and
BoTorch are visualized for a randomly selected repetition in each example, see
Figure 11. As it can be observed, LMGPCA alternates between all the sources:
in each example the majority of the samples are from the LF sources which are
queried based on exploration (see Eq. 30 and Eq. 32) while the expensive HF
source is typically used with much lower frequency. The sampling orders are
particularly interesting for Rosenbrock and Borehole where, unlike BoTorch
which struggles to alternate between the sources, LMGPCA effectively uses all
sources during the optimization.

Figure 11: Source-wise sampling orders: A repetition is randomly selected from
each example to visualize the sampling orders made by LMGPCA and BoTorch. The
horizontal axis enumerates the number of BO iterations while the vertical axis
denotes the example and the MF BO method used for optimization. There are
$2,2,3\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ 4$ data sources in
each example from top to bottom. Unlike BoTorch, LMGPCA balances exploration
and exploitation throughout the optimization process.

###  4.2 Real-world Datasets

In this Section, we study two materials design problems where the goal is to
identify the composition that optimizes the property of interest. Unlike the
examples in Section 4.1, these two problems are noisy and have categorical
inputs. We compare the performance of LMGPCA only against the SF methods since
BoTorch does not accommodate categorical inputs. In both examples, the HF and
LF data are obtained via simulations (based on the density functional theory)
which have different fidelity levels.

The first problem is bi-fidelity and aims to find the member of the
nanolaminate ternary alloy (NTA) family with the largest bulk modulus [64].
The HF and LF datasets each have $224$ samples, one response, and $10$
features ($7$ quantitative and $3$ categorical where the latter have $10$,
$12$, and $2$ levels). In our studies, we presume a cost ratio of $10$ to $1$
between the sources and proceed as follows: we exclude the composition with
the largest bulk modulus from the HF dataset, take $20$ and $10$ samples from,
respectively, the HF and LF datasets (SF methods only use HF samples), and
then initiate the BO methods. We repeat this process $15$ times for each BO
method to quantify its robustness to the random initial data.

Our second problem concerns hybrid organic–inorganic perovskite (HOIP)
crystals where the goal is to find the compound with the smallest inter-
molecular binding energy. There are one HF and two LF datasets in this problem
and their sampling costs are $15$, $10$, and $5$, respectively101010We assign
these sampling costs randomly to the LF sources as we do not know which one is
more accurate.. The three datasets have similar dimensionalities ($1$ output
and $3$ categorical inputs with $10$, $3$, and $16$ levels) but are of
different sizes: the HF dataset has $480$ samples, while the first and second
LF sources have $179$ and $240$ samples, respectively. We apply the three BO
methods to this problem as follows: we first exclude the best compound from
the HF dataset and build the initial MF data by randomly sampling from the
three datasets. Then, we launch the BO process where LMGPEI and LMGPPI only
use the HF samples. We set the size of the initial datasets to $(15,20,15)$
for the HF and LF sources, respectively, and repeat the BO process $15$ times
to assess the repetition-size variability.

(a) NTA (initial data)

(b) HOIP (initial data)

(c) NTA (all data)

(d) HOIP (all data)

Figure 12: Fidelity manifolds: The manifolds in the top row are learned via
the initial data and indicate that the LF sources should not be used in BO for
HOIP because the latent points corresponding to them are positioned far from
the point encoding the HF source. The manifolds in the second row are built
using the entire MF data in each example. The similarity between the two
fidelity manifolds of each example indicates that LMGP can effectively learn
source-wise discrepancies via small data.

Per Step $0$ in Algorithm 2, we train an LMGP to the initial data in each
problem to determine whether any of the LF sources must be excluded from BO.
As demonstrated in Figure 12(a), the LF and HF sources in NTA are highly
correlated since their corresponding latent points are very close in the
learnt fidelity manifold of LMGP. However, the latent points in Figure 12(b)
are quite distant and hence we exclude both LF sources from the BO process. It
is noted that $(1)$ we provide these manifolds for a randomly selected
repetition in each example since they insignificantly change across the
repetitions (most changes are due to rotation and translation of all the
latent points which do not affect the relative distances), and $(2)$ even
though small initial data is used in training the LMGPs, the resulting
manifolds provide trustworthy representations of the relative fidelities. To
test this second point, we fit an LMGP to the entire data in each example and
visualize the resulting manifolds, see Figures 12(c) and 12(d). As it can be
observed, while the manifolds do not match exactly, the relative distances
between the latent points are similar.

Figure 13 summarizes the convergence histories by tracking the best HF
estimate found by each method (i.e., $y_{l}^{*}$ in Eq. 31) as a function of
the accumulated costs. As expected, LMGPCA outperforms the two SF methods in
NTA but not in HOIP. In Figure 13(a) we observe that LMGPEI and LMGPPI cannot
find the optimum compound before the convergence criteria terminate the
optimization process. However, these two methods perform quite well in HOIP
and converge to a value that is very close to the ground truth (the small
difference can be eliminated by relaxing the convergence metrics). As
expected, LMGPCA finds a sub-optimal compound in HOIP since the highly biased
LF sources steer the search process in the wrong direction.

(a) NTA

(b) HOIP

Figure 13: Convergence history: The plots illustrate the best HF sample (i.e.,
$y_{l}^{*}$ in Eq. 31) found by each method as a function of sampling costs
accumulated during the BO iterations (the cost of initial data is included).
As expected, LMGPCA outperforms the single-fidelity methods only in NTA since
the LF sources of HOIP have major discrepancies. The solid thick curves
indicate the average behavior across the $20$ repetitions.

(a) NTA

(b) HOIP

Figure 14: Accumulated costs before improvements plateau: The box-plots
illustrate the accumulated costs up to and including the iteration at which
the best HF sample is first obtained (i.e., these box-plots do not consider
termination criteria). On the right axis, the converged solution (averaged
across the $20$ repetitions) and ground truth are demonstrated via,
respectively, the blue marker and the horizontal dashed line.

Similar to Section 4.1 we also provide the accumulated cost up to and
including the iteration at which each method finds its best compound (which
may not correspond to the ground truth) in each example, see Figure 14. In the
case of NTA, LMGPCA outperforms both LMGPEI and LMGPPI in terms of both
accuracy (i.e., finding the ground truth - compare the blue dots to the
horizontal dashed line) and consistency (i.e., showing small variations across
the repetitions - compare the box heights). In the case of HOIP, however,
LMGPCA provides lower accuracy than the SF methods since it is using highly
biased LF sources. Even though LMGPCA is more robust to variations in the
initial data, the lower accuracy does not justify its use for HOIP.

(a) NTA

(b) HOIP

(c) NTA

(d) HOIP

(e) HOIP

Figure 15: Sampling behavior of single- and multi-fidelity BO: In these
examples LMGPCA takes no more iterations (i.e., total number of additional
samples of any fidelity) than either LMGPPI or LMGPEI. In both examples, most
of these samples are from LF datasets which is desirable in MF BO as long as
the LF sources are sufficiently correlated with the HF source.

We now investigate the resource allocation behavior of LMGPCA. As shown in
Figures 15(a) and 15(b), LMGPCA takes equal or fewer iterations to converge
(note that since one sample is taken per iteration, this means that LMGPCA
takes fewer overall samples). In the case of NTA, this behavior is desirable
especially since most samples are taken from the LF source which is cheaper to
query, see Figure 15(c). However, in the case of HOIP, this seemingly
desirable behavior results in convergence to an incorrect solution. Hence, we
emphasize the importance of Step $0$ in Algorithm 2: while LMGPCA can
effectively allocate resources based on the initial dataset sizes and data
collection costs (see Figures 15(d) and 15(e)), highly biased LF sources can
steer the search in the wrong direction and, in turn, result in convergence to
an incorrect solution.

##  5 Conclusion

We introduce a multi-fidelity cost-aware framework for Bayesian optimization
of expensive black-box functions. Compared to single-source BO, our framework
provides improved accuracy and convergence rate by leveraging inexpensive LF
sources during the optimization. Unlike existing MF BO techniques, our method
accommodates an arbitrary number of LF sources and can effectively balance
exploration and exploitation regarding both the search space and source
utilization. We demonstrate these benefits on both analytic and engineering
examples and argue that they are the results of our new acquisition function
as well as integrating LMGPs with BO.

One of the major outcomes of our work is determining (only via the initial
data) if using LF sources in BO improves the performance. Currently, we make
this decision by inspecting the learnt fidelity manifold of LMGP: if the point
representing an LF source is far from the point which encodes the HF source,
then that LF source should not be used in MF BO. This distance is directly
related to the global correlation between an LF and the HF sources and we use
this relation to judge whether the discrepancy is large enough. While this
simple approach works quite well, it may provide sub-optimal results and hence
we plan to improve it in two major directions. Firstly, we envision developing
a local metric which enables LF sources to contribute to BO even if they are
only correlated with the HF source on a small portion of the search space.
Secondly, we plan to integrate the fidelity metrics with the AFs to scale the
information values based on the sample fidelity. With these additions, all LF
sources are kept in the loop since they may provide locally useful
predictions.

Our new AF does not have any calibration parameters but one can certainty
scale its individual components to prioritize (based on, e.g., prior
knowledge) sampling from specific sources. There is also potential in
designing new utility functions that, in addition to (in lieu of) expected and
probability of improvement, are inspired by other AFs such as upper confidence
bound. The examples of this papers do not explore these options since we
observe high performance (which is much better than the competing methods). In
addition, in our studies we use a very simple mechanism for encoding the
fidelity via LMGP and assume the data collection costs are given and fixed but
these choices and parameters can be adjusted based on the application. We plan
to study these directions in our future works.

## Acknowledgments

We appreciate the support from National Science Foundation (award numbers
OAC-2211908 and OAC-2103708) and the Early Career Faculty grant from NASA’s
Space Technology Research Grants Program (award number 80NSSC21K1809).

## Code and Data Availability

The datasets corresponding to analytic and engineering examples are available
on on this GitLab repository. The source codes for reproducing the results of
BoTorch are also available here.

## Appendices

##  Appendix A Formulation of EI and PI

We first derive the AF for PI and then follow a similar procedure for EI. We
insert PI’s utility function in Eq. 12:

$\begin{split}\alpha_{\mathrm{PI}}(\boldsymbol{x})=\mathbb{E}[I(\boldsymbol{x})\mid\mathcal{D}]=\int_{-\infty}^{\infty}\operatorname{Pr}(I(\boldsymbol{x}))I(\boldsymbol{x})\mathrm{d}x\end{split}$
(A-1)

Eq. 13 demonstrates that $I(\boldsymbol{x})$ is zero for
$y(\boldsymbol{x})<y^{*}$, so:

$\begin{split}\alpha_{\mathrm{PI}}(\boldsymbol{x})=\mathbb{E}[I(\boldsymbol{x})\mid\mathcal{D}]=\int_{y^{*}}^{\infty}\operatorname{Pr}(I(\boldsymbol{x})>0)\mathrm{d}y\end{split}$
(A-2)

Hence, to calculate $\alpha_{PI}$ we only need to find
$\operatorname{Pr}(I(\boldsymbol{x})>0)$ which is a function of the random
variable $y(\boldsymbol{x})$. In a GP, the response $y(\boldsymbol{x})$
follows a normal distribution with mean $\mu{(\boldsymbol{x})}$ and variance
$\sigma^{2}(\boldsymbol{x})$:

$\begin{split}y(\boldsymbol{x})\sim\mathcal{N}\left(\mu(\boldsymbol{x}),\sigma^{2}(\boldsymbol{x})\right)\end{split}$
(A-3)

We now apply the reparameterization trick to $y(\boldsymbol{x})$ to calculate
PI. Considering $z\sim\mathcal{N}(0,1)\text{, then
}y(\boldsymbol{x})=\mu(\boldsymbol{x})+\sigma(\boldsymbol{x})z$ is a normal
distribution with mean $\mu{\boldsymbol{(x)}}$ and variance
$\sigma^{2}(\boldsymbol{x})$. Then:

$\begin{split}\operatorname{Pr}(I(\boldsymbol{x})>0)\Leftrightarrow\operatorname{Pr}(y^{*}<y(\boldsymbol{x}))=\operatorname{Pr}(\frac{y^{*}-\mu(\boldsymbol{x})}{\sigma(\boldsymbol{x})}<\frac{y(\boldsymbol{x})-\mu(\boldsymbol{x})}{\sigma(\boldsymbol{x})})\end{split}$
(A-4)

defining $z_{0}=\frac{y^{*}-\mu(\boldsymbol{x})}{\sigma(\boldsymbol{x})}$
which follows $\mathcal{N}\left(0,\sigma^{2}\right)$ simplifies the above as:

$\begin{split}\operatorname{Pr}(\frac{y^{*}-\mu(\boldsymbol{x})}{\sigma(\boldsymbol{x})}<z)=1-\operatorname{Pr}(z\leq\frac{y^{*}-\mu(\boldsymbol{x})}{\sigma(\boldsymbol{x})})=1-\Phi(z)=\Phi(-z)=\Phi(\frac{\mu(\boldsymbol{x})-y^{*}}{\sigma(\boldsymbol{x})})\end{split}$
(A-5)

where $\Phi$ is the CDF of the standard normal variable [35, 36].

In the case of EI, we follow similarly and use the reparameterization trick to
define $y(\boldsymbol{x})=\mu(\boldsymbol{x})+\sigma(\boldsymbol{x})z$ to
rewrite Eq. 15 as
$I(\boldsymbol{x})=\mu(\boldsymbol{x})+\sigma(\boldsymbol{x})z-y^{*}$ where
$z$ is the standard normal random variable. We now insert this utility
function into Eq. 12:

$\begin{split}\alpha_{EI}(\boldsymbol{x})\equiv\mathbb{E}[I(\boldsymbol{x})|D]=\int_{-\infty}^{\infty}\max(\mu(\boldsymbol{x})+\sigma(\boldsymbol{x})z-y^{*},0)\phi{(z)}\mathrm{d}z\end{split}$
(A-6)

To eliminate the $max$ operator and simplify the integration, we divide
$\mu(\boldsymbol{x})+\sigma(\boldsymbol{x})z-y^{*}$ into two negative and
positive components by finding the switch point:

$\begin{split}y(\boldsymbol{x})=y^{*}\Rightarrow\mu(\boldsymbol{x})+\sigma(\boldsymbol{x})z=y^{*}\Rightarrow
z_{0}=\frac{y^{*}-\mu(\boldsymbol{x})}{\sigma(\boldsymbol{x})}\end{split}$
(A-7)

choosing $z_{0}$ as our switch point converts the integration to:

$\begin{split}\alpha_{EI}(\boldsymbol{x})=\underbrace{\int_{-\infty}^{z_{0}}(\mu(\boldsymbol{x})+\sigma(\boldsymbol{x})z-y^{*})\phi(z)\mathrm{d}z}_{\text{zero
since
}z<z_{0},I(\boldsymbol{x})=0}+\int_{z_{0}}^{\infty}(\mu(\boldsymbol{x})+\sigma(\boldsymbol{x})z-y^{*})\phi(z)\mathrm{d}z\end{split}$
(A-8)

substituting $\phi(z)=\frac{1}{\sqrt{2\pi}}\exp(-z^{2}/2)$ in Eq. A-8 the
integrating provides the AF:

$\begin{split}\begin{aligned}
\alpha_{EI}(\boldsymbol{x})&=\int_{z_{0}}^{\infty}\left(\mu(\boldsymbol{x})-y^{*}\right)\phi(z)\mathrm{d}z+\int_{z_{0}}^{\infty}\sigma(\boldsymbol{x})z\frac{1}{\sqrt{2\pi}}e^{-z^{2}/2}\mathrm{\leavevmode\nobreak\
d}z\\\
&=\left(\mu(\boldsymbol{x})-y^{*}\right)\underbrace{\int_{z_{0}}^{\infty}\phi(z)\mathrm{d}z}_{1-\Phi\left(z_{0}\right)\equiv
1-\mathrm{CDF}\left(z_{0}\right)}+\frac{\sigma(\boldsymbol{x})}{\sqrt{2\pi}}\int_{z_{0}}^{\infty}ze^{-z^{2}/2}\mathrm{\leavevmode\nobreak\
d}z\\\
&=\left(\mu(\boldsymbol{x})-y^{*}\right)\left(1-\Phi\left(z_{0}\right)\right)-\frac{\sigma(\boldsymbol{x})}{\sqrt{2\pi}}\int_{z_{0}}^{\infty}\left(e^{-z^{2}/2}\right)^{\prime}\mathrm{d}z\\\
&=\left(\mu(\boldsymbol{x})-y^{*}\right)\left(1-\Phi\left(z_{0}\right)\right)-\frac{\sigma(\boldsymbol{x})}{\sqrt{2\pi}}\left[e^{-z^{2}/2}\right]_{z_{0}}^{\infty}\\\
&=\left(\mu(\boldsymbol{x})-y^{*}\right)\underbrace{\left(1-\Phi\left(z_{0}\right)\right)}_{\Phi\left(-z_{0}\right)}+\sigma(\boldsymbol{x})\phi\left(z_{0}\right)\end{aligned}\end{split}$
(A-9)

or:

$\begin{split}\alpha_{EI}(\boldsymbol{x})=(\mu(\boldsymbol{x})-y^{*})\Phi(\frac{\mu(\boldsymbol{x})-y^{*}}{\sigma(\boldsymbol{x})})+\sigma(\boldsymbol{x})\phi(\frac{\mu(\boldsymbol{x})-y^{*}}{\sigma(\boldsymbol{x})})\end{split}$
(A-10)

##  Appendix B Fidelity Kernels of Single-Task Multi-Fidelity GP

Single-task multi-fidelity GP (STGP) uses two fidelity features $(1)$ the data
fidelity parameter, $x_{a}$, which distinguishes between different fidelity
sources, and $(2)$ iteration fidelity parameter, $x_{b}$, which is optional
and usually exists in hyperparameter tuning problem. These two features are
used in $e_{i}(\cdot)$ which are user-defined functions that model the cross-
source correlations in Eq. 20. The formulation of these functions is as
follows:

$\begin{split}e_{1}(x_{a},x^{\prime}_{a})=(1-x_{a})(1-x^{\prime}_{a})(1+x_{a}x^{\prime}_{a})^{p}\end{split}$
(B-1)

where $p$ is the degree of polynomial (which needs to be estimated) and has a
Gamma prior. $e_{3}$ is defined similarly but for the second fidelity:

$\begin{split}e_{3}(x_{b},x^{\prime}_{b})=(1-x_{b})(1-x^{\prime}_{b})(1+x_{b}x^{\prime}_{b})^{p}\end{split}$
(B-2)

Finally, $e_{2}$ is the interaction term with four deterministic terms and one
polynomial kernel:

$\begin{split}e_{2}([x_{a},x_{b}]^{T},[x_{a}^{\prime},x_{b}^{\prime}]^{T})=(1-x_{b})(1-x^{\prime}_{b})(1-x_{a})(1-x^{\prime}_{a})(1+[x_{a},x_{b}]^{T}[x_{a}^{\prime},x_{b}^{\prime}]^{T})^{p}\end{split}$
(B-3)

##  Appendix C Table of Numerical Examples

Table 1 lists the analytic functions studied in section 4.1. The error of each
LF source with respect to the corresponding HF source is calculated via
relative root mean squared error (RRMSE):

$RRMSE=\sqrt{\frac{(\boldsymbol{y}_{l}-\boldsymbol{y}_{h})^{T}(\boldsymbol{y}_{l}-\boldsymbol{y}_{h})}{10000\times
var(\boldsymbol{y}_{h})}}$ (C-1)

where $\boldsymbol{y}_{l}$ and $\boldsymbol{y}_{h}$ are vectors of size
$10000\times 1$ that store random samples taken from the LF and HF sources,
respectively.

Name |  Source ID |  Formulation |  $n$ |  RRMSE | Cost
---|---|---|---|---|---
Double-well Potential |  HF |  $0.6x^{4}-0.3x^{3}-3x^{2}+2x$ |  $5$ |  $-$ | $1000$
2-6 |  LF |  $0.6x^{4}-0.3x^{3}-3x^{2}-1.2x$ |  $0$ |  $1.14$ | $1$
Rosenbrock |  HF |  $(1-x_{1})^{2}+100(x_{2}-x_{1}^{2})^{2}-456.3$ |  5 |  $-$ | $1000$
2-6 |  LF |  $(1-x_{1})^{2}+100$ |  $10$ |  $1.42$ | $1$
Borehole |  HF |  $\frac{2\pi T_{u}(H_{u}-H_{l})}{\ln(\frac{r}{r_{w}})(1+\frac{2LT_{u}}{\ln(\frac{r}{rw})r_{w}^{2}k_{w}}+\frac{T_{u}}{T_{l}})}$ |  $5$ |  $-$ | $1000$
2-6 |  LF1 |  $\frac{2\pi T_{\mathrm{u}}(H_{\mathrm{u}}-0.8H_{l})}{\ln(\frac{r}{r_{w}})(1+\frac{1LT_{\mathrm{u}}}{\ln(\frac{r}{r_{w}})r_{w}^{2}k_{w}}+\frac{T_{u}}{T_{l}})}$ |  $5$ |  $4.40$ | $100$
2-6 |  LF2 |  $\frac{2\pi T_{u}(H_{u}-H_{l})}{\ln(\frac{r}{r_{w}})(1+\frac{8LT_{u}}{\ln(\frac{r}{rw})r_{w}^{2}k_{w}}+0.75\frac{T_{u}}{T_{l}})}$ |  $50$ |  $1.54$ | $10$
2-6 |  LF3 |  $\frac{2\pi T_{\mathrm{u}}(1.09H_{u}-H_{l})}{\ln(\frac{4r}{r_{w}})(1+\frac{3LT_{u}}{\ln(\frac{r}{r_{w}})r_{w}^{2}k_{w}}+\frac{T_{u}}{T_{l}})}$ |  $5$ |  $1.30$ | $100$
2-6 |  LF4 |  $\frac{2\pi T_{\mathrm{u}}(1.05H_{u}-H_{l})}{\ln(\frac{2r}{r_{w}})(1+\frac{3LT_{iu}}{\ln(\frac{r}{\tau w})r_{w}^{2}k_{\mathrm{W}}}+\frac{T_{\mathrm{u}}}{T_{l}})}$ |  $50$ |  $1.3$ | $10$
Wing |  HF |  $0.36s_{w}^{0.758}w_{fw}^{0.0035}(\frac{A}{\cos^{2}(\Lambda)})^{0.6}q^{0.006}\times\lambda^{0.04}(\frac{100t_{c}}{\cos(\Lambda)})^{-0.3}(N_{z}W_{dg})^{0.49}+s_{w}w_{p}$ |  $5$ |  $-$ | $1000$
2-6 |  LF1 |  $0.36s_{w}^{0.758}w_{fw}^{0.0035}(\frac{A}{\cos^{2}(\Lambda)})^{0.6}q^{0.006}\times\lambda^{0.04}(\frac{100t_{c}}{\cos(\Lambda)})^{-0.3}(N_{z}W_{dg})^{0.49}+w_{p}$ |  $5$ |  $0.19$ | $100$
2-6 |  LF2 |  $0.36s_{w}^{0.8}w_{fw}^{0.0035}(\frac{A}{\cos^{2}(\Lambda)})^{0.6}q^{0.006}\times\lambda^{0.04}(\frac{100t_{c}}{\cos(\Lambda)})^{-0.3}(N_{z}W_{dg})^{0.49}+w_{p}$ |  $10$ |  $1.14$ | $10$
2-6 |  LF3 |  $0.36s_{w}^{0.9}w_{fw}^{0.0035}(\frac{A}{\cos^{2}(\Lambda)})^{0.6}q^{0.006}\times\lambda^{0.04}(\frac{100t_{c}}{\cos(\Lambda)})^{-0.3}(N_{z}W_{dg})^{0.49}$ |  $50$ |  $5.75$ | $1$
Table 1: List of analytic functions: The examples have a diverse degree of
dimensionality, number of sources, and complexity. $n$ denotes the number of
initial samples and the relative root mean squared error (RRMSE) of an LF
source is calculated by comparing its output to that of the HF source at
$10000$ random points, see Eq. C-1. For Borehole, LF3 and LF4 become the first
and second LF sources, respectively, once LMGP identifies that the listed LF1
and LF2 in this table are highly biased.

##  Appendix D Effect of Highly Biased Low-fidelity Sources

In Section 4.1, we exclude two LF sources from Borehole due to their high
discrepancy with respect to the HF source, see Figure 6. Below, we summarize
the performance of LMGPCA on this problem without removing these two sources.
As it can be observed in Figure 16(a), the optimization is terminated based on
the second convergence metric which caps the maximum number of iterations
without improvement in the best HF sample (i.e., $y_{l}^{*}$ in Eq. 31).

(a) Convergence history

(b) Frequency of source utilization

Figure 16: Effect of highly biased and inexpensive low-fidelity sources: While
LMGPCA effectively samples from all sources (considering their costs and
contribution to the initial data), it converges to an incorrect solution since
LF1 and LF2 are highly biased. The initial data are not included in 16(b).

##  Appendix E List of Abbreviations and Symbols

Abbreviation | Explanation
---|---
AF | Acquisition Function
BO | Bayesian optimization
BoTorch | Multi-fidelity BO with BoTorch
Co-K | Co-Kriging
EI | Expected Improvement
GP | Gaussian Process
HF | High-fidelity
HOIP | Hybrid Organic–Inorganic Perovskite
KG | Knowledge Gradient
LF | Low-Fidelity
LMGP | Latent Map Gaussian Process
LMGPCA | Proposed MFCA BO approach
LMGPEI | Single-fidelity BO whose emulator and AF are LMGP and EI, respectively
LMGPPI | Single-fidelity BO whose emulator and AF are LMGP and PI, respectively
MF | Multi-Fidelity
MFCA | Multi-Fidelity Cost-Aware
MLE | Maximum Likelihood Estimation
MLEI | Most Likely Expected Improvment
MSE | Mean Squared Error
NTA | Nanolaminate Ternary Alloy
PI | Probability of Improvement
RRMSE | Relative Root Mean Squared Error
SF | Single-Fidelity
STGP | Single-Task Multi-Fidelity Gaussian
Symbol | Description
---|---
$\boldsymbol{A}$ | Rectangular matrix that maps $\boldsymbol{\zeta(t)}$ to $\boldsymbol{z(t)}$
$c(\boldsymbol{x},\boldsymbol{x}^{\prime})$ | Covariance function
$ds$ | Number of data sources
$dt$ | Dimension of categorical inputs
$dz$ | Dimension of the latent map
$dx$ | Dimension of numerical inputs
$\boldsymbol{h}(s)=[h_{1},\cdots,h_{dh}]^{T}$ | Latent representation of data source $s$
$I(\boldsymbol{x})$ | Utility function
${}^{\prime}\boldsymbol{j}^{\prime}$ | Categorical vector of size ${n_{j}\times 1}$ whose elements are all set to ${}^{\prime}j^{\prime}$
$\mathcal{N}\left(\mu(\boldsymbol{x}),\sigma^{2}(\boldsymbol{x})\right)$ | Normal distribution with mean $\mu(\boldsymbol{x})$ and standard deviation $\sigma^{2}(\boldsymbol{x})$
$n_{j}$ | Number of samples obtained from $s(j)$ (i.e., source $j$)
$l_{i}$ | Number of distinct levels in ${i}^{th}$ categorical input
$\boldsymbol{R}$ | Correlation matrix
$r(.,.)$ | Parametric correlation function
$s=\\{^{\prime}1^{\prime},\cdots,^{\prime}ds^{\prime}\\}$ | Categorical variable whose $j^{th}$ element corresponds to data source $j$
$\boldsymbol{t}$ | Categorical inputs (all except source indicator)
$\boldsymbol{U}_{j}$ |  $n_{j}\times(dx+dt)$ feature matrix obtained from $s(j)$
$\boldsymbol{u}$ | Mixed inputs
$\boldsymbol{x}$ | Input vector
$y(\boldsymbol{x})$ | Output/response
$\boldsymbol{y}_{j}$ |  ${n_{j}\times 1}$ vector of responses obtained from $s(j)$
$\boldsymbol{z(t)}$ | Points on latent map corresponding to combination $\boldsymbol{t}$ of the categorical variables
$\alpha(\boldsymbol{x})$ | Acquisition Function
$\boldsymbol{\zeta}(\boldsymbol{t})$ | Unique prior vector representation of $\boldsymbol{t}$
$\xi(\boldsymbol{x})$ | Zero-mean GP
$\sigma^{2}$ | Variance of process
$\Phi(z)$ | Cumulative density function (CDF)
$\phi(z)$ | Probability density function (PDF)
$\boldsymbol{\Omega}$ | $\operatorname{diag}(\boldsymbol{\omega})$
$\boldsymbol{\omega}$ | Scale parameters
$\otimes$ | Kronecker product

## References

  * [1] Ryan-Rhys Griffiths and José Miguel Hernández-Lobato  “Constrained Bayesian optimization for automatic chemical design using variational autoencoders”  In _Chemical science_ 11.2 Royal Society of Chemistry, 2020, pp. 577–586
  * [2] Y. Zhang, D.. Apley and W. Chen  “Bayesian Optimization for Materials Design with Mixed Quantitative and Qualitative Variables”  In _Sci Rep_ 10.1, 2020, pp. 4924  DOI: 10.1038/s41598-020-60652-9
  * [3] Yiqun Wang, Akshay Iyer, Wei Chen and James M. Rondinelli  “Featureless adaptive optimization accelerates functional electronic materials design”  In _Applied Physics Reviews_ 7.4, 2020, pp. 041403
  * [4] Anh Tran, Julien Tranchida, Tim Wildey and Aidan P. Thompson  “Multi-fidelity machine-learning with uncertainty quantification and Bayesian optimization for materials design: Application to ternary random alloys”  In _The Journal of Chemical Physics_ 153.7, 2020, pp. 074705  DOI: 10.1063/5.0015672
  * [5] Elvis Osamudiamhen Ebikade et al.  “Active learning-driven quantitative synthesis–structure–property relations for improving performance and revealing active sites of nitrogen-doped carbon for the hydrogen evolution reaction”  In _Reaction Chemistry & Engineering_ 5.12 Royal Society of Chemistry, 2020, pp. 2134–2147
  * [6] Jasper Snoek, Hugo Larochelle and Ryan P Adams  “Practical bayesian optimization of machine learning algorithms”  In _Advances in neural information processing systems_ 25, 2012
  * [7] Benjamin Burger et al.  “A mobile robotic chemist”  In _Nature_ 583.7815 Nature Publishing Group, 2020, pp. 237–241
  * [8] Takamitsu Matsubara et al.  “Data-efficient human training of a care motion controller for human transfer assistant robots using bayesian optimization”  In _2016 6th IEEE International Conference on Biomedical Robotics and Biomechatronics (BioRob)_ , 2016, pp. 606–611  IEEE
  * [9] Carl Edward Rasmussen  “Gaussian processes for machine learning”, 2006
  * [10] Matthew Plumlee and Daniel W. Apley  “Lifted Brownian Kriging Models”  In _Technometrics_ 59.2, 2017, pp. 165–177  DOI: 10.1080/00401706.2016.1211555
  * [11] I. Hassaninia, R. Bostanabad, W. Chen and H. Mohseni  “Characterization of the Optical Properties of Turbid Media by Supervised Learning of Scattering Patterns”  In _Sci Rep_ 7.1, 2017, pp. 15259  DOI: 10.1038/s41598-017-15601-4
  * [12] Siyu Tao et al.  “Enhanced Gaussian Process Metamodeling and Collaborative Optimization for Vehicle Suspension Design Optimization”  In _ASME 2017 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference_ 2B American Society of Mechanical Engineers
  * [13] Yichi Zhang, Siyu Tao, Wei Chen and Daniel W. Apley  “A Latent Variable Approach to Gaussian Process Modeling with Qualitative and Quantitative Factors”  In _Technometrics_ 62.3, 2019, pp. 291–302  DOI: 10.1080/00401706.2019.1638834
  * [14] Olivier Roustant et al.  “Group kernels for Gaussian process metamodels with categorical inputs”  In _SIAM/ASA Journal on Uncertainty Quantification_ 8.2, 2020, pp. 775–806
  * [15] Qiong Zhang et al.  “Mixed-input Gaussian process emulators for computer experiments with a large number of categorical levels”  In _Journal of Quality Technology_ , 2020, pp. 1–11  DOI: 10.1080/00224065.2020.1778431
  * [16] Rohit Tripathy, Ilias Bilionis and Marcial Gonzalez  “Gaussian processes with built-in dimensionality reduction: Applications to high-dimensional uncertainty propagation”  In _Journal of Computational Physics_ 321, 2016, pp. 191–223
  * [17] Songqing Shan and G. Wang  “Metamodeling for high dimensional simulation-based design problems”  In _Journal of Mechanical Design_ 132.5, 2010
  * [18] Liwei Wang et al.  “Scalable Gaussian Processes for Data-Driven Design Using Big Data With Categorical Factors”  In _Journal of Mechanical Design_ 144.2, 2021  DOI: 10.1115/1.4052221
  * [19] , 2021, pp. 3133–3141  URL: http://proceedings.mlr.press/v130/stanton21a.html
  * [20] Jacob R. Gardner et al.  “Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration”  In _arXiv preprint arXiv:1809.11165_ , 2018
  * [21] James Hensman, Nicolo Fusi and Neil D Lawrence  “Gaussian processes for big data”  In _arXiv preprint arXiv:1309.6835_ , 2013
  * [22] Ramin Bostanabad  “Quantification of Microstructure Induced Uncertainty in Multiscale Materials with Random Processes”, 2019
  * [23] Shan Ba and V. Joseph  “Composite Gaussian process models for emulating expensive functions”  In _The Annals of Applied Statistics_ 6.4, 2012, pp. 1838–1860  DOI: 10.1214/12-aoas570
  * [24] Maximilian Balandat et al.  “BoTorch: A framework for efficient Monte-Carlo Bayesian optimization”  In _Advances in Neural Information Processing Systems_ 33, 2020
  * [25] James T. Wilson, Riccardo Moriconi, Frank Hutter and Marc Peter Deisenroth  “The reparameterization trick for acquisition functions”  In _arXiv preprint arXiv:1712.00424_ , 2017
  * [26] Nicholas Oune and Ramin Bostanabad  “Latent map Gaussian processes for mixed variable metamodeling”  In _Computer Methods in Applied Mechanics and Engineering_ 387 Elsevier, 2021, pp. 114128
  * [27] Ramin Bostanabad et al.  “Leveraging the nugget parameter for efficient Gaussian process modeling”  In _International journal for numerical methods in engineering_ 114.5 Wiley Online Library, 2018, pp. 501–516
  * [28] R. Bostanabad et al.  “Leveraging the nugget parameter for efficient Gaussian process modeling”  In _International Journal for Numerical Methods in Engineering_ 114.5, 2018, pp. 501–516  DOI: 10.1002/nme.5751
  * [29] Jianjun Wang, Yizhong Ma, Linhan Ouyang and Yiliu Tu  “A new Bayesian approach to multi-response surface optimization integrating loss function with posterior probability”  In _European Journal of Operational Research_ 249.1 Elsevier, 2016, pp. 231–237
  * [30] Raul Astudillo and Peter Frazier  “Bayesian optimization of composite functions”  In _International Conference on Machine Learning_ , 2019, pp. 354–363  PMLR
  * [31] Michael A Gelbart, Jasper Snoek and Ryan P Adams  “Bayesian optimization with unknown constraints”  In _arXiv preprint arXiv:1403.5607_ , 2014
  * [32] Jialei Wang, Scott C Clark, Eric Liu and Peter I Frazier  “Parallel Bayesian global optimization of expensive functions”  In _Operations Research_ 68.6 INFORMS, 2020, pp. 1850–1865
  * [33] Yunxiang Zhang, Xiangyu Zhang and Peter I Frazier  “Two-step Lookahead Bayesian Optimization with Inequality Constraints”  In _arXiv preprint arXiv:2112.02833_ , 2021
  * [34] Hao Wang, Bas Stein, Michael Emmerich and Thomas Back  “A new acquisition function for Bayesian optimization based on the moment-generating function”  In _2017 IEEE International Conference on Systems, Man, and Cybernetics (SMC)_ , 2017, pp. 507–512  IEEE
  * [35] Donald R Jones  “A taxonomy of global optimization methods based on response surfaces”  In _Journal of global optimization_ 21.4 Springer, 2001, pp. 345–383
  * [36] H-M Gutmann  “A radial basis function method for global optimization”  In _Journal of global optimization_ 19.3 Springer, 2001, pp. 201–227
  * [37] Matthias Schonlau, William J Welch and Donald R Jones  “Global versus local search in constrained optimization of computer models”  In _Lecture notes-monograph series_ JSTOR, 1998, pp. 11–25
  * [38] Wenlong Lyu et al.  “An efficient bayesian optimization approach for automated optimization of analog circuits”  In _IEEE Transactions on Circuits and Systems I: Regular Papers_ 65.6 IEEE, 2017, pp. 1954–1967
  * [39] Peter I Frazier, Warren B Powell and Savas Dayanik  “A knowledge-gradient policy for sequential information collection”  In _SIAM Journal on Control and Optimization_ 47.5 SIAM, 2008, pp. 2410–2439
  * [40] Maximilian Balandat et al.  “BoTorch: a framework for efficient Monte-Carlo Bayesian optimization”  In _Advances in neural information processing systems_ 33, 2020, pp. 21524–21538
  * [41] Guanghui Zhu and Ruancheng Zhu  “Accelerating Hyperparameter Optimization of Deep Neural Network via Progressive Multi-Fidelity Evaluation”  In _Pacific-Asia Conference on Knowledge Discovery and Data Mining_ , 2020, pp. 752–763  Springer
  * [42] Marius Lindauer et al.  “BOAH: A tool suite for multi-fidelity bayesian optimization & analysis of hyperparameters”  In _arXiv preprint arXiv:1908.06756_ , 2019
  * [43] Hyunghun Cho et al.  “Basic enhancement strategies when using Bayesian optimization for hyperparameter tuning of deep neural networks”  In _IEEE Access_ 8 IEEE, 2020, pp. 52588–52608
  * [44] Jian Wu, Saul Toscano-Palmerin, Peter I Frazier and Andrew Gordon Wilson  “Practical multi-fidelity Bayesian optimization for hyperparameter tuning”  In _Uncertainty in Artificial Intelligence_ , 2020, pp. 788–798  PMLR
  * [45] Tharathep Phiboon et al.  “Experiment and computation multi-fidelity multi-objective airfoil design optimization of fixed-wing UAV”  In _Journal of Mechanical Science and Technology_ 35.9 Springer, 2021, pp. 4065–4072
  * [46] Qi Sun et al.  “Correlated multi-objective multi-fidelity optimization for HLS directives design”  In _ACM Transactions on Design Automation of Electronic Systems (TODAES)_ 27.4 ACM New York, NY, 2022, pp. 1–27
  * [47] Syrine Belakaria, Aryan Deshwal and Janardhan Rao Doppa  “Multi-fidelity multi-objective Bayesian optimization: An output space entropy search approach”  In _Proceedings of the AAAI Conference on artificial intelligence_ 34.06, 2020, pp. 10035–10043
  * [48] Johan Dahlin, Thomas Bo Schön and Mattias Villani  “Approximate inference in state space models with intractable likelihoods using Gaussian process optimisation”, 2014
  * [49] Mahdi Imani, Seyede Fatemeh Ghoreishi, Douglas Allaire and Ulisses M Braga-Neto  “MFBO-SSM: Multi-fidelity Bayesian optimization for fast inference in state-space models”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 33.01, 2019, pp. 7858–7865
  * [50] Edwin V Bonilla, Kian Chai and Christopher Williams  “Multi-task Gaussian process prediction”  In _Advances in neural information processing systems_ 20, 2007
  * [51] David A Stainforth et al.  “Uncertainty in predictions of the climate response to rising levels of greenhouse gases”  In _Nature_ 433.7024 Nature Publishing Group, 2005, pp. 403–406
  * [52] Weizhao Zhang et al.  “A numerical Bayesian-calibrated characterization method for multiscale prepreg preforming simulations with tension-shear coupling”  In _Composites Science and Technology_ 170 Elsevier, 2019, pp. 15–24
  * [53] Siyu Tao et al.  “Input Mapping for Model Calibration with Application to Wing Aerodynamics”  In _AIAA Journal_ 57.7, 2019, pp. 2734–2745  DOI: 10.2514/1.J057711
  * [54] S Ashwin Renganathan, Vishwas Rao and Ionel M Navon  “CAMERA: A Method for Cost-aware, Adaptive, Multifidelity, Efficient Reliability Analysis”  In _arXiv preprint arXiv:2203.01436_ , 2022
  * [55] Bobak Shahriari et al.  “Taking the human out of the loop: A review of Bayesian optimization”  In _Proceedings of the IEEE_ 104.1 IEEE, 2015, pp. 148–175
  * [56] Shion Takeno et al.  “Multi-fidelity Bayesian optimization with max-value entropy search and its parallelization”  In _International Conference on Machine Learning_ PMLR, pp. 9334–9345
  * [57] Niranjan Srinivas, Andreas Krause, Sham M Kakade and Matthias W Seeger  “Information-theoretic regret bounds for gaussian process optimization in the bandit setting”  In _IEEE transactions on information theory_ 58.5 IEEE, 2012, pp. 3250–3265
  * [58] Daniel Russo and Benjamin Van Roy  “Learning to optimize via posterior sampling”  In _Mathematics of Operations Research_ 39.4 INFORMS, 2014, pp. 1221–1243
  * [59] José Miguel Hernández-Lobato, Matthew W Hoffman and Zoubin Ghahramani  “Predictive entropy search for efficient global optimization of black-box functions”  In _Advances in neural information processing systems_ 27, 2014
  * [60] Rémi Pautrat, Konstantinos Chatzilygeroudis and Jean-Baptiste Mouret  “Bayesian optimization with automatic prior selection for data-efficient direct policy search”  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 7571–7578  IEEE
  * [61] Wolfgang Ponweiser, Tobias Wagner, Dirk Biermann and Markus Vincze  “Multiobjective optimization on a limited budget of evaluations using model-assisted S-Metric selection”  In _International Conference on Parallel Problem Solving from Nature_ , 2008, pp. 784–794  Springer
  * [62] Yichi Zhang, Daniel W Apley and Wei Chen  “Bayesian optimization for materials design with mixed quantitative and qualitative variables”  In _Scientific reports_ 10.1 Nature Publishing Group, 2020, pp. 1–13
  * [63] Jia Wu et al.  “Hyperparameter optimization for machine learning models based on Bayesian optimization”  In _Journal of Electronic Science and Technology_ 17.1 Elsevier, 2019, pp. 26–40
  * [64] Chan-Uk Yeom and Keun-Chang Kwak  “Performance evaluation of automobile fuel consumption using a fuzzy-based granular model with coverage and specificity”  In _Symmetry_ 11.12 MDPI, 2019, pp. 1480