{'fraction_non_alphanumeric': 0.037073024889752594, 'fraction_numerical': 0.011772180282532328, 'mean_word_length': 4.42328739359546, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transfer Learning with Gaussian Processes for Bayesian Optimization

Bayesian optimization is a powerful paradigm to optimize black-box functions based on scarce and noisy data. Its data efficiency can be further improved by transfer learning from related tasks. While recent transfer models meta-learn a prior based on large amount of data, in the low-data regime methods that exploit the closed-form posterior of Gaussian processes (GPs) have an advantage. In this setting, several analytically tractable transfer-model posteriors have been proposed, but the relative advantages of these methods are not well understood. In this paper, we provide a unified view on hierarchical GP models for transfer learning, which allows us to analyze the relationship between methods. As part of the analysis, we develop a novel closed-form boosted GP transfer model that fits between existing approaches in terms of complexity. We evaluate the performance of the different approaches in large-scale experiments and highlight strengths and weaknesses of the different transfer-learning methods.


INTRODUCTION
Bayesian optimization (BO) is an elegant and powerful approach to black-box optimization. It has been successfully applied to several challenging black-box optimization problems, ranging from hyperparameter optimization (Snoek et al., 2012) to materials design (Zhang et al., 2020) and controller tuning (Calandra et al., 2016). One major advantage of BO is sample efficiency: it is specifically tailored to expensive black- box functions, where each function evaluation is costly, e.g., when laborious experiments on physical systems are involved. For such applications, the benefit of fewer function evaluations outweighs the increase in computational cost to make informed decisions about where to measure next. A key ingredient to BO's sample efficiency is a probabilistic surrogate model of the objective. In the low-data regime, the most commonly employed model is a Gaussian process (GP) (Rasmussen and Williams, 2006), as it provides closed-form posteriors with accurate uncertainty estimates to guide the search for the global optimum.
Many black-box optimization problems are not oneoff tasks, but rather several related instances of the tasks can be encountered. Here, the data efficiency of optimization can be further improved by transferring knowledge from related tasks. Especially in the lowdata regime that BO is specialized on, transfer learning provides great value. To this end, Swersky et al. (2013);Joy et al. (2016); Shilton et al. (2017); Golovin et al. (2017) propose several approaches to transfer learning for BO. These transfer learning models combine related task data and the data of the current task into a joint model that guides BO in the search for the current task's global optimum. This setting implies a certain asymmetry: while we may use data from related tasks, we are only interested in informative models for the current task. Further, the related task data are considered as given and additional samples cannot be acquired. One common approach to account for this asymmetry, often referred to as hierarchical model, is to model the difference between the current and related tasks. However, when combining different task data into a single GP model, the computational complexity scales cubically with the number of data points. This complexity can easily become prohibitive if data from multiple related tasks are available. Several other transfer models with tractable posteriors and more favourable scaling have been proposed, e.g., an ensemble of GPs (Feurer et al., 2018) or a hierarchical model for the mean prior (Golovin et al., 2017). However, their relation to the joint models with full cu-bic complexity, their relative strengths and weaknesses, and underlying assumptions are not well understood.

RELATED WORK AND CONTRIBUTIONS
The problem of speeding up Bayesian optimization by using additional information besides the objective function evaluations is a focus topic in the black-box optimization community. A large body of literature promotes the idea to meta-learn models (Finn et al., 2017;Perrone et al., 2018a;Flennerhag et al., 2019), whole optimizers (Chen et al., 2017) or acquisition functions (Volpp et al., 2020) for BO. These methods are powerful, but require large amounts of training data to extract the characteristics of a class of tasks.
In the low-data regime, methods that rely on the uncertainty estimates of a GP posterior to guide the search outperform large-scale meta-learning. Several GP-based transfer learning approaches have been proposed, which can mostly be summarized into two categories that either (i) build a global GP model, or ( In multi-fidelity BO, which only differs from the transfer learning setting by removing the constraint that related task data cannot be acquired during optimization, Forrester et al. (2007); Marco et al. (2017) successfully employ GP models. All these approaches ultimately compute a GP model on a data set of the size of the joint related and target task data sets and suffer from the cubic complexity. On the other hand, Golovin et al. (2017), Feurer et al. (2018 and Wistuba et al. (2016) maintain separate GPs for each data set to avoid the accumulated scaling. In the hierarchical model presented in Golovin et al. (2017), this reduction in computational complexity stems from neglecting the uncertainty of the source models for the transfer, which is detrimental to the optimization process. Instead, RGPE (Feurer et al., 2018) builds a ranking-based ensemble of the separate GPs, which involves careful tuning of the weight estimator and also does not result in a Bayesian treatment of the uncertainty. In light of the described properties of these two classes of GP-based transfer learning models, the need for models that can bridge the gap between fully con-sidering source uncertainty and not considering it at all becomes evident. Such an intermediate regime may be particularly relevant for optimizing monetarily expensive processes with relatively little duration such as destructive measurements. In the engineering community, efficient hierarchical models have been developed for the special case where related tasks are evaluated on decreasing subsets of inputs, see e.g. Kennedy and O'Hagan (2000); Le Gratiet and Garnier (2014). However, this setting is not suitable for transfer learning from generic historical data.
While computationally efficient approximations to GPs exist that could be applied to any GP model (including those proposed in this paper), e.g. by Lázaro-Gredilla et al. (2010), this is an orthogonal research direction that is agnostic to the structure of the underlying data and treats every data point the same. Contextual optimization (Krause and Ong, 2011) is another related line of research in which context variables influence the process to-be-optimized. By contrast, transfer learning leverages discrete historical tasks as context.
Contributions In this paper, we provide a unifying view on GP based transfer learning models and shed light on their relation and assumptions. This unified view enables a thorough analysis of the methods: their computational complexity, their modelling assumptions and their respective advantages for different transfer scenarios. In addition, we fill the gaps in this unified framework. Firstly, we present a modular version of Hierarchical GP that we name Sequential Hierarchical GP with improved scaling while maintaining competitive performance. Secondly, we develop a new transfer learning model that is a middle ground between existing approaches in terms of computational complexity and optimization performance. This method is conceptually inspired by boosting architectures, and yet results in an analytically tractable posterior, which we name Boosted Hierarchical GP.
A comprehensive empirical evaluation of all methods supports our analytical insights and provides valuable guidance on the different methods' comparative advantages.

PROBLEM STATEMENT
Our goal is to find the optimum of a target black-box function f t : D → R, based on noisy observations D t = {x n , y n } Nt n=1 where each observation y n = f t (x n ) + ε n is corrupted by i.i.d. zero-mean Gaussian noise, ε n ∼ N (0, σ 2 t ). We model our prior belief over the target function f t with a GP (Rasmussen and Williams, 2006), f t ∼ GP (m, k) , with mean function m(·) and kernel function k(·, ·). Conditioned on the data, the posterior distribution is a GP again with the posterior mean and variance at a query point x * where X t = (x 1 , . . . , x Nt ) and y t = (y 1 , . . . , y Nt ) is the vector of corresponding, noisy observations. Given this belief over the function f t , BO aims to actively select a new query point x Nt+1 that is informative about the optimum of f t . This requires trading off exploitation and exploration and is usually accomplished through an acquisition function α that depends on the posterior, Common choices for α are the expected improvement (Jones et al., 1998) and the upper confidence bound (Srinivas et al., 2010), among several others. The dataefficiency of these methods crucially depends on how fast the posterior distribution collapses around the true target function f t .
To enable faster optimization, transfer learning additionally exploits data from n s ≥ 1 related source tasks based on functions f s : D → R. While all methods discussed in the following apply to this general setting, we focus on n s = 1 for ease of exposition. In this case, we have N s additional data points D s = {X s , y s } that can be used to improve the target model. Further, we assume without loss of generality that the source is modelled with a GP, f s ∼ GP (0, k s ), with zero mean prior and arbitrary kernel function k s .

GAUSSIAN PROCESS TRANSFER MODELS
In this section, we provide a unified overview of existing GP models for transfer learning that allow for a closedform posterior distribution.
Multi-Task GP (MTGP) Our starting point is the kernel function k of the joint model of source and target. Let (x, i), (x , j) be two points from tasks i, j ∈ {s, t}. We assume that k is a sum of separable kernels where the dirac-delta function δ ij is equal to one if i = j and zero otherwise, and k ν are arbitrary kernel functions (Álvarez et al., 2012). The positive semidefinite matrices W ν are often referred to as coregionalization matrices, since their diagonal (off-diagonal) entries control correlations within (between) data sets. When counting the hyperparameters, we assume that each kernel function has at least one hyperparameter. For n s > 1 sources, this model has O(n 3 s ) scalar hyperparameters from the matrices W ν , O(n s ) scalar noise hyperparameters σ i , and, in addition, the hyperparameters of the kernels k ν . The computational complexity of training this model with a given set of hyperparameters scales as O(N 3 ), where N is the total number of data points from all tasks. This makes it expensive and challenging to optimize the hyperparameters of the model. In the following, we introduce several simplifications from the literature that have fewer hyperparameters and/or better scaling properties.
Multi-Task-Single-k GP (MTKGP) A common heuristic is to assume that the source and target functions share a common structure and consider k s = k t . The joint kernel contains only one coregionalization matrix, which reduces the number of scalar hyperparamters to O(n 2 s ), while the computational complexity is the same as (2).

Weighted Source GP (WSGP)
Another common simplification is to set correlations between different source data sets to zero The hyperparameters w st quantify how much the source is correlated with the target. Each coregionalization matrix has at most one parameter, which reduces the total parameter number to O(n s ).
Hierarchical GP (HGP) The asymmetric setting in transfer learning, i.e., modelling the target but not the source tasks, motivates a hierarchical approach to model the differences of target to source data with an additive kernel defined by As in WSGP, the number of scalar parameters is of order O(n s ). If optimized jointly, the computational complexity is the same as of MTGP.
Mean Hierarchical GP (MHGP) Golovin et al. (2017) propose a simpler way to transfer information from a pretrained source model to the target and only use the posterior mean of the source model as prior mean for the target. All correlations of the source model are neglected, which leads to The kernel is block diagonal in the different tasks. The training complexity of MHGP is therefore the complexity of training the target model plus the additional cost of evaluating the source posterior mean at the target points. Note that Golovin et al. (2017) additionally combine source and target uncertainty heuristically, which we do not consider for ease of exposition.

INTERMEDIATE-COMPLEXITY TRANSFER MODELS
The Bayesian methods presented above require an expensive optimization of the joint likelihood. In this section, we introduce our novel models, Sequential Hierarchical GP (SHGP) and Boosted Hierarchical GP (BHGP), which leverage the asymmetric setting of transfer learning to lower the complexity.
Sequential Hierarchical GP (SHGP) The starting point is the HGP model from (4) in which both the source, p(θ s |D), and target, p(θ t |D), model parameters are influenced by all task data, D = D s ∪ D t . Here, the quantities θ s and θ t parametrize the prior distribution of the GP. The asymmetric nature of transfer learning, in which the target data are accumulated during optimization while the source data remain invariant, motivates the use of a modular Bayesian approach in which the target data do not influence the source model (Bayarri et al., 2009). Training such a model involves (i) finding a suitable p(θ s |D s ) by optimizing the partial likelihood, L s , and (ii) finding p(θ t |D) by optimizing the total likelihood while keeping L s fixed to the value from the first step. This sequential training procedure is the telltale of SHGP. Under more restrictive assumptions and scope, a similar approach was proposed in Kennedy and O'Hagan (2000); Le Gratiet and Garnier (2014). The modular nature of SHGP leads to a large decrease in training complexity compared to HGP as shown later. The price to pay for the complexity reduction is the inability of the source model to react to the target data, which may become an issue in case p(θ s |D s ) is misspecified.
We focus, for concreteness, on hierarchical models in this paper but note that such modularizations, which decrease the training complexity, can, in principle, be conducted for other Bayesian designs from Sec. 4, too. Exploring such approaches is left for future work.
Boosted Hierarchical GP (BHGP) In SHGP, the source model parameters, θ s , do not depend on the target data but the source posterior distribution does. We propose BHGP, which makes also the latter independent of the target data. This ad-hoc assumption takes BHGP outside the Bayesian realm but provides an elegant connection to boosting, a well-known approach in machine learning that combines an ensemble of weak learners into a strong one (Schapire, 2003). We study this connection in Sec. 6. In contrast to HGP, this leads to an asymmetric treatment of observed tar-get data and query points. In Sec. 6.2 we show that this approach results in the kernel Note that BHGP introduces an additive term Σ boost * = Σ s * , * + α * ,t Σ s t,t α T * ,t − α * ,t Σ s t, * − Σ s * ,t α T * ,t to the covariance of the query points, which resembles a robustified version of the target model. Here α * ,t = k t (x * , X t ) k t (X t , X t ) + σ 2 t I −1 , and Σ s t,t , Σ s t, * , Σ s * ,t , Σ s * , * are the blocks of the posterior covariance matrix of the source evaluated at target and query points. The number of hyperparameters is equal to MHGP, while the computational complexity is the same as of SHGP.

THEORETICAL FOUNDATIONS
In Secs. 4 and 5 we introduced several transfer-learning models that rely on the hierarchical kernel architecture. In the following, we present insights on the connection between these methods and their design choices.
Mean Hierarchical GP is the simplest of the models. It trains a GP on the source data and propagates the posterior mean to be the prior mean of the target model. This approach has the lowest computational complexity, see Table 1 (MHGP). Neglecting the transfer of uncertainty from the source to target comes at a cost, as shown in the experimental section (Sec. 7), since well-calibrated uncertainty estimates are at the core of BO's sample efficiency and neglecting the uncertainty may be detrimental to the optimization. So far we have adopted a Bayesian view when presenting the transfer models in Secs. 4 and 5. This view is, however, not particularly useful for the non-Bayesian nature of BHGP. We therefore switch to an alternative and unifying view in which uncertainty propagates via the stochastic realizations of the models. In particular, we study the creation of an ensemble of target models based on functional samples from the source posterior. Averaging over the ensemble would then lead to the desired target model. In this section, we show that both our proposed transfer-learning models, SHGP and BHGP, are closed-form solutions to such averaging procedures.
Performing model averaging over the ensemble of target priors results in the target model of Hierarchical GP, as we show in Sec. 6.1. In light of the ensemble averaging, this kernel naturally lends itself to a sequential optimization of hyperparameters resulting in our proposed Sequential Hierarchical GP. The key advantages Table 1: Computational complexity for one source task. Training the source model has a complexity of O(N 3 s ). Querying the target has a complexity of O(N 2 t + N s ) for MHGP and O((N t + N s ) 2 ) for the other techniques. In Appendix B we discuss the generalization to multiple sources.

Models
Abbreviation HPs Joint HPO Training of target of this approach are the accurate uncertainty estimates with reduced computational complexity, see Table 1 (HGP vs SHGP), and the sequential optimization of weakly correlated subsets of hyperparameters, which stabilizes training and may improve performance, see Before diving into the theoretical analysis, we illustrate the fundamental differences between the three hierarchical models in Fig. 1. In the plots, we train the techniques on data generated from an Alpine function family, f (x; c) = x sin(x + π) + cx, with one input, x ∈ (−10, 10), and one parameter defining the family, c ∈ R. The source data are generated uniformly within x ∈ (−10, 0) with c = 1/2, which is why the source posterior distribution is uncertain on the right, x ∈ (0, 10), see Fig. 1(a). The target data are generated with c = −1/2 for X t = (1, 2, 3, 4). The posterior of Mean Hierarchical GP underestimates the uncertainty of the target function on the right, since its uncertainty originates solely from the target data, see Fig. 1(b). Boosted Hierarchical GP alleviates this shortcoming on the right, while on the left it has the same uncertainty since the source model is confident. The posterior of Sequential Hierarchical GP is Bayesian and correctly captures the variation of the target function on the right, see Fig

Bayesian model averaging
In this section, we show that Bayesian model averaging over prior mean functions of the target, which are distributed according to the posterior of the source, leads to the HGP model defined by (4).
Then, the joint model of f s and f bayes t is a GP with zero mean prior and the kernel (4).
We provide a closed form solution for the marginalization of (7) over the source posterior in Appendix A. Since both integrands are Gaussian, the integral yields a Gaussian distribution for p(f bayes t | D t ). By conditioning the joint model from Proposition 1 on the source data, we can derive the following prior for the target model.
We provide the detailed derivations in Appendix A. The training complexity scales as O( .5 for the proof. We therefore save, during training, the steep complexity contribution of O(N 3 s ) of the conventional Bayesian methods.

Posterior prediction averaging
In an alternative approach we take inspiration from the well-known principle of boosting (Schapire, 2003) and average over the posterior distributions of the target models. In contrast to Bayesian model averaging in Sequential Hierarchical GP, this approach neglects the implicit dependency of the source model on the (7). This distribution is analytically tractable and equals the posterior distribution of a GP with kernel (6), i.e., the Boosted Hierarchical GP.
Proposition 2. Let f s , f t be as in Proposition 1 and Then, the joint model of f s and f boost t is a GP with zero mean and the covariance function defined in (6).
The proof follows the same ideas as for Corollary 1 and can be found in Appendix A. Since the prior for the target of BHGP coincides with the prior of MHGP, the hyperparameter optimization is the same. The training complexity is the same as for SHGP,

EXPERIMENTS
We provide an experimental evaluation of the transferlearning algorithms discussed in Table 1. In addition, we evaluate three baselines: GP-based BO (GPBO) without any source data, RGPE (Feurer et al., 2018) in which the target function is modelled as a weighted sum of the predictions of all task GPs, and GC3P (Salinas et al., 2020) in which a Gaussian Copula regression with a parametric prior is employed to scale to large data. We implement the models using GPy (GPy, since 2012) and run BO with Emukit (Paleyes et al., 2019), licensed under BSD 3 and Apache 2.0, respectively. For GC3P we use the publicly available code (Salinas et al., 2020). GP hyperparameters are optimized by maximizing the likelihood for the observed data. We focus on maximum likelihood for simplicity but, in high dimensions, it may be preferable to employ a Bayesian treatment of hyperparameters. Code to reproduce our results can be found on Github 1 .

Synthetic Function Families
The synthetic function families that are derived from conventional benchmark functions by placing probability distributions on their parameters. We consider a broad spectrum of families, where the optimum of tasks varies locally for some and in the entire domain for others: the one-dimensional Alpine (Momin and Yang, 2013) function, and the multi-dimensional Branin, Hartmann3 and Hartmann6. Details about the function families are provided in Appendix C.1. Performance is reported via mean simple regret over 50 runs together with standard error of the mean.

Single-source benchmark functions
We start the experimental analysis in the simplest setting with a single source task. This fundamental regime provides valuable insight into the trade-off between performance and scalability of the transfer learning algorithms. Applications with scarce historical data where the transfer efficiency is important may particularly benefit from this analysis.
The algorithms are benchmarked on three function families, see Fig. 2. Results on two further function families are available in Appendix D.1. A notable difference in the convergence on the different benchmarks can be observed. This is likely to be caused by the different particularities of each function family, which we discuss in Appendix C.1. Despite these differences, there are a few important generalizing features that characterize the algorithms: (i) The Bayesian techniques SHGP, HGP, and WSGP outperform the non-Bayesian algorithms. This is not surprising, since Bayesian models compute highly informed probabilistic predictions. The best and most consistent performers are the HGP and WSGP algorithms. We attribute this to the flexibility of their design. During training, all kernel's hyperparameters are jointly optimized allowing for superior model quality. (ii) Our modular Bayesian technique, SHGP, is competitive and an attractive alternative with lower computational complexity. The other non-Bayesian methods trail behind, which is likely due to their inability to reliably propagate model uncertainty. Among these, our technique, BHGP, performs consistently better than other non-Bayesian techniques  (left), three-dimensional Hartmann3 (middle), and six-dimensional Hartmann6 (right). The source data are sampled randomly from the source function and contain 20N dim points with N dim being the input dimension of the benchmark. We add i.i.d. observational noise of standard deviation σ s = σ t = 0.1 for Hartmann3, Hartmann6, and σ s = σ t = 1.0 for Branin during data generation. and provides a reasonable compromise between computational complexity and efficiency. (iii) Surprisingly, the most advanced and expensive techniques, MTKGP and MTGP, perform worse than the other Bayesian techniques. This is likely caused by the challenging training procedure, where a large number of hyperparameters are optimized.
In addition to the asymptotic complexities in Table 1, we demonstrate the lower computational complexity of our methods in a direct runtime analysis of model training in Fig. 3(a). SHGP and BHGP are orders of magnitude faster than the full Bayesian methods for moderately large source data sets. The runtimes reported in Fig. 3(a) are representative of the entire optimization runtime since model training has steeper complexity than acquisition-function optimization, see Appendix D.2 for empirical evidence.
Observational noise and propagation of uncertainty Beyond these generic trends, the performance of the techniques depends on other factors like the amount of source data and the magnitude of the observational noise. Such a study is presented in Fig. 3(b, c); more results are available in Appendix D.3. We distinguish three data regimes: (i) In the limit of scarce source data, which is insufficient to describe the global shape of the source function, the kernel methods with joint HPO, HGP and WSGP, outperform the other methods significantly. The reason is that they model all data jointly, in contrast to the methods trained sequentially that end up with source models of poor quality. (ii) The case of moderate amount of source data, which is sufficient to describe the source function probabilistically but not deterministically, is discussed in the previous section. (iii) In the limit of lots of source data, which is sufficient to deterministically describe the function, all transfer-learning methods converge to a similar performance. Here, propagation of uncertainty is not required and the non-Bayesian approaches are appealing due to their favorable scaling.
The observational noise affects the boundary between the aforementioned data regimes. Enhanced noise leads to an increased amount of data required by a model of fixed quality.

Multi-source benchmark functions
For the more general case of multiple source data sets, the algorithms are benchmarked on two function families: (i) the six-dimensional Hartmann6 family with three source data sets in which the functions are sampled randomly as in Sec. 7.1, and (ii) the one-dimensional Alpine family with five source data sets, where the source and target functions are fixed as in (Feurer et al., 2018). MTGP is not benchmarked here because of its complexity. The performance of the algorithms on the three-source Hartmann6 is consistent with the single-source benchmarks, compare Figs. 2 and 4. HGP and WSGP still perform best, while our less complex SHGP provides a competitive compromise. The picture changes in the Alpine benchmark, where the more structured approaches to HPO show superior performance. Our techniques, SHGP and BHGP, outperform others despite the lower computational complexity. This is likely due to the stable,   Figure 4: Performance on the multi-source benchmarks: the three-source Hartmann6 (left), five-source Alpine (middle), and ten-source OpenML SVM benchmark (right). The source data are sampled randomly and contain 60, 20, and 60 points per task for Hartmann6, Alpine, and SVM, respectively. I.i.d. observational noise of standard deviation σ s = σ t = 0.1 is added to the source and target data of the synthetic benchmarks. hierarchical training procedure, where at most three hyperparemeters are optimized at a time. By contrast, HGP optimizes 18 hyperparameters jointly, WSGP 23, and MTKGP 33. Algorithms employing hierarchical hyperparameter optimization are therefore ever more appealing for increasing number of source tasks.

Meta-Learning Surrogate Benchmarks
In the following we study the performance of the transfer-learning techniques on hyperparameter optimization problems in machine learning. We follow Perrone et al.  2020), we carry out a discrete optimization at the observed evaluations to avoid potential problems related to the bias of a fitted surrogate model.
Since the data sets contain over fifty thousand evaluations, which is too much for the closed-form Bayesian techniques to handle, we perform downsampling and choose 11 tasks at random. The data of the target task are used as the ground truth for the discrete optimization, while the data of the source tasks are downsampled to 60 points per task. We present results averaged over 500 independent runs. For better comparability, we employ a rescaled version of the regret function called average distance to the global minimum (Wistuba et al., 2016), where the regret values on each data set are rescaled between zero and one before the averaging procedure. The results are shown in the right panel of Fig. 4.
Despite the markedly different nature of the data compared to the synthetic benchmarks, the performance of the different transfer-learning techniques exhibits somewhat similar trends. The performance is dominated by WSGP, which achieves achieve a regret of 0.01 more than three times faster than classical GPBO. Interestingly, HGP performs worse. We hypothesize this to be caused by the incompatibility between the hierarchical nature of HGP and the structure of the SVM data set. Our techniques, SHGP and BHGP, provide, as before, a tradeoff between complexity and performance, and achieve the same regret level about twice as fast as GPBO.

DISCUSSION AND CONCLUSION
In this paper, we have presented a unified view on transfer learning methods based on Gaussian processes with tractable model posteriors. One end of this spectrum is populated by models that maintain a full covariance matrix of all source and target data and perform proper Bayesian inference on all data jointly. While these are powerful models with well-calibrated uncertainty estimates, they quickly become computationally infeasible due to their cubic complexity in the cumulated number of source and target data points. The other side of the spectrum is occupied by heuristics that aggregate the predictions of separate, independently trained, task models. These methods exhibit a much more favourable computational complexity, but unfortunately fail to properly propagate uncertainty between the individual task models. To amend this unsatisfactory binary choice between computational feasibility and uncertainty propagation, we proposed two novel transfer learning models, Sequential Hierarchical GP and Boosted Hierarchical GP, that both provide a compromise between these extremes: they add some computational complexity to account for the uncertainty of the involved models with some slight approximations as compared to the full covariance approaches. Our analysis is supported by comprehensive experiments that pinpoint strengths, weaknesses, and trade-offs of these transfer learning methods. The benchmarks demonstrate the appeal of Sequential Hierarchical GP and Boosted Hierarchical GP as robust and competitive techniques.

Limitations
The low-data regime is challenging for any method. In particular, the empirical success of Gaussian processes hinges on the smoothness assumptions encoded in the kernel. At the same time, it is not obvious how nonsmooth functions can be optimized efficiently. Another limitation is that they do not scale to large scale data without approximations, which is not a problem in our setting since we specifically focus on the low-data regime. We are not aware of any negative societal impacts of our work, since we focus on improving the data efficiency of optimization methods.

Supplementary Material: Transfer Learning with Gaussian Processes for Bayesian Optimization Overview
In the appendix we provide the detailed proofs for all claims in the paper, complexity analysis for all discussed methods for multiple sources, ablation studies, and the detailed hyperparameter configurations used in the paper. An overview is shown below.

A EQUIVALENCES
In the following sections, we provide the proofs of the claims in the main paper. For ease of reference, we restate those claims here. We begin with deriving closed form solutions for the prior and posterior of Bayesian model averaging. and Then, the joint model of f s and f bayes t is a GP with zero mean prior and the kernel (4).
Then under the assumptions of Proposition 1 it holds that Subsequently, we show that posterior prediction averaging leads to a analytically tractable distribution which equals the Boosted Hierarchical GP (6).
Proposition 2. Let f s , f t be as in Proposition 1 and Then, the joint model of f s and f boost t is a GP with zero mean and the covariance function defined in (6). Further, The analysis of both ensemble averaging techniques introduced in Sec. 6 is based on marginalization of the models over the posterior of the source at target data points X t and query points x * . Samples from the source posterior can be obtained from the posterior mean E (f s | D s ) and covariance of the source cov (f s | D s ) posterior mean and covariance of the source model as follows.
Let L be the Cholesky decomposition of posterior covariance of the source LL T = cov f s |X t , D s evaluated at target data points and query pointsX t = X t ∪ x * . Let ε be a standard normal vector of size N t + 1 then is a sample from the source posterior at target data and query points. The following commonly known result will be useful for deriving the desired Gaussian distributions. For the sake of completeness, we provide a proof of Lemma 1 in Appendix A.4.
Lemma 1. Let ε be a standard normal random vector of size N × 1.
We begin with the proof of Corollary 1, which is a simple application of Lemma 1 to the ensemble of priors used in Bayesian model averaging. We derive Proposition 1 from Corollary 1 by showing that conditioning the joint model from (4) leads to the same posterior. The proof of Proposition 2 requires a more involved application of Lemma 1 to the ensemble of posteriors, that we propose in posterior prediction averaging. A detailed derivation of Lemma 1 can be found in Appendix A.4.

A.1 Proof of Corollary 1
In Bayesian model averaging, we build an ensemble of priors for the target model from multiple functional samples of the source posterior. For a fixed functional sample f s,ε (X t ) the prior of the target is a multi-variate normal with mean E f s |X t , D s + Lε and covariance matrix k t (X t , X t ). Consequently, Corollary 1 follows from Lemma 1 and LL T = cov f s |X t , D s .

A.2 Proof of Proposition 1
To show Proposition 1, we start with the model defined in (4) and condition it on the source data. The joint kernel of the Hierarchical GP kernel is given by and after conditioning this model on source data D s , we obtain the prior of the target model by standard GP inference. Then, As a consequence, the posterior of the Hierarchical GP model given the source data has mean and covariance We can rewrite Together with Corollary 1, we obtain that conditioning the Hierarchical GP model defined in (4) on source data leads to the prior target model of Bayesian model averaging. Proposition 1 follows because under the same modelling assumption for the source, the prior is unique.

A.3 Proof of Proposition 2
In this subsection, we compute the distribution obtained from directly averaging over the posteriors on the target data induced by the mean priors sampled from the source posterior (and, thus, ignoring the implicit dependency of the source model on the target data). Denote by L s, * and L s,t the first N * and remaining N t rows of L s . For a fixed sample f s,ε (X t , x * ), the posterior of the target GP at x * is a multivariate normal with mean and covariance matrix The posterior mean can be rewritten as with L boost = L s, * + α * ,t L s,t . Recall, that Σ s t,t , Σ s t, * , Σ s * ,t , Σ s * , * denote the blocks of the posterior covariance matrix of the source evaluated at target and query points. Then it holds that From this, we obtain With this observation, the posterior distribution stated in Proposition 2 can be derived from Lemma 1.
Finally, it remains to show that the joint model defined in (6) leads to the same posterior distribution. Since the source kernel is only evaluated at source points, conditioning on the source first yields Note that the additive term Σ boost * only contributes to the covariance at query points. That is, it will only contribute an additive term to the posterior covariance. Using standard formulas for conditional distributions in multi-variate Gaussians, we obtain that f boost

A.4 Proof of Lemma 1
Let p (x; µ + Lε, Σ) denote the density of a multi-variate normal distribution with mean µ + Lε and covariance matrix Σ. Further, denote by |Σ| the determinant of Σ. Then, Similarly as above, denote by p (x; 0, 1) denote the density of a spherical multi-variate Gaussian. Then from (22), we obtain From this the result follows by integrating out ε

B GENERALIZATION TO MULTIPLE SOURCES
In this appendix, we discuss the generalization of Sec. 4 to multiple sources. We discuss the computational complexity to predict mean and covariance of the target model at a new query point x * . We collect all calculations that can be done without knowledge of the query point as "training complexity" and the remainder as "prediction complexity". In this section we assume n s > 1 sources. For easier generalization we denote the sources by ν = 1, 2, ..., n s and the target by ν = n s + 1. Where appropriate, the analysis is structured in four parts: Training of the first source, training of the n > 1th source assuming all "lower" sources n < n have already been trained, training of the target assuming all sources have been trained, and prediction from the target model. Our basis operation for evaluating the computational complexity is scalar multiplication. Sometimes, we will also present our results using N t = N ns+1 and N s = ns ν=1 N ν .

B.1 Multi-Task GP (MTGP)
The joint kernel for n s sources and the target in Multi-Task GP between two points (x, i), (x , j) from tasks i, j ∈ {1, 2, ..., n s + 1} is a sum of separable kernels given by The number of hyperparameters is O(n 3 s ). Training and prediction complexity are O((N t + N s ) 3 ) and O((N t + N s ) 2 ), as for a simple GP with N t + N s data points.

B.2 Multi-Task-Single-k GP (MTKGP)
The only difference to Multi-Task GP is that k ν = k s , which reduces the number of hyperparameters to O(n 2 s ).

B.3 Weighted Source GP (WSGP)
The joint kernel for n s sources and the target in Weighted Source GP is given by Since all the blocks between different sources are zero, inverting the kernel matrix K is computationally simpler than for the more general kernels. We define the source source-block (which is diagonal in i, j) A = [K] i≤ns,j≤ns , the target block D = [K] ns+1,ns+1 , and the connecting blocks B = [K] i≤ns,j=ns+1 and C = [K] i=ns+1,j≤ns . With this definition, we see that using block-inversion (Lu and Shiou, 2002).
we can use the block-diagonal structure of A for an efficient inversion of the kernel matrix. The total complexity is O(N 3 t + N 2 t N s + N t N 2 s + ns ν=1 N 3 ν ). Note how due to the blockdiagonal structure of A, the term ( ns ν=1 N ν ) 3 for the inversion of a matrix of the size of A collapses to the much more favourable ns ν=1 N 3 ν . Prediction at a query point is of order O((N t + N s ) 2 ).

B.4 Hierarchical GP (HGP)
The joint kernel for n s sources and the target in HGP is defined by (25) with [W ν ] i≥ν,j≥ν = 1 and 0 otherwise .
The number of hyperparameters is hence the number of hyperparameters in the kernel functions, which is of order O(n s ). When optimizing the hyperparameters of all sources and the target jointly, the training and prediction complexity are O((N t + N s ) 3 ) and O((N t + N s ) 2 ), as for a simple GP with N t + N s data points.

B.5 Sequential Hierarchical GP (SHGP)
SHGP carries out a sequential training, one source at at time. In the following we analyse the complexity of these sequential steps individually.
Training of the first source The first source is a simple GP, with training complexity O(N 3 1 ) and prediction complexity O(N 2 1 ). Training of the nth source In order to train the nth source, we need to evaluate the posterior mean and covariance of the n − 1th source at the data points X n . This requires evaluating the posterior means of all underlying sources at the same points, resulting in O(N n n−1 ν=1 N ν ) multiplications. In addition we need to invert the covariance matrix Σ n (X n , X n ) + σ 2 n I which is of order O(N 3 n ). In order to obtain the covariance matrix, we need to evaluate the posterior covariance of the n − 1th source X n Σ n (X n , X n ) = k n (X n , X n ) + Σ post n−1 (X n , X n ), Σ post n−1 (X n , X n ) = Σ n−1 (X n , X n ) − Σ n−1 (X n , X n−1 ) Σ n−1 (X n−1 , X n−1 ) + σ 2 n−1 I −1 Σ n−1 (X n−1 , X n ) Inverting Σ n−1 (X n−1 , X n−1 ) + σ 2 n−1 I is part of the training of the n − 1 source.
The later has complexity O( Hence, we actually need the complexity of Σ post ). All together we have a training complexity of Training of the target Training of the target is analogous to training the nth source. Setting n = n s + 1, we have complexity O(N 3 t + N 2 t N s + N t N 2 s ). Prediction from the target The prediction at a new query point x * is dominated by evaluating the covariance matrix. We can read-off this complexity from (31) by ignoring the N 3 n term and setting n = n s + 2 and N ns+2 = 1, giving O((N t + N s ) 2 ).

B.6 Mean Hierarchical GP (MHGP)
The joint kernel for n s sources and the target in Mean Hierarchical GP is defined by (25) with The number of hyperparameters is hence the number of hyperparameters in the kernel functions, which is of order O(n s ).
Training of the first source The first source is a simple GP, with training complexity O(N 3 1 ) and prediction complexity O(N 2 1 ). Training of the nth source In order to train the nth source, we need to evaluate the posterior mean of the n − 1th source at the data points of the nth source. This in turn requires the evaluating the posterior means of all underlying sources at the same points, all together O(N n n−1 ν=1 N ν ) multiplications. Except for this, the complexity is the same as in a simple GP with N n data points, hence all together we get a training complexity of O(N 3 n + N n n−1 ν=1 N ν ). Training of the target Training of the target is analogous to training the nth source. Setting n = n s + 1, we have complexity O(N 3 t + N t N s ). Prediction from the target In order to predict mean of the target model at a new query point x * , we need to evaluate the posterior means of all sources at x * , requiring O(N s ) multiplications. For predicting the covariance, the complexity is the same as for a simple GP O(N 2 t ). The overall prediction complexity is hence O(N 2 t + N s ).

B.7 Boosted Hierarchical GP (BHGP)
To the best of our knowledge in Boosted Hierarchical GP there is no unique generalization to n s > 1 sources, because of the heuristic aggregation of uncertainty. Below we outline several avenues that one could consider, which differ in terms of complexity. The most straightforward approach is to do everything as in Mean Hierarchical GP, and add a boosting term to the posterior prediction of the target task. The kernel for that is simply given by In this single layer boosting, we have [W ν ] ij = δ iν δ jν , k ns+2 = k ns+1 +Σ boost ns+1 , and σ ns+2 = 0. Here ν, i, j = n s +2 describes the query points and In that case, the boosting is only done in the last layer. The complexity of Boosted Hierarchical GP has two contributions. One, the complexity of Mean Hierarchical GP and the complexity of calculating the boosting term. For the boosting term, we need to calculate Σ post ns (X ns+1 , X ns+1 ) which is of order O N 2 ns+1 N ns + N ns+1 N 2 ns , this can be done during training, since it does not depend on the query point. During prediction we need to evaluate α ns which is of order O N 2 ns+1 because the matrix inversion can be done as part of the training. We also need to to evaluate Σ post ns (x * , x * ) and Σ post ns (x * , X ns+1 ) which is of order O N ns+1 N ns + N 2 ns . To conclude the calculation of the boosting term, we need to evaluate (34) which is of order O N 2 ns+1 . All together we have training complexity O ns+1 ν=1 N 3 ν + N 2 ns+1 N ns + N ns+1 N 2 ns and prediction complexity O N 2 ns+1 + N ns+1 N ns + N 2 ns + ns ν=1 N ν ).
An alternative approach which we follow in our experiments is to add such boosting terms at each intermediate prediction in the hierarchy by substituting Σ post ns → Σ post ns + Σ boost ns recursively in (34). Note that Σ post 1 = 0. Note that the likelihood of the boosted model is the same as of Mean Hierarchical GP. Thus, this change only affects the prediction. However, before being able to make predictions, we now have to pre-calculate additional terms. In particular, this requires to recursively evaluate the posterior covariances, very similarly as in (30). This leads to an additional training complexity of O(N 2 t N s + N t N 2 s ). During prediction, the additional cost comes from the boost term which involves recursive evaluation of the posterior covariance terms containing the query point. This increases the total prediction cost to O((N t + N s ) 2 ).

C DETAILS ON EXPERIMENTS
All the experiments use default BO parameters with Upper Confidence Bound (UCB) and exploration coefficient β = 3 as acquisition function. The source and target functions are sampled randomly from the function family, and the source data are sampled randomly from each source task, see Appendix C.1 for details. BO is used to optimize the target function from scratch, i.e., without initial samples. A certain amount of i.i.d. zero-mean Gaussian observational noise is added during the data generation process with a standard deviation specified in each figure caption.
All GP models are based on the squared-exponential kernel with automatic relevance determination. The GP hyperparameters are optimized by maximizing the likelihood function for the observed data using the L-BFGS-B optimizer with 10 initial guesses, x 0 , where each instance is sampled from This generic optimization strategy is possible since we normalize the observed data to zero-mean, unit-variance. All other details regarding training and prediction are kept fixed to GPy's defaults (GPy, since 2012). The only exception is the OpenML SVM benchmark, where, for complexity reasons, we only use 3 initial guesses and a maximum of 20 gradient-ascent steps for the optimization of the likelihood function of HGP and WSGP.

C.1 Details on the Synthetic Function Families
The Forrester Family The Forrester function is a one-dimensional synthetic function with one global minimum, one local minimum, and a zero-gradient inflection point. The function is given by The original function is given by a = 1, b = 0, c = 0. In this paper, a family of functions is formed by choosing the following probability distributions for the parameters (a, b, c): where U(·, ·) denotes the uniform distribution. The Forrester family is therefore defined over a three-dimensional uniform distribution.
For generating the data of n s source data sets, we draw n s random tasks using (38), and sample a given number of points per task using a random distribution. In Fig. 2, 20 points per task are sampled.
The Alpine Family The Alpine function is a one-dimensional synthetic function with one global minimum and three local minima. The function is given by f (x; s) = x sin(x + π + s) + 0.1x, x ∈ [−10, 10].
The original function is given by s = 0. In this paper, the Alpine function is used for the multi-source experiments in Fig. 4. There, s = 0 is used for the target function, while the five source functions are given by s = kπ/12, k = 1, ..., 5. Twenty random points are sampled for each source function. Note that this benchmark is identical to the one presented in (Feurer et al., 2018, Sec. 5.1).
In the single-source experiments presented in Fig. 2, the target function is also set to s = 0, while the source function is sampled randomly from the five choices above.
The Branin family is therefore defined over a six-dimensional uniform distribution. For generating the data of n s source data sets, we draw n s random tasks using (41), and sample a given number of points per task using a random distribution. In Fig. 2, 40 points per task are sampled.
The Hartmann3 family is therefore defined over a four-dimensional uniform distribution. For generating the data of n s source data sets, we draw n s random tasks using (43), and sample a given number of points per task using a random distribution. In Fig. 2, 60 points per task are sampled.
The Hartmann6 Family The Hartmann6 family is a sum of four six-dimensional Gaussian distributions and is defined by The original Hartmann6 function is given by α = (1.0, 1.2, 3.0, 3.2) T . In this paper, a family of functions is formed by choosing the following probability distributions for the parameters α = (α 1 , α 2 , α 3 , α 4 ) T : The Hartmann6 family is therefore defined over a four-dimensional uniform distribution. For generating the data of n s source data sets, we draw n s random tasks using (45), and sample a given number of points per task using a random distribution. In Fig. 2 and Fig. 4, 120 and 60 points per task are sampled, respectively.

D.1 Results on one-dimensional benchmark functions
Results on the single-source one-dimensional Forrester and Alpine function families are presented in Fig. 5. The function families are described in Appendix C.1. These results are consistent with the discussion and analysis from Sec. 7.1.

D.2 Runtime Analysis
Here we present runtime studies for the entire optimization procedure, which includes (i) training the models including optimizing their hyperparameters, and (ii) optimization of the acquisition function. We time the runs by optimizing the Hartmann6 benchmark function. We consider one single source task with the same configuration as in Fig. 2. Results are shown in Fig. 6 and are consistent with training runtimes reported in Fig. 3. This reinforces that model training is more expensive than acquisition-function optimization since the former has steeper complexity than the latter for the Bayesian models, see Tables 1 and 2.

D.3 Impact of Observational Noise and Propagation of Uncertainty: Ablation Studies
We perform a systematic experimental study on the impact of observational noise and amount of source data on the algorithm performance. The three regimes discussed in Sec. 7.1 can be clearly observed for the Hartmann3 function family in Fig. 3.
Here, we showcase the dynamics of these three regimes for the challenging Hartmann6 function family in Figs. 7-11. The first regime, in which the amount of source data is insufficient to build a model that faithfully describes the global shape of the source function, is presented in Fig. 7, where only 30 data points are considered. Here, algorithms that model the data jointly (HGP, WSGP) have a clear advantage. A crossover to the intermediate regime, in which the source data describe the source function much better but with significant uncertainty, happens gradually in Figs. 8-10, where the source data are increased from 60 to 250 points. Here, our relatively lightweight SHGP is competitive and on par in terms of performance with the more general HGP. The third regime is reached in Fig. 11, where the performance of all methods except GPBO and RGPE converge to similar performance.
The observational noise adds a further level of complexity to these dynamics. In essence, the observational noise controls the boundary between the three regimes: the more noise is present, the more data points are required to describe the function for a fixed model quality. In other words, observational noise leads to enhanced model uncertainty, whose propagation is key to obtaining a high-quality model. It is therefore expected that (i) the performance gap between Bayesian and non-Bayesian techniques becomes ever bigger for increased observational noise, and (ii) our developed techniques, SHGP and BHGP, particularly excel for elevated observational noise       Figure 11: Performance on the single-source Hartmann6 function family for 500 points in the source data set. We vary the amount of observational noise, σ s = σ t , from 0.001 to 0.1. The left panel depicts the mean ± standard error of the mean. The right panel shows the median and 20% quantile range.