{'fraction_non_alphanumeric': 0.03595163744094804, 'fraction_numerical': 0.010681399631675874, 'mean_word_length': 4.501850220264317, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-Fidelity Cost-Aware Bayesian Optimization

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.


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 lowand 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  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 (x1 , x2) and one categorical input (t) that has two levels (A and Z). In each scenario, x2 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.
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 tobe-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.

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.

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: where x = [x 1 , x 2 , . . . , x dx ] T is the input vector, y(x) is the output/response, m(x) = [m 1 (x), . . . , m dβ (x)] are a set of parametric basis functions, β = [β 1 , β 2 , ..., β dβ ] T are unknown coefficients, and ξ(x) is a zero-mean GP whose covariance function or kernel is: where σ 2 is the variance of the process and r(., .) is a parametric correlation function. A common choice for r(., .) is: where ω = [ω 1 , . . . , ω dx ] T , −∞ < ω i < ∞ are the scale parameters and Ω = diag(ω). 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 t = [t 1 , . . . , t dt ] T where variable t i has l i distinct levels. For example, t 1 = {M ale, F emale} and t 2 = {P ersian, American, Spanish} are two categorical inputs where l 1 = 2 and l 2 = 3. To handle mixed inputs, 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 space 1 . This mapping function can be incorporated into any standard correlation function such as the Gaussian which is reformulated as follows for mixed inputs: or equivalently, where · 2 denotes the Euclidean 2-norm and z(t) = [z 1 (t), . . . , z dz (t)] 1×dz is the location in the learned manifold corresponding to the specific combination of the categorical variables denoted by 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 transformation 2 that maps these unique vectors into a compact manifold (aka latent space) with dimensionality dz: where t denotes a specific combination of the categorical variables, z(t) is the 1 × dz posterior latent representation of t, ζ(t) is the unique prior vector representation of t, and A is a rectangular matrix that maps ζ(t) to z(t). There are various methods for constructing the prior vectors ζ and we refer the reader to [26] for more details. In this paper, we use grouped one-hot encoding which makes ζ(t) and A to be of sizes 1 × dt i=1 l i and dt i=1 l i × dz, respectively. For instance, in the above example the grouped one-hot encoded version of the combination t = [F emale, American] T is ζ(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 (β, A, ω, and σ 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(x), that is: or equivalently: where log(·) is the natural logarithm, | · | denotes the determinant operator, y = [y 1 , . . . , y n ] T is the n × 1 vector of outputs in the training data, R is the n × n correlation matrix with the (i, j) th element R ij = r(x i , x j ) for i, j = 1, . . . , n , and M is the n × dβ matrix with the (i, j) th element M ij = m j (x i ) for i = 1, . . . , n and j = 1, . . . , dβ.
Once the hyperparameters are estimated, the conditional distribution formulas are used to predict the response distribution at the arbitrary point x * . The mean and variance of this normal distribution are: where E denotes expectation, m(x * ) = [m 1 (x * ), . . . , m dβ (x * )], r(x * ) is an n × 1 vector with the i th element r(x i , x * ), and g(x * ) = m(x * ) − M T R −1 r(x * ).

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., E[I(x)] where I(x) is the user-defined utility function. These queries are used in the auxiliary optimization problem 3 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., E[I(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. 4. Query f (x) at x k+1 to obtain y k+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 myopic 4 , KG is look-ahead. Any AF calculates the expected value of a user-defined utility function conditioned on the available data D, that is: 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 (x), this AF uses the following utility function: where y(x) is the emulator-based prediction at x. Based on Eq. 13, if y(x) is less than y * , the point x has zero utility. Assuming a GP is used for emulation, y(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 P I (x) is [35,36]: where µ(x) and σ(x) are defined in Eq. 10 and Eq. 11, and Φ(z) is the cumulative density function (CDF) of the standard normal random variable z. Eq. 14 clearly indicates that α P I (x) favors exploitation because Φ(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(x), I P I (x) is either 0 or 1), EI rewards large improvements over y * by adopting the following utility function: 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): where µ(x) and σ(x) are given in Eq. 10 and Eq. 11, respectively, and φ(z) is the probability density function (PDF) of z. Eq. 16 clearly demonstrates that α EI (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 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 x at iteration k + 1 as: where max µ k (x) denotes the maximum mean prediction of the GP trained on D k . The expectation operation in Eq. 17 appears due to the fact that y k+1 is not observed yet and α KG (x) is relying on the predictive distribution provided by the GP that is trained on D k . This expectation cannot be calculated analytically and hence a Monte Carlo estimate is used in practice: where max µ k+1 m (x) is calculated by first drawing a sample at x from the GP that is trained on D k and then retraining the GP on D k ∪ (x, y m ) where y m is response of the drawn sample. In practice, a small value must be chosen for M since maximizing α KG (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.

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 unbalanced 5 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): where ⊗ denotes the Kronecker product and Σ is a symmetric positive-definite matrix of size 3 × 3. This reformulation assumes that all the responses (regardless of the source) follow a multi-variate normal distribution and that the matrix Σ ij captures the overall correlation between sources i and j 6 . While this method can fuse any number of data sources, it fails to accurately capture cross-source correlations since the matrix Σ 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 [20,24]. 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: where c i (x, x ) are Matern kernels 7 that characterize the spatial correlations across the numerical inputs and e i (·) are user-defined functions that model the cross-source correlations. e 1 (x a , x a ) and e 3 (x b , x b ) are bias kernels that aim to take the discrepancies among the sources into account while 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 (·) 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 [24,25,55,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 GPpredictive 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, traceaware 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].

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.

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 = { 1 , · · · , ds } whose j th element corresponds to data source j for j = 1, · · · , ds. After this augmentation, the inputs and outputs of all the datasets are concatenated as (following the notation of Figure 2): 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), U j and y j are, respectively, the n j × (dx + dt) feature matrix and the n j × 1 vector of responses obtained from s(j), and j is a categorical vector of size n j × 1 whose elements are all set to j . Once the {U , y} dataset is built, it is directly fed into LMGP to build an MF emulator. 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., 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 Σ 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.
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 /0*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.
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 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): where h(s) = [h 1 , · · · , 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 A h : In our studies, we always design the prior by one-hot-encoding the categorical variable s = { 1 , · · · , ds } that identifies the data source and estimate all of LMGP's hyperparameters via MLE: where the elements of R are now built using Eq. 22.
We now use Eq. 22 to explain the relation between the latent fidelity representations, i.e., 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 is: So sources with similar fidelities which provide highly correlated responses, must have similar latent representations, i.e., h i (s) h i (s ). This relation is illustrated in Figure 6 where sources with similar fidelities are encoded by close-by points in the manifold.

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 4a 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 4b for a one-dimensional example. 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 samples 8 during optimization.
To develop the AF for the j th LF source with j = [1, · · · , ds] and j = l where l denotes the HF source, we follow Section 2.2.1 and define the improvement (for a maximization problem) as y j (x)-y * j where y j (x) denotes the LMGP-based prediction at 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: where µ j (x) and σ j (x) are defined in Eq. 10 and Eq. 11, respectively. In Eq. 26 the (µ j (x) − y * j ) and σ j (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: which is used for the j th LF source in our framework. We obtain the corresponding AF by substituting I LF (x; j) in Eq. 12: The integral is zero for y j (x) < y * j so we find the corresponding switch point in terms of z: Inserting Eq. 29 in Eq. 28 yields: 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: We use α HF (x; l) and α LF (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: and where O(j) is the cost associated with taking a single sample from source j and x k+1 is the point that source j k+1 ∈ {1, · · · , 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 sourcewise 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 optimality 9 ). 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, φ(z) and Φ(z) have comparable values up to E[z] = 0 but the ratio Φ(z)/φ(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 σ j (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 Φ(z)/φ(z) can be quite large. Regarding the σ j (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.

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.

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): • LMGP CA : Our proposed MFCA BO approach.
• LMGP EI : Single-fidelity BO whose emulator and AF are LMGP and EI, respectively.
• LMGP PI : 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 LMGP CA 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 3. Solve the auxiliary optimization problem: 4. Query f (x; j) at point x k+1 from source j k+1 to obtain y k+1 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 LMGP CA , LMGP EI , and LMGP PI 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 LMGP EI and LMGP PI 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.

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). 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 LMGP CA 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. 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. 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 LMGP CA consistently outperforms all other methods across the four examples. In particular, LMGP CA demonstrates the advantage of leveraging inexpensive LF sources in BO by accelerating the convergence without sacrificing the accuracy (compare LMGP CA to LMGP EI and LMGP PI for any of the examples in Figure 7). Additionally, unlike BoTorch, the performance of LMGP CA is robust to the input dimensionality and sampling costs. For instance, BoTorch estimates the optimum as y * l = 10.02 in Figure 7c 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 7d. These issues are resolved in LMGP CA by the intricate interplay between Φ(z) and φ(z) as explained in Section 3.2.
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 &RVW î true optimum (compare the blue dots to the horizontal dashed line), LMGP CA outperforms all other methods and is followed by LMGP EI and then LMGP PI (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.
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 LMGP CA require more iterations as they aim to leverage cheap LF sources to reduce the overall costs. LMGP CA needs fewer iterations for convergence than BoTorch except in Figure 9b. This behavior is a result of the termination criteria: the ground truth of Rosenbrock is −456.3 while LMGP CA 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.
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 LMGP CA and  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.
BoTorch in Figure 10 which demonstrates that LMGP CA 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, LMGP CA 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 LMGP CA 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 LMGP CA 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 10d where, unlike BoTorch, LMGP CA 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 LMGP CA in balancing exploration and exploitation. To this end, source-wise sampling orders made by LMGP CA and BoTorch are visualized for a randomly selected repetition in each example, see Figure 11. As it can be observed, LMGP CA alternates between all the sources: 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).
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, LMGP CA effectively uses all sources during the optimization.

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 LMGP CA 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, respectively 10 . The three datasets have similar 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. 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 LMGP EI and LMGP PI 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.

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 12a, 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 12b 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 12c and 12d. As it can be observed, while the manifolds do not match exactly, the relative distances between the latent points are similar.  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.
(i.e., y * l in Eq. 31) as a function of the accumulated costs. As expected, LMGP CA outperforms the two SF methods in NTA but not in HOIP. In Figure 13a we observe that LMGP EI and LMGP PI 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, LMGP CA finds a sub-optimal compound in HOIP since the highly biased LF sources steer the search process in the wrong direction.
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, LMGP CA outperforms both LMGP EI and LMGP PI 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, LMGP CA provides lower accuracy than the SF methods since it is using highly biased LF sources. Even  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.
though LMGP CA is more robust to variations in the initial data, the lower accuracy does not justify its use for HOIP.
We now investigate the resource allocation behavior of LMGP CA . As shown in Figures 15a and 15b, LMGP CA takes equal or fewer iterations to converge (note that since one sample is taken per iteration, this means that LMGP CA 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 15c. 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 LMGP CA can effectively allocate resources based on the initial dataset sizes and data collection costs (see Figures 15d and 15e), highly biased LF sources can steer the search in the wrong direction and, in turn, result in convergence to an incorrect solution.

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.

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.

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: Eq. 13 demonstrates that I(x) is zero for y(x) < y * , so: Hence, to calculate α P I we only need to find Pr(I(x) > 0) which is a function of the random variable y(x). In a GP, the response y(x) follows a normal distribution with mean µ(x) and variance σ 2 (x): We now apply the reparameterization trick to y(x) to calculate PI. Considering z ∼ N (0, 1), then y(x) = µ(x) + σ(x)z is a normal distribution with mean µ(x) and variance σ 2 (x). Then: defining z 0 = y * −µ(x) σ(x) which follows N 0, σ 2 simplifies the above as: where Φ 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(x) = µ(x) + σ(x)z to rewrite Eq. 15 as I(x) = µ(x) + σ(x)z − y * where z is the standard normal random variable. We now insert this utility function into Eq. 12: To eliminate the max operator and simplify the integration, we divide µ(x) + σ(x)z − y * into two negative and positive components by finding the switch point: choosing z 0 as our switch point converts the integration to: substituting φ(z) = 1 √ 2π exp(−z 2 /2) in Eq. A-8 the integrating provides the AF: or:

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 (·) which are userdefined functions that model the cross-source correlations in Eq. 20. The formulation of these functions is as follows: 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: Finally, e 2 is the interaction term with four deterministic terms and one polynomial kernel: 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): RRM SE = (y l − y h ) T (y l − y h ) 10000 × var(y h ) (C-1) where y l and y h are vectors of size 10000 × 1 that store random samples taken from the LF and HF sources, respectively.

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 LMGP CA on this problem without removing these two sources. As it can be observed in Figure 16a, 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).  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 16b. 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   Dimension of numerical inputs h(s) = [h 1 , · · · , h dh ] T Latent representation of data source s I(x) Utility function j Categorical vector of size n j × 1 whose elements are all set to j N µ(x), σ 2 (x) Normal distribution with mean µ(x) and standard deviation σ 2 (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 R Correlation matrix r(., .) Parametric correlation function s = { 1 , · · · , ds } Categorical variable whose j th element corresponds to data source j t Categorical inputs (all except source indicator) U j n j × (dx + dt) feature matrix obtained from s(j) u Mixed inputs x Input vector y(x) Output/response y j n j × 1 vector of responses obtained from s(j) z(t) Points on latent map corresponding to combination t of the categorical variables α(x) Acquisition Function ζ(t) Unique prior vector representation of t ξ(x) Zero-mean GP σ 2 Variance of process Φ(z) Cumulative density function (CDF) φ(z) Probability density function (PDF) Ω diag(ω) ω Scale parameters ⊗ Kronecker product