# HyperBO+: Pre-training a universal prior for Bayesian optimization with
hierarchical Gaussian processes

Zhou Fan
Harvard University
<EMAIL_ADDRESS>
&Xinran Han
Harvard University
<EMAIL_ADDRESS>
Zi Wang
Google Research, Brain Team
<EMAIL_ADDRESS>

###### Abstract

Bayesian optimization (BO), while proved highly effective for many black-box
function optimization tasks, requires practitioners to carefully select priors
that well model their functions of interest. Rather than specifying by hand,
researchers have investigated transfer learning based methods to automatically
learn the priors, e.g. multi-task BO (Swersky et al.,, 2013), few-shot BO
(Wistuba and Grabocka,, 2021) and HyperBO (Wang et al.,, 2022). However, those
prior learning methods typically assume that the input domains are the same
for all tasks, weakening their ability to use observations on functions with
different domains or generalize the learned priors to BO on different search
spaces. In this work, we present HyperBO+: a pre-training approach for
hierarchical Gaussian processes that enables the same prior to work
universally for Bayesian optimization on functions with different domains. We
propose a two-step pre-training method and analyze its appealing asymptotic
properties and benefits to BO both theoretically and empirically. On real-
world hyperparameter tuning tasks that involve multiple search spaces, we
demonstrate that HyperBO+ is able to generalize to unseen search spaces and
achieves lower regrets than competitive baselines.111This is an extended
version of the workshop paper.

##  1 Introduction

Bayesian optimization (BO) is a fundamental and compact framework to study
intelligent decision making systems: how to make sequential decisions on data
acquisition by reasoning about the Bayesian beliefs on a function, with a goal
to optimize this function in a sample-efficient manner. Each decision requires
selecting inputs to the function and acquiring their corresponding function
evaluations.

Despite the fact that reasoning with Bayesian beliefs can be hard and costly,
BO has found its ways towards real-world deployment because, (1) using a
_Gaussian process (GP)_ as the Bayesian belief makes inference much cheaper
and easier in BO, and (2) global optimization problems on expensive black-box
functions need strong sample efficiency and have relatively low requirement on
the speed of decision making for data acquisition. BO with GP beliefs have
demonstrated numerous successes on such global optimization problems,
including hyperparameter tuning in machine learning (Snoek et al.,, 2012;
Wistuba and Grabocka,, 2021; Eggensperger et al.,, 2021; Pineda-Arango et
al.,, 2021; Cowen-Rivers et al.,, 2022), action selection in robotics (Drieß
et al.,, 2017; Wang et al.,, 2017), search in design spaces of drugs (Pyzer-
Knapp,, 2018) and chemicals (Shields et al.,, 2021; Griffiths and Hernández-
Lobato,, 2020), and even designing aircrafts (Lam et al.,, 2018). However, we
can not credit these accomplishments solely to BO, whose efficacy relies
heavily on expert knowledge that specifies the initial beliefs on their
functions of interest.

In past literature, researchers have also emphasized the importance of
employing a well-specified _GP prior_ in BO with both empirical and
theoretical evidences (Bogunovic and Krause,, 2021; Wang et al.,, 2018; Kim et
al.,, 2019). Using a poorly designed prior, often due to the lack of
quantitative understandings of function landscapes, can be detrimental to the
performance of BO in applications such as hyperparameter tuning (Schulz et
al.,, 2016).

How do we obtain well-specified priors for BO? While it is difficult for us to
think quantitatively and somehow come up with a reasonable GP out of our mind
(Wilson et al.,, 2015), we can sort our past experience with related functions
and partition the observed data on those functions to employ “prior learning”
with i.i.d. sets of non-i.i.d. points (Baxter,, 1996; Minka and Picard,, 1997;
Wang et al.,, 2018). HyperBO (Wang et al.,, 2022), in particular, showed that
we can _pre-train_ Gaussian process priors on related functions efficiently
and effectively for BO applications.

Wang et al., (2022), however, assumed that HyperBO must pre-train on functions
with the same domain, i.e. the same BO search space, making it impossible to
pre-train on observed data from functions with different domains, let alone
generalizing to optimization problems with different search spaces.
Essentially, the pre-trained GP prior in HyperBO is restricted to modeling
functions that have a fixed input dimension - hence the prior is not operable
when we optimize a new function with a different input dimension.

In this work, we aim to pre-train a universal prior that works unanimously for
BO with any search spaces. Our new method, HyperBO+, is based on the
observations that practitioners commonly use stationary kernels in real-world
applications (Balandat et al.,, 2020; Golovin et al.,, 2017; Cowen-Rivers et
al.,, 2022) and that the same prior distribution can be shared by some
parameters, e.g. length-scales of each input dimension. Even if the function
of interest is non-stationary or the output is non-Gaussian, one can employ
input or output warping strategies to ensure GPs have the model capacity
(Cowen-Rivers et al.,, 2022).

In HyperBO+, we partition training data, in the form of past observations, by
relevant search spaces and functions in each search space. We then pre-train a
hierarchical GP to obtain our universal prior which models functions with
different domains. Note that we only need to run the pre-training algorithm
once and we can use the same pre-trained hierarchical GP for BO on any
functions, as long as their search spaces are relevant to the search spaces in
our training data. Theoretically, we show that our pre-trained prior is
consistent with the ground-truth distribution of the functions, a fundamental
milestone for BO to go beyond the Bayesian cradle. Empirically, we demonstrate
the promising results of HyperBO+ to generalize over different search spaces
and its competitive performance on classic multi-task problems.

To the best of our knowledge, our work is the first transfer learning approach
that is capable of doing optimization on unrestricted search spaces in the
literature of BO with GP-based priors. Our contributions include (1) a new
two-step pre-training method for hierarchical GPs that can model functions
with different domains; (2) HyperBO+, a novel transfer learning BO method that
generalizes to universal search spaces; (3) theoretical and empirical analyses
of the asymptotic properties of our pre-training method; and (4) robust
results showing the capability of HyperBO+ to generalize to unseen search
spaces and its superior BO performance through experiments on synthetic
datasets and real-world hyperparameter tuning tasks.

##  2 Related work

Our work is inspired by HyperBO (Wang et al.,, 2022), which showed pre-
training helps learning a better GP prior and improves BO performance on
different functions with the same domain. Similar ideas (Wistuba and
Grabocka,, 2021; Perrone et al.,, 2018; Wang et al.,, 2018; Volpp et al.,,
2020) have also been proposed for transfer learning and meta learning in
Bayesian optimization but they can only transfer the knowledge among functions
with the same domain.

Recently, Chen et al., (2022) introduced OptFormer, a transformer based
transfer learning method for hyperparameter tuning on universal search spaces,
a goal shared by our work. However, due to its core idea of proposing
hyperparameters in an end-to-end fashion, OptFormer must be trained on
millions of BO trajectories, requires giant transformer models, demands
expensive TPU hardware and works only with continuous domains. On the
contrary, our model does not have any constraints on how the data is
collected, uses simple GP models with a small number of parameters, trains
with much shorter time on CPUs and works with both continuous and discrete
domains. Due to practical reasons, we do not compare to OptFormer empirically.

Notably, our work focuses on pre-training a universal prior, a critical and
specific component of BO. Our pre-training method can be naturally coupled
with other important components or heuristics to complete a BO software, e.g.
acquisition functions, input and output warping, cross validation etc. Whereas
frameworks such as BoTorch (Balandat et al.,, 2020), Vizier (Golovin et al.,,
2017) and OptFormer (Chen et al.,, 2022) are hyperparameter tuning software
packages. The pre-training method in HyperBO+ can be directly incorporated in
BoTorch and Vizier by replacing their default hierarchical GP model.

The theoretical soundness of our approach relies on the quality of the
estimated GP parameters during pre-training. The asymptotic behavior for type
II maximum likelihood (ML-II) estimators of covariance parameters under a
single GP has commonly been studied in two asymptotic frameworks (Bevilacqua
et al.,, 2019; Zhang and Zimmerman,, 2005). Under the fixed domain setting,
observations are sampled in a bounded set and thus become increasingly dense
with more samples. Previous work (Kaufman and Shaby,, 2013) shows that only
the microergodic parameters can be estimated consistently and provides
rigorous proof for isotropic Matérn covariance functions under dimension $\leq
3$. A covariance parameter is microergodic if, when taken two different
values, the corresponding Gaussian measures are orthogonal and non-
microergodic when the two are equivalent (Stein,, 1999).

In the increasing domain setting, observations are collected with a minimal
spacing in between data points. Previous results from Mardia and Marshall,
(1984) and Stein, (1999) showed that MLE for covariance parameters are
consistent and asymptotically normal in general with mild regularity
conditions for increasing domains. According to Bachoc, (2014), the key
difference with increasing domain setting is that there is vanishing
dependence between observations that are far apart. Thus, larger sample size
is more informative of the covariance structure. Our theoretical analysis
shows that our problem formulation belongs to the increasing domain setting
and previous theoretical guarantees can be used to ensure the performance of
our pre-training approach.

##  3 Problem formulation

Our goal is to optimize unseen black-box functions by pre-training on existing
data from functions in multiple search spaces, which can have different
numbers of dimensions for their respective domains.

We use the term super-dataset to denote all datapoints collected across
multiple search spaces, while a dataset includes the data from a single search
space, and a sub-dataset means the collection of datapoints from a single
function within a search space.

More formally, we define a super-dataset as $D=\\{D_{i}\\}_{i=1}^{N}$. Each
dataset $D_{i}$ consists of observations on a collection of black-box
functions
$F_{i}=\\{f_{ij}:\mathcal{X}_{i}\rightarrow\mathbb{R}\\}_{j=1}^{M_{i}}$ where
functions in $F_{i}$ share the same compact search space
$\mathcal{X}_{i}\in\mathbb{R}^{d_{i}}$. Let
$D_{i}=\\{D_{ij}\\}_{j=1}^{M_{i}}$, where each sub-dataset
$D_{ij}=\\{(x^{ij}_{k},y^{ij}_{k})\\}_{k=1}^{L_{ij}}$. $L_{ij}$ is the number
of observations on function $f_{ij}$ perturbed by i.i.d. additive Gaussian
noise, i.e. $y^{ij}_{k}\sim{\mathcal{N}}(f_{ij}(x^{ij}_{k}),\sigma_{i}^{2})$.

For each $i=1,...,N$, we assume all functions in $F_{i}$ are i.i.d. function
samples from the same GP: $\mathcal{GP}_{i}=\mathcal{GP}(\mu_{i},k_{i})$.

$\theta_{i}$$\mu_{i}$$f_{ij}$$k_{i}$$a$$y_{k}^{ij}$$x_{k}^{ij}$$\sigma_{i}$$L_{ij}$$M_{i}$$N$
Figure 1: Graphical model for a hierarchical GP.

For each function set $F_{i}$ and its corresponding $\mathcal{GP}_{i}$ with
mean function $\mu_{i}:\mathcal{X}_{i}\rightarrow\mathbb{R}$ and kernel
$k_{i}:\mathcal{X}_{i}\times\mathcal{X}_{i}\rightarrow\mathbb{R}$, we denote
the parameters of $\mu_{i},k_{i}$ by $\theta_{i}$.

Our major assumption is that the GP parameters and noise standard deviation
$\\{(\theta_{i},\sigma_{i})\\}_{i=1}^{N}$ are i.i.d. samples from a
distribution, as in $(\theta_{i},\sigma_{i})\sim p((\theta,\sigma);a)$, which
is parameterized by $a$. The distribution $p((\theta,\sigma);a)$ is a
universal prior for all search spaces. Fig. 1 illustrates the graphical model
of our hierarchical GP.

For simplicity, we use a length-scale based kernel, e.g. Matérn kernel, and
assume that each dimension of the GP parameters and the noise variance are all
independent. We further assume that all length-scale parameters share the same
prior distribution. Thus, we can estimate the distribution $p(\cdot;a)$ using
the super-dataset and this learned universal prior can be used to guide BO on
new black-box functions, either from seen search spaces or from unseen search
spaces.

#### Evaluation metrics.

For an unseen blackbox “testing function” $f$, we run BO on the function for
$T$ steps and accumulate the observations $\\{x_{t},y_{t}\\}$. We evaluate the
optimization performance of a given method on an individual testing function
using the normalized simple regret

$\displaystyle
r_{T}=\frac{\mathrm{max}_{x\in\mathcal{X}}f(x)-\mathrm{max}_{t\in[T]}y_{t}}{\mathrm{max}_{x\in\mathcal{X}}f(x)-\mathrm{min}_{x\in\mathcal{X}}f(x)}.$
(1)

We use the mean of normalized simple regret on all testing functions as the
overall evaluation metric.
Additionally, we also report the negative log-likelihood (NLL) of the compared
models (either a hierarchical GP or a vanilla GP) on the training and testing
datasets. For a vanilla GP parameterized by $(\theta,\sigma)$, the NLL on a
dataset $D$ is defined as $\mathrm{NLL}(\theta,\sigma)=-\log
p(D|\theta,\sigma)$. For a hierarchical GP parameterized by $a$, the NLL is
defined as follows and is computed via sampling in practice.

$\begin{split}&\mathrm{NLL}(a)=-\log p(D\mid a)\\\
&=-\log\left(\int_{(\theta,\sigma)}p(D\mid\theta,\sigma)p((\theta,\sigma);a)\mathrm{d}(\theta,\sigma)\right)\\\
&\approx-\log\left(\frac{1}{Q}\sum_{q=1}^{Q}p(D\mid\theta_{q},\sigma_{q})p((\theta_{q},\sigma_{q});a)\right)\end{split}$
(2)

where $\\{(\theta_{q},\sigma_{q})\\}_{q=1}^{Q}$ are i.i.d. samples from
$p((\theta,\sigma);a)$.

##  4 Our method

The HyperBO+ framework we propose consists of mainly two phases: (1) Training:
estimate the universal prior $a$ from the super-dataset
$D=\\{D_{i}\\}_{i=1}^{N}$ with a two-step approach. (2) Optimization: running
BO with the hierarchical GP parameterized by the learned $a$ on testing
functions.

###  4.1 Two-step pre-training

#### Estimating GP parameters of each search space.

For each function collection $F_{i}$ with domain $\mathcal{X}_{i}$, we can
infer its GP parameters $\theta_{i}$ and noise standard deviation $\sigma_{i}$
by minimizing the negative log-likelihood of the dataset as in the original
HyperBO (Wang et al.,, 2022):

$\displaystyle\begin{split}&L(\theta,\sigma\mid D_{i})=-\log
p(D_{i}\mid\theta,\sigma)\\\ &=-\sum_{j=1}^{M_{i}}\log
p(D_{ij}\mid\theta,\sigma)\\\
&=\sum_{j=1}^{M_{i}}\left(\frac{1}{2}Y_{ij}^{\intercal}K^{-1}Y_{ij}+\frac{1}{2}\log|K|+\frac{L_{ij}}{2}\log
2\pi\right)\end{split}$ (3)

where
$\bm{x}^{ij}=[x_{k}^{ij}]_{k=1}^{L_{ij}},\bm{y}^{ij}=[y_{k}^{ij}]_{k=1}^{L_{ij}}$,
$Y_{ij}=\left(\bm{y}^{ij}-\mu(\bm{x}^{ij})\right)$,
$K=k(\bm{x}^{ij})+\sigma^{2}\bm{I}$, and $\mu$ and $k$ are the mean and
covariance functions parameterized by $\theta$. Specifically, we minimize the
negative log-likelihood to get the maximum likelihood estimator (MLE) of
$\theta_{i},\sigma_{i}$ as

$\displaystyle(\hat{\theta}_{i},\hat{\sigma}_{i})_{ML}=\operatorname{arg\,min}_{\theta,\sigma}L(\theta,\sigma,D_{i}).$
(4)

#### Estimate the universal prior.

Using the estimated $\\{(\hat{\theta}_{i},\hat{\sigma}_{i})\\}_{i=1}^{N}$ from
all datasets, we can use the the maximum likelihood estimator for the
universal prior parameter $a$ as
$\hat{a}=\operatorname{arg\,max}_{a}p(\\{(\hat{\theta}_{i},\hat{\sigma}_{i})\\}_{i=1}^{N};a)$.
For instance, we can choose the universal prior to be a normal distribution,
i.e.

$p(\\{(\hat{\theta}_{i},\hat{\sigma}_{i})\\}_{i=1}^{N};a)={\mathcal{N}}(a_{0},\textbf{I}a_{1})$

where $a=[a_{0},a_{1}]\in\mathbb{R}^{2}$, and directly estimate the mean and
variance as the parameter $a$. For length-scales, we treat all estimated
length-scales as i.i.d. samples from the same distribution.

Algorithm 1 HyperBO+ pre-training and Bayesian optimization with acquisition
function $ac(\cdot)$.

1:function HyperBO+ ($f,D$)

2: for $D_{i}\in D$ do

3: $\hat{\theta}_{i},\hat{\sigma}_{i}\leftarrow\textsc{Pre-Train}(D_{i})$

4: end for

5:
$\hat{a}\leftarrow\textsc{MLE}(\\{\hat{\theta}_{i},\hat{\sigma}_{i}\\}_{i=1}^{N})$

6: $D_{f}\leftarrow\emptyset$

7: for $t=1,\cdots,T$ do

8:
$x_{t}\leftarrow\underset{x\in\mathfrak{X}}{\operatorname{arg\,max}}{\,ac_{t}\left(x;\hat{a})\right)}$
(Eq. 5)

9: $y_{t}\leftarrow$ Observe$\left(f(x_{t})\right)$

10: $D_{f}\leftarrow D_{f}\cup\\{(x_{t},y_{t})\\}$

11: end for

12: return $D_{f}$

13:end function

###  4.2 Bayesian optimization

HyperBO+ models functions via a hierarchical GP with pre-trained universal
prior parameter $\hat{a}$. At step $t$ of Bayesian optimization when
optimizing a testing function $f$, we compute the acquisition function for
HyperBO+ by first computing the posterior distribution of GP parameters and
noise stand deviation $(\theta,\sigma)$ themselves, then taking the average of
acquisition function values according to $(\theta,\sigma)$-s sampled from this
posterior distribution:

$\displaystyle
ac_{t}(x;\hat{a})=\sum_{r=1}^{R}\left[ac_{t}(x;\theta_{r},\sigma_{r})p((x_{k},y_{k})_{k=1}^{t}|\theta_{r},\sigma_{r})\right]$
(5)

where $(\theta_{1},\sigma_{1}),\dots,(\theta_{R},\sigma_{R})$ are i.i.d.
samples from the prior distribution $p((\theta,\sigma);\hat{a})$. By injecting
the term $p((x_{k},y_{k})_{k=1}^{t}|\theta_{r},\sigma_{r})$ to the acquisition
function value of each sample, we convert sampling from the prior distribution
to sampling from the posterior distribution
$p((\theta,\sigma)(x_{k},y_{k})_{k=1}^{t};\hat{a})$ according to Bayes’ Rule
$p((\theta,\sigma)|(x_{k},y_{k})_{k=1}^{t};\hat{a})=p((x_{k},y_{k})_{k=1}^{t}|\theta_{r},\sigma_{r})p((\theta_{r},\sigma_{r});\hat{a})/p((x_{k},y_{k})_{k=1}^{t})$.
Here we can ignore the denominator $p((x_{k},y_{k})_{k=1}^{t})$ since the BO
selection is invariant under any positive affine transformation of the
acquisition function. Note that posterior sampling could be a preferred option
than re-weighting with likelihood. We find our re-weighting strategy to be
sufficient since the samples from the pre-trained prior are often
representative of the samples from the posterior.

##  5 Theoretical analysis

Figure 2: Converting a series of sub-datasets equivalently to samples from a
single function with increasing observation locations.

Showing the asymptotic behaviors of the estimated GP parameter
$\hat{\theta}_{i}$ as the number of datasets and sub-datasets approaches
infinity is non-trivial. These asymptotics differs from standard analyses of
estimators where samples are drawn i.i.d. from the same distribution.
Oftentimes, the determinant $|K|$ and gram matrix inverse $K^{-1}$ in the
likelihood function has no explicit expressions (Bachoc,, 2021) thus making it
even more challenging to analyze the properties of the estimated GP parameter
based on type II maximum likelihood.

We recall that in the increasing domain setting, there is a minimal spacing
between any given pairs of observations and thus make the sampling domain
unbounded as the number of observations increases. That is, there exists a
fixed $\Delta>0$ such that

$\displaystyle\text{inf}_{i,j\in\mathbb{N},i\neq j}\|x_{i}-x_{j}\|\geq\Delta.$
(6)

In previous work showing the consistency of ML estimators for GP parameters,
they consider a a single function sampled from a zero-mean GP with $n$
observations to minimize the NLL

$\displaystyle
L_{n}(\theta)=\frac{1}{n}\log(|K_{\theta}|)+\frac{1}{n}y^{T}K_{\theta}^{-1}y$

for
$\hat{\theta}_{ML}\in\operatorname{arg\,min}_{\theta\in\Theta}L_{n}(\theta)$,
where $K_{\theta}$ is the covariance matrix parameterized by $\theta$. The
authors also showed that though the MLE given finite observations may not be
unique, it converges to the unique global minimizer with probability
converging to one as the number of observations goes to infinity. Extensions
of such results on increasing domain for GP with non-zero mean function were
discussed in Bachoc et al., (2020).

We now analyze the asymptotic behavior of the first step in §4.1 \- GP
parameter estimation for each search space. We make the following assumptions:

  1. (1)

Each GP is zero-mean with known noise variance $\sigma^{2}$. Its ground-truth
kernel belongs to the Matérn covariance family with known smoothness term
$\nu$ (e.g. $\nu=3/2$) and the parameter $\theta$ consists of the length-scale
only.

  2. (2)

The kernel for each GP is fixed and can be parameterized by a unique
$\theta_{0}\in\Theta$ for some compact $\Theta$. Moreover, the setup is well-
specified, where the true covariance parameters $\theta_{0}$ do belong to the
search space $\Theta$ used for estimation.

In our problem setup, a dataset $D_{i}=\\{D_{ij}\\}_{j=1}^{M_{i}}$
($i=1,\dots,N$) contains a finite number of observations in each of its sub-
datasets. For the following analysis we assume that there exists a fixed
$\delta>0$ such that $\|x^{ij}_{k}-x^{ij}_{k^{\prime}}\|\geq\delta,(k\neq
k^{\prime})$ for all sub-datasets $D_{ij}$. As the number of sub-datasets
$M_{i}\rightarrow\infty$, we can show that it is equivalent to the increasing
domain setting within a single function from GP and prior work guarantee good
asymptotic properties of the ML estimation for covariance parameters (Mardia
and Marshall,, 1984; Stein,, 1999).

######  Lemma 1.

Given assumption (1), for any $i\in\\{1,\cdots,N\\}$, $\mathcal{GP}_{i}$ and
its corresponding _dataset_ $D_{i}=\\{D_{i1},\cdots,D_{iM_{i}}\\}$, there
exists a sub-dataset $\bar{D}$ observed on a function
$f^{\prime}\sim\mathcal{GP}_{i}$, such that $\bar{D}$ satisfies the increasing
domain characterization in Eq. 6, and

$\displaystyle L(\theta\mid D_{i})\equiv L(\theta\mid\\{\bar{D}\\}),$ (7)

where $L(\theta\mid\mathcal{D})=-\sum_{\tilde{D}\in\mathcal{D}}\log
p(\tilde{D}\mid\theta).$

Intuitively, this lemma shows that observations from multiple independently
generated functions from a fixed GP can be viewed as being sampled from one
function, with infinitely large interval between sub-datasets’ observations.
We illustrate this process with two sub-datasets on a 1-dimensional domain in
Fig. 2. A detailed proof is presented in §B.

We introduce our assumption (3) that there is sufficient information in
sampling locations $x_{1},x_{2},\dots,x_{n}$ of the augmented sub-dataset
$\bar{D}$ to distinguish covariance functions $k_{\theta}$ with
$k_{\theta_{0}}$ (Bachoc,, 2014). That is, for $\theta\in\Theta$

$\displaystyle\leavevmode\resizebox{436.8687pt}{}{$\liminf_{n\rightarrow\infty}\inf_{\|\theta-\theta_{0}\|\geq\epsilon}\frac{1}{n}\sum_{i\neq
j}(k_{\theta}(x_{i}-x_{j}))-k_{\theta_{0}}(x_{i}-x_{j}))^{2}>0$}.$ (Asymptotic
Identifiability)

######  Theorem 2.

Given assumptions (1)-(3), for dataset $D_{i}$ with $M_{i}$ sub-datasets
generated from the same mean and covariance function, as
$M_{i}\rightarrow\infty$, we have
$\hat{\theta}_{ML}\rightarrow^{p}\theta_{0}.$

######  Theorem 3.

Given assumptions (1)-(3), as the number of datasets, and sub-datasets
$N,M_{i}\rightarrow\infty,\forall i$, the maximum likelihood estimator for the
prior distribution $a$ is consistent.

Theorem 2 and 3 are easy to show and the proofs can be found in §B.

######  Remark.

(Ye and Chen, (2017) and Rice, (2006))
(1). When the prior is assumed to be from a Gamma distribution
$\Gamma(\alpha,\beta)$, the MLE $\hat{\alpha},\hat{\beta}$ is consistent and
asymptotically normal as the number of samples $n\rightarrow\infty$ with
distribution
$\mathcal{N}((\alpha,\beta),\frac{1}{n}\mathcal{I}(\alpha,\beta)^{-1})$, where
$\mathcal{I}$ is the Fisher information matrix.
(2). Similar results hold when the prior is a normal distribution parametrized
by $\mu,\sigma$: as $n\rightarrow\infty$, we have
$\hat{\mu}_{ML},\hat{\sigma}_{ML}\rightarrow\mu,\sigma$.

A key observation from our theoretical analyses is that with sufficient
observations in each dataset $D_{i}$, increasing the number of sub-datasets
can effectively improve the estimation of the covariance parameters and thus
the prior. Fig. 3 demonstrates the asymptotic behavior empirically for fitting
the length-scale parameter of a single GP using an increasing number of sub-
datasets. Note that while our analysis considers only the length-scale
estimation, empirical results show strong evidence for the consistency of
function mean and signal/noise variance as well, which is illustrated in Fig.
11 in the supplement. Empirical asymptotic behavior of the two-step pre-
training method is demonstrated in §6.2.

Figure 3: Estimated length-scale parameter (distribution plotted over 50
random runs) of a 1-dimensional GP as the number of sub-datasets increases.
Each sub-dataset has 25 observations. The variance of estimated length-scale
parameter decreases as the number of sub-datasets increases.

##  6 Experiments

We demonstrate the performance of HyperBO+ on two datasets: a synthetic Super-
dataset and HPO-B (Pineda-Arango et al.,, 2021), a collection of real-world
hyperparameter tuning tasks that involves multiple search spaces. We generate
the synthetic super-dataset with a known ground-truth prior distribution over
GP parameters, which is helpful for comparing with the ground-truth and
empirically demonstrating asymptotic behaviors of the two-step pre-training
method. The same cannot be done for real-world datasets since there is not a
known ground-truth prior.

Synthetic Super-dataset is generated over multiple search spaces following the
generative process illustrated in Fig. 1. The super-dataset includes 20
datasets (each with its own search space) with 10 sub-datasets in each
dataset. Each sub-dataset includes noisy observations at 300 random input
locations in its respective search space. The dimensions of search spaces are
between 2 and 5.

HPO-B Super-dataset (Pineda-Arango et al.,, 2021) is a large-scale benchmark
for hyperparameter optimization. As a super-dataset, it consists of 16
different search spaces and more than 6 million evaluations in total. The
dimensions of these search spaces vary from 2 to 18.

###  6.1 Experiment setups

To test the capability of HyperBO+ to generalize to new tasks in both seen and
unseen search spaces, we design experiments for the two super-datasets and
compare HyperBO+ (and its variants) with 5 competitive baselines including
HyperBO (Wang et al.,, 2022), which is known to be the state-of-the-art prior
pre-training method for BO.

For the Synthetic Super-dataset, we design the following two setups as two
different ways of splitting training and testing data.

Setup A is designed to demonstrate the ability of HyperBO+ to generalize its
learned prior to unseen search spaces, potentially with different input
dimensions. We split the Synthetic Super-dataset into training datasets and
testing datasets. In particular, we use 80% of the super-dataset for training
and the remaining 20% for testing. To test BO performances of different
methods, we report their average normalized simple regrets on sub-datasets in
the testing datasets. We also report the NLLs of different GP-based methods on
training datasets and testing datasets.

Setup B aims to to test the ability of HyperBO+ to generalize to new functions
in seen search spaces so that we can compare its performance with HyperBO. We
split each dataset in the Synthetic Super-dataset into training sub-datasets
and testing sub-datasets. We use 80% of each dataset as the training set of
sub-datasets and the remaining 20% for testing. We report the average
normalized simple regrets on all testing sub-datasets and the NLLs of
different GP-based methods on training and testing sub-datasets. HyperBO
individually pre-trains one GP per search space, while HyperBO+ learns a
universal hierarchical GP prior by pre-training on the collection of training
sub-datasets of all search spaces and uses the same pre-trained prior to run
BO on testing sub-datasets for all search spaces.

HPO-B Setup HPO-B Super-dataset comes with a pre-specified train/test split:
each dataset (search space) contains multiple training sub-datasets and
multiple testing sub-datasets. For BO performance, we report the average
normalized simple regrets on testing sub-datasets in all search spaces. We
also report the NLLs of different GP-based methods on training sub-datasets
and testing sub-datasets. HyperBO works by individually pre-training a single
GP per search space. On the other hand, HyperBO+ learns a universal
hierarchical GP prior by pre-training on the collection of all training sub-
datasets.

Baselines: (1) Random sampling for optimization. (2) A hand-specified (and
potentially misspecified) hierarchical GP prior, fixed over all search spaces.
(3) A non-informative hierarchical GP prior, fixed over all search spaces. (4)
HyperBO (Wang et al.,, 2022). Note that HyperBO is only applicable to Setup B.
(5) The ground-truth hierarchical GP prior, which is available to Synthetic
Super-dataset. Please see §C for the exact configurations, including GP prior
parameters, of all baselines.

Variants of HyperBO+: In addition to the standard HyperBO+, we also test the
performance of the following two variants. (1) Instead of fitting a
parameterized universal prior with MLE, xHyperBO+ uses the uniform
distribution over the finite set of GP parameter values
$\\{(\hat{\theta}_{i}),\hat{\sigma_{i}}\\}_{i=1}^{N}$ as an explicit prior
distribution (i.e. crossing out the MLE step at Line 5 of Alg. 1). (2) For the
HPO-B Super-dataset, instead of pre-training HyperBO+ on training sub-datasets
of all search spaces, zHyperBO+ is pre-trained with the training sub-datasets
that do not belong to the search space of the testing sub-dataset (i.e. zero
knowledge of its own search space). zHyperBO+ is the same as the standard
HyperBO+ in Setup A of the Synthetic Super-dataset.

We fit the GP parameters in each search space by minimizing Eq. 3 using L-BFGS
(Liu and Nocedal,, 1989) for the Synthetic Super-dataset, and using the Adam
optimizer (Kingma and Ba,, 2015; Wistuba and Grabocka,, 2021) for HPO-B. All
the GP-based methods use a constant mean function and a Matérn kernel with
smoothness parameter $\nu=3/2$. More experiment details and analyses can be
found in §C and §D.

###  6.2 Empirical asymptotic behaviors

Figure 4: Estimated shape parameter $\hat{\alpha}$ (y-axis) of Gamma
distribution prior for the length-scale GP parameter w.r.t. the number of
training datasets (x-axis). We show the mean and violin plots for
$\hat{\alpha}$ over 5 random seeds. Figure 5: We plot the probability density
functions of the Gamma distributions that model length-scales. The pre-trained
Gamma distributions with 16 training datasets are more stable than those with
2 training datasets and match well with the ground-truth. Figure 6: NLLs of
the testing datasets on pre-trained priors as the number of training datasets
increases. We show the mean and violin plots for the NLLs over 5 random seeds.

To empirically demonstrate asymptotic properties of our two-step pre-training
method analyzed in §5, we present results with increasing number of training
datasets used in Setup A on Synthetic Super-dataset. The testing super-dataset
is fixed, while the training datasets are uniformly randomly sampled from all
available training datasets. We repeat the experiments 5 times with different
random seeds.

Fig. 4 shows the estimated shape parameter $\alpha$ of the Gamma distribution
prior $\Gamma(\alpha,\beta)$ which models the length-scale parameter of the GP
for each training dataset. We illustrate how well the pre-trained Gamma prior
matches the ground-truth by plotting the probability density functions in Fig.
5.

We observe that as the number of training datasets increases, the variance of
the estimated $\hat{\alpha}$ gradually decreases and the mean becomes closer
to the ground-truth in Fig. 4. Similarly, Fig. 5 also shows that more training
datasets help the stability of pre-training. This is consistent with our
theoretical analysis of asymptotic properties of the two-step pre-training
method. The mean value of $\alpha$ is still a bit off from the ground-truth,
which is expected given that there are only 16 samples for the second step of
pre-training. However, this estimation quality is good enough for HyperBO+ to
outperform baselines on this Synthetic Super-dataset as shown below in §6.3.

In Fig. 6, we show the NLL of the testing datasets on the estimated
hierarchical GP prior as a function of the number of training datasets. Both
the mean and variance of the NLL on testing dataset drop as the number of
training datasets increases, which is consistent with our asymptotic analysis.

###  6.3 Bayesian optimization results

Figure 7: Average normalized simple regret of each method during BO in both
Setup A and Setup B for the Synthetic Super-dataset. The averages are taken
over 5 random seeds, and the highlighted areas show
$\mathrm{mean}\pm\mathrm{std}$ for each method. Results for acquisition
functions PI and EI are shown here. Figure 8: Average normalized simple regret
of each method during BO for the HPO-B Super-dataset. The averages are taken
over 5 random seeds, and the highlighted areas show
$\mathrm{mean}\pm\mathrm{std}$ for each method. Results for acquisition
functions PI and EI are shown here. Figure 9: NLL evaluations of
training/testing data on HyperBO+ and HyperBO for the HPO-B Super-dataset.
Figure 10: The estimated length-scale values for each dimension of every
search space and the learned Gamma prior distribution of length-scale for
HPO-B Super-dataset. Table 1: NLL evaluations of the training and testing data
on different models and normalized simple regret of running BO with each
method (PI acquisition function) in both Setup A and Setup B for the Synthetic
Super-dataset. ($\mathrm{mean}\pm\mathrm{std}$)

Methods Setup A Setup B Train NLL Test NLL BO Regret Train NLL Test NLL BO
Regret Non-informative $122.2\pm 0.4$ $106.7\pm 0.2$ $0.192\pm 0.01$ $128.8\pm
0.3$ $127.0\pm 0.7$ $0.180\pm 0.01$ Hand-specified $99.9\pm 0.4$ $52.2\pm 0.6$
$0.228\pm 0.03$ $90.5\pm 0.6$ $91.5\pm 1.0$ $0.207\pm 0.02$ HyperBO N/A N/A
N/A $48.4\pm 0.4$ $50.1\pm 0.9$ $0.014\pm 0.00$ Ground-truth $62.2\pm 0.5$
$22.4\pm 1.1$ $0.022\pm 0.01$ $53.9\pm 0.4$ $55.5\pm 1.0$ $0.014\pm 0.00$
HyperBO+ $62.2\pm 0.5$ $22.8\pm 1.0$ $0.032\pm 0.01$ $54.2\pm 0.4$ $56.1\pm
1.0$ $0.015\pm 0.00$ xHyperBO+ $61.9\pm 0.2$ $22.2\pm 0.7$ $0.030\pm 0.00$
$53.4\pm 0.3$ $55.1\pm 0.6$ $0.013\pm 0.00$

Table 2: NLL evaluations of the training and testing data on different models
and normalized simple regret of running BO with each method (PI acquisition
function) for the HPO-B Super-dataset. ($\mathrm{mean}\pm\mathrm{std}$)

Methods Train NLL Test NLL BO Regret Non-informative $90.5\pm 0.0$ $96.7\pm
0.4$ $0.088\pm 0.01$ Hand-specified $-23.9\pm 0.0$ $-26.3\pm 0.1$ $0.091\pm
0.01$ HyperBO $-197.1\pm 0.8$ $-180.7\pm 2.5$ $0.044\pm 0.01$ HyperBO+
$-261.4\pm 0.4$ $-272.9\pm 1.9$ $0.038\pm 0.00$ zHyperBO+ $-256.0\pm 0.7$
$-268.3\pm 2.1$ $0.041\pm 0.00$ xHyperBO+ $-263.2\pm 0.5$ $-273.5\pm 2.1$
$0.041\pm 0.01$

Fig. 7 shows the BO performances of compared methods in both Setup A and Setup
B for the Synthetic Super-dataset, and Fig. 8 shows the results for the HPO-B
Super-dataset. The budget for BO is 50, and there are 5 initial observations
that are randomly sampled for each of the 5 random seeds. The acquisition
function used for all GP-based methods is Probability of Improvement (PI)
(Kushner,, 1964) with target value $\max(y_{t})+0.1$. As shown by Wang and
Jegelka, (2017), PI can obtain high BO performance by setting good target
values. Results with other types of acquisition functions are included in the
supplement.

In Setup A for the Synthetic Super-dataset (Fig. 7(a)), HyperBO+ outperforms
the non-informative hierarchical GP, the hand-specified hierarchical GP, and
the random baseline, demonstrating its ability to learn a universal prior from
training datasets and generalizing to tasks in new search spaces that are
unseen during pre-training. HyperBO+ achieves good BO performance similar to
the ground-truth model. This shows that the two-step pre-training method is
very effective in learning the prior parameters.

In Setup B for the Synthetic Super-dataset (Fig. 7(b)), HyperBO+ again
achieves lower regrets than the non-informative hierarchical GP, the hand-
specified hierarchical GP, and the random baseline. HyperBO+ and HyperBO
achieve BO performances similar to using the ground-truth prior. xHyperBO+
achieves a slightly lower regret than the standard HyperBO+ in both setups.

For the HPO-B Super-dataset (Fig. 8), HyperBO+ achieves lower regrets than the
non-informative hierarchical GP, the hand-specified hierarchical GP, and the
random baseline. Notably, HyperBO+ outperforms HyperBO on HPO-B, which is
rather surprising given that HyperBO specifically pre-trains a GP that targets
a single search space. In contrast, HyperBO+ learns one universal prior on
training sub-datasets from all search spaces, and the learned universal prior
would have to automatically capture the shared characteristics of testing
functions in different search spaces.

We conjecture that the performance difference is because HyperBO overfits or
fits an incorrect GP for some search spaces which damages its performance on
the testing sub-datasets from those search space. Meanwhile, HyperBO+ uses a
single learned prior for BO in all search spaces, which to some extent
regularizes an incorrect GP learned on a particular search space with GPs
learned from other search spaces. Fig. 9 shows that HyperBO can obtain high
test NLLs despite having much lower train NLLs in some search spaces of HPO-B.

Notably, zHyperBO+ outperforms the baselines with only slightly worse
performance than the standard HyperBO+ even though it is not trained on the
training sub-datasets of the search space it is tested in. The hyperparameters
in tuning tasks of HPO-B can be entirely irrelevant, e.g. hyperparameters for
SVM and Random Forest. As a result, the training data for zHyperBO+ may not
even be equipped to inform new tasks for testing. The good performance of
zHyperBO+ further demonstrates the capability of HyperBO+ to generalize across
search spaces in real-world problems.

Interestingly, xHyperBO+ achieves a similar and slightly worse performance
than the standard HyperBO+. Recall that xHyperBO+ does not use the MLE step to
estimate parameterized priors, and instead, just use the uniform distribution
over the parameters of mean and kernel functions. Fig. 10 shows the histogram
of the length-scale values learned from the pre-training phase, and the
probability density of the universal Gamma prior on length-scales learned by
HyperBO+ in Line 5 of Alg. 1. Despite the fact that our universal Gamma prior
assigns low probability at length-scale=$9$, which is inconsistent with the
histogram, using the pre-trained Gamma prior as part of the standard HyperBO+
still performs very well or even slightly better than xHyperBO+. The results
here also demonstrate that some training data can have negative impacts on
testing tasks, but HyperBO+ is able to mitigate the negative transfer learning
issues.

###  6.4 Evaluations of NLLs

Tab. 1 shows the NLL evaluations of training/testing data on the compared
models together with the BO regret in both Setup A and Setup B for the
Synthetic Super-dataset, while Tab. 1 shows the results for the HPO-B super-
dataset.

The results show that HyperBO+ achieves a much lower NLL than the hand-
specified hierarchical GP on both super-datasets, and achieves a similar NLL
as the ground-truth hierarchical GP in both setups on the Synthetic Super-
dataset. This further demonstrates the effectiveness of the two-step pre-
training of HyperBO+. The two variants of HyperBO+, zHyperBO+ and xHyperBO+,
both achieve a similar NLL to the standard HyperBO+.

In Setup B, the original HyperBO gets a lower NLL than HyperBO+ and the
ground-truth hierarchical GP on the Synthetic Super-dataset, which makes sense
as HyperBO specifically pre-trains a GP that targets a single search space.
However, HyperBO has higher NLLs than HyperBO+ on the HPO-B super-dataset.
This is in line with our observation of BO performances on HPO-B. The reason
for this could to be similar to what we conjectured in §6.3 \- that is,
HyperBO overfits in some of the search spaces of HPO-B, while HyperBO+ is
capable of regularizing incorrect GP parameters learned in these parameters
with GPs learned in other search spaces because it uses the same universal
prior for all search spaces.

##  7 Conclusions

We present HyperBO+, the first GP-based transfer learning BO method that works
on universal search spaces. The key idea is to pre-train a universal
hierarchical GP prior on training data partitioned over search spaces and
functions within each search space. We demonstrate the appealing asymptotic
properties of our pre-training method and show HyperBO+ outperforms
competitive baselines through experiments on challenging hyperparameter tuning
tasks.

### Acknowledgments

We thank Jasper Snoek and Eytan Bakshy for helpful conversations and/or
feedback. Our work also benefited from Microsoft Azure credits provided by the
Harvard Data Science Initiative, as well as Google Cloud Platform Credit Award
provided by Google.

## References

  * Bachoc, (2014) Bachoc, F. (2014).  Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes.  Journal of Multivariate Analysis, 125:1–35.
  * Bachoc, (2021) Bachoc, F. (2021).  Asymptotic analysis of maximum likelihood estimation of covariance parameters for Gaussian processes: an introduction with proofs.  In Advances in Contemporary Statistics and Econometrics, pages 283–303. Springer.
  * Bachoc et al., (2020) Bachoc, F., Betancourt, J., Furrer, R., and Klein, T. (2020).  Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes.  Electronic Journal of Statistics, 14(1):1962–2008.
  * Balandat et al., (2020) Balandat, M., Karrer, B., Jiang, D. R., Daulton, S., Letham, B., Wilson, A. G., and Bakshy, E. (2020).  BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization.  In Advances in Neural Information Processing Systems (NeurIPS).
  * Baxter, (1996) Baxter, J. (1996).  A Bayesian/information theoretic model of bias learning.  In Conference on Learning Theory (COLT), New York, New York, USA.
  * Bevilacqua et al., (2019) Bevilacqua, M., Faouzi, T., Furrer, R., and Porcu, E. (2019).  Estimation and prediction using generalized wendland covariance functions under fixed domain asymptotics.  The Annals of Statistics, 47(2):828–856.
  * Bogunovic and Krause, (2021) Bogunovic, I. and Krause, A. (2021).  Misspecified Gaussian process bandit optimization.  Advances in Neural Information Processing Systems (NeurIPS).
  * Chen et al., (2022) Chen, Y., Song, X., Lee, C., Wang, Z., Zhang, Q., Dohan, D., Kawakami, K., Kochanski, G., Doucet, A., Ranzato, M., et al. (2022).  Towards learning universal hyperparameter optimizers with transformers.  In Advances in Neural Information Processing Systems (NeurIPS).
  * Cowen-Rivers et al., (2022) Cowen-Rivers, A. I., Lyu, W., Tutunov, R., Wang, Z., Grosnit, A., Griffiths, R. R., Maraval, A. M., Jianye, H., Wang, J., Peters, J., et al. (2022).  Hebo: Pushing the limits of sample-efficient hyper-parameter optimisation.  Journal or Artificial Intelligence Research (JAIR), 74:1269–1349.
  * Drieß et al., (2017) Drieß, D., Englert, P., and Toussaint, M. (2017).  Constrained Bayesian optimization of combined interaction force/task space controllers for manipulations.  In International Conference on Robotics and Automation (ICRA), pages 902–907. IEEE.
  * Eggensperger et al., (2021) Eggensperger, K., Müller, P., Mallik, N., Feurer, M., Sass, R., Klein, A., Awad, N., Lindauer, M., and Hutter, F. (2021).  HPOBench: A collection of reproducible multi-fidelity benchmark problems for HPO.  In Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2).
  * Golovin et al., (2017) Golovin, D., Solnik, B., Moitra, S., Kochanski, G., Karro, J. E., and Sculley, D. (2017).  Google Vizier: A service for black-box optimization.  In ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD).
  * Griffiths and Hernández-Lobato, (2020) Griffiths, R.-R. and Hernández-Lobato, J. M. (2020).  Constrained Bayesian optimization for automatic chemical design using variational autoencoders.  Chemical Science, 11(2):577–586.
  * Karvonen and Oates, (2022) Karvonen, T. and Oates, C. J. (2022).  Maximum likelihood estimation in Gaussian process regression is ill-posed.  arXiv preprint arXiv:2203.09179.
  * Kaufman and Shaby, (2013) Kaufman, C. and Shaby, B. A. (2013).  The role of the range parameter for estimation and prediction in geostatistics.  Biometrika, 100(2):473–484.
  * Kim et al., (2019) Kim, B., Wang, Z., Kaelbling, L. P., and Lozano-Pérez, T. (2019).  Learning to guide task and motion planning using score-space representation.  International Journal of Robotics Research (IJRR), 38(7):793–812.
  * Kingma and Ba, (2015) Kingma, D. P. and Ba, J. (2015).  Adam: A method for stochastic optimization.  In Bengio, Y. and LeCun, Y., editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings.
  * Kushner, (1964) Kushner, H. J. (1964).  A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise.  Journal of Fluids Engineering, 86(1):97–106.
  * Lam et al., (2018) Lam, R., Poloczek, M., Frazier, P., and Willcox, K. E. (2018).  Advances in Bayesian optimization with applications in aerospace engineering.  In AIAA Non-Deterministic Approaches Conference, page 1656.
  * Liu and Nocedal, (1989) Liu, D. C. and Nocedal, J. (1989).  On the limited memory bfgs method for large scale optimization.  Math. Program., 45(1–3):503–528.
  * Mardia and Marshall, (1984) Mardia, K. V. and Marshall, R. J. (1984).  Maximum likelihood estimation of models for residual covariance in spatial regression.  Biometrika, 71(1):135–146.
  * Minka and Picard, (1997) Minka, T. P. and Picard, R. W. (1997).  Learning how to learn is learning with point sets.  Technical report, MIT Media Lab.
  * Perrone et al., (2018) Perrone, V., Jenatton, R., Seeger, M., and Archambeau, C. (2018).  Scalable hyperparameter transfer learning.  In Advances in Neural Information Processing Systems (NeurIPS), pages 6846–6856.
  * Pineda-Arango et al., (2021) Pineda-Arango, S., Jomaa, H. S., Wistuba, M., and Grabocka, J. (2021).  HPO-B: A large-scale reproducible benchmark for black-box HPO based on openml.  Neural Information Processing Systems (NeurIPS) Track on Datasets and Benchmarks.
  * Pyzer-Knapp, (2018) Pyzer-Knapp, E. O. (2018).  Bayesian optimization for accelerated drug discovery.  IBM Journal of Research and Development, 62(6):2–1.
  * Rice, (2006) Rice, J. A. (2006).  Mathematical Statistics and Data Analysis.  Cengage Learning.
  * Schulz et al., (2016) Schulz, E., Speekenbrink, M., Hernández-Lobato, J. M., Ghahramani, Z., and Gershman, S. J. (2016).  Quantifying mismatch in Bayesian optimization.  In NeurIPS Workshop on Bayesian Optimization: Black-box Optimization and Beyond.
  * Shields et al., (2021) Shields, B. J., Stevens, J., Li, J., Parasram, M., Damani, F., Alvarado, J. I. M., Janey, J. M., Adams, R. P., and Doyle, A. G. (2021).  Bayesian reaction optimization as a tool for chemical synthesis.  Nature, 590(7844):89–96.
  * Snoek et al., (2012) Snoek, J., Larochelle, H., and Adams, R. P. (2012).  Practical Bayesian optimization of machine learning algorithms.  In Advances in Neural Information Processing Systems (NeurIPS).
  * Stein, (1999) Stein, M. L. (1999).  Interpolation of Spatial Data: Some Theory for Kriging.  Springer Science & Business Media.
  * Swersky et al., (2013) Swersky, K., Snoek, J., and Adams, R. P. (2013).  Multi-task Bayesian optimization.  In Advances in Neural Information Processing Systems (NeurIPS).
  * Volpp et al., (2020) Volpp, M., Fröhlich, L. P., Fischer, K., Doerr, A., Falkner, S., Hutter, F., and Daniel, C. (2020).  Meta-learning acquisition functions for transfer learning in Bayesian optimization.  In International Conference on Learning Representations (ICLR).
  * Wackerly et al., (2014) Wackerly, D., Mendenhall, W., and Scheaffer, R. L. (2014).  Mathematical Statistics with Applications.  Cengage Learning.
  * Wang et al., (2022) Wang, Z., Dahl, G. E., Swersky, K., Lee, C., Mariet, Z., Nado, Z., Gilmer, J., Snoek, J., and Ghahramani, Z. (2022).  Pre-trained Gaussian processes for Bayesian optimization.  arXiv preprint arXiv:2109.08215.
  * Wang and Jegelka, (2017) Wang, Z. and Jegelka, S. (2017).  Max-value entropy search for efficient Bayesian optimization.  In International Conference on Machine Learning (ICML).
  * Wang et al., (2017) Wang, Z., Jegelka, S., Kaelbling, L. P., and Lozano-Pérez, T. (2017).  Focused model-learning and planning for non-Gaussian continuous state-action systems.  In International Conference on Robotics and Automation (ICRA).
  * Wang et al., (2018) Wang, Z., Kim, B., and Kaelbling, L. P. (2018).  Regret bounds for meta Bayesian optimization with an unknown Gaussian process prior.  In Advances in Neural Information Processing Systems (NeurIPS).
  * Wilson et al., (2015) Wilson, A. G., Dann, C., Lucas, C., and Xing, E. P. (2015).  The human kernel.
  * Wistuba and Grabocka, (2021) Wistuba, M. and Grabocka, J. (2021).  Few-shot Bayesian optimization with deep kernel surrogates.  In International Conference on Learning Representations (ICLR).
  * Ye and Chen, (2017) Ye, Z.-S. and Chen, N. (2017).  Closed-form estimators for the gamma distribution derived from likelihood equations.  The American Statistician, 71(2):177–181.
  * Zhang and Zimmerman, (2005) Zhang, H. and Zimmerman, D. L. (2005).  Towards reconciling two asymptotic frameworks in spatial statistics.  Biometrika, 92(4):921–936.

##  Appendix A Discussions

#### Limitations and future work.

HyperBO+ learns a prior distribution on Gaussian processes (GPs) from existing
data to improve the performance of Bayesian optimization (BO) on new task.
However, there are other components that can also be meta-learned, such as
acquisition functions (Volpp et al.,, 2020), to maximize the effectiveness of
the BO system. One direction of future work is jointly pre-train all
components of BO to allow more flexibility and further improve the
performance. In addition, HyperBO+ currently assumes a stationary kernel
function and constant mean function when fitting the GP parameters. Possible
future work includes relaxing this assumption and incorporating architecture
search for the kernel function and mean function to enrich the space of
hierarchical GP priors.

#### Societal impact of this work.

Hyperparameter tuning for machine learning (ML) models, especially deep
learning models, can be very costly if we repeatedly evaluate a large number
of hyperparameters. Each single evaluation of a hyperparameter value requires
training and testing a new instantiation of the model. Our framework HyperBO+,
with its superior BO performance discussed in §6.3, can help to reduce the
number of evaluations needed for hyperparameter tuning tasks and thus reduce
their computational cost and carbon footprint potentially by a large margin.

##  Appendix B More details of theoretical analysis

#### Proof of Lemma 1

###### Proof.

Consider functions $f_{ij}$ whose observations
$D_{ij}=\\{(x_{k}^{ij},\,y_{k}^{ij})\\}_{k=1}^{L_{ij}}$ is in dataset $D_{i}$,
where $j\in\\{1,\dots,M_{i}\\}$.

Since the domain $\mathcal{X}_{i}$ of $f_{ij}$ is compact, and we define
$Q^{\prime}=\sup_{x,x^{\prime}\in\mathcal{X}_{i}}\|x^{\prime}-x\|+1\ll\infty$.
Since the stationary covariance function only concerns the relative location
between two points, we can shift all inputs jointly while the NLL
$L(\theta\mid\\{D_{ij}\\})$ stays the same. More formally, let
$\bar{D}_{j}=\\{(x_{k}^{ij}+1_{d_{i}}(Q+Q^{\prime})j,\,y_{k}^{ij})\\}_{k=1}^{L_{ij}}$
for some $Q>0$, where $d_{i}$ is the input dimension and $1_{d_{i}}$ is a
$d_{i}-$dimensional vector filled with ones. Because the distance between any
pairs of inputs stays the same, we have
$L(\theta\mid\\{D_{ij}\\})=L(\theta\mid\\{\bar{D}_{j}\\})=-\log
p(D_{ij}\mid\theta)$.

We define an augmented sub-dataset
$\bar{D}=\bar{D}_{1}\cup\cdots\cup\bar{D}_{M_{i}}$.

As $Q\rightarrow\infty$, we can show that the sum of NLLs over a set of sub-
datasets and the NLL of the single augmented sub-dataset are equivalent.
Without loss of generality, consider the Matérn covariance function with
smoothness term $\nu=1/2$ and variance $\sigma^{2}$, length-scale $\rho$:
$C_{1/2}(d)=\sigma^{2}\exp(-\frac{d}{\rho})$. Since the covariance
exponentially decays with the distance $d$ between points, for any
$\epsilon>0$, we can find a $Q$ such that $C(d)<\epsilon$ for $d>Q$.

For the augmented sub-dataset
$\bar{D}=\bar{D}_{1}\cup\cdots\cup\bar{D}_{M_{i}}$, we have
$\|x^{\prime}-x\|>Q$ for any $x\in\bar{D}_{j}$ and
$x^{\prime}\in\bar{D}_{j^{\prime}}$, $j\neq j^{\prime}$. As
$\epsilon\rightarrow 0$, the covariance matrix for this augmented sub-dataset
becomes block diagonal. As a result, we can express its NLL as a sum of NLL
for $D_{i1},\dots,D_{iM_{i}}$ respectively, i.e.
$L(\theta\mid\\{\bar{D}\\})=-\sum_{j=1}^{M_{i}}\log p(\bar{D}_{j}\mid\theta)$.
This gives us the equivalence in Eq. 7.

Note that the distance between inputs in $\bar{D}$ will always be at least
$\min_{j}\min_{k\neq k^{\prime}}\|x^{ij}_{k}-x^{ij}_{k^{\prime}}\|>\delta$,
since any two datapoints in each sub-dataset $D_{ij}$ are at least $\delta$
apart by construction in problem setup. As $M_{i}\rightarrow\infty$, the
augmented sub-dataset $\bar{D}$ is unbounded and satisfies the increasing
domain characterization in Eq. 6. ∎

Another interesting related work is from Karvonen and Oates, (2022), where the
authors pointed out that the MLE for the length-scale parameter is indeed ill-
posed when the data is noiseless and can be unstable to small perturbations.
However, they found that regularization with small additive Gaussian noise
does guarantee well-posedness and this is equivalent to data corrupted by
additive Gaussian noise. This is identical to our setting with the synthetic
data and we found that in practice an accurate estimation for noise allows
reliable inference for covariance parameters.

#### Proof of Theorem 2

###### Proof.

By Lemma 1, our equivalence setup satisfies the conditions in Theorem 1 from
Bachoc, (2021). Therefore, the conclusion holds that
$\hat{\theta}_{ML}\rightarrow^{p}\theta_{0}$. ∎

#### Proof of Theorem 3

###### Proof.

The asymptotic properties for the second training step follows from standard
statistics literature (Wackerly et al.,, 2014) since we assume each datasets’
covariance parameters are drawn i.i.d. from some prior distribution. ∎

##  Appendix C More details of the experiment setups

Our code for the experiments is built upon the codebase of HyperBO (Wang et
al.,, 2022), and is available at https://github.com/Evensgn/hyperbo.

The Synthetic Super-dataset includes 20 datasets (search spaces) with 10 sub-
datasets in each dataset. Each sub-dataset includes noisy observations at 300
input locations in its respective search space. The dimensions of the search
spaces are randomly sampled between 2 and 5. As all the datasets are i.i.d.
samples, we use the first $16$ datasets to be training datasets and the
remaining $4$ datasets to be testing datasets in Setup A. Similarly, we use
the first $8$ sub-datasets of each dataset to form the training data and use
the remaining $2$ sub-datasets of each dataset to form the testing data in
Setup B.

HPO-B Super-dataset includes 16 datasets (search spaces), while each dataset
contains different numbers of sub-datasets and has its pre-specified
train/test split of sub-datasets. We use this pre-specified train/test split
for the experiments on HPO-B.

In the experiments, all GP-based methods use a constant mean function and a
Matérn kernel with smoothness parameter $\nu=3/2$. In our setting, each
dataset $D_{i}$ ($i=1,\dots,N$) corresponds to $\mathcal{GP}_{i}$
parameterized by $\theta_{i}$, which includes the following parameters:
constant mean (the value of the constant mean function), length-scale (which
has the same number of dimensions as the search space), signal variance. In
addition, $D_{i}$ is also parameterized by the noise variance $\sigma_{i}$.

For the prior distribution types of each of these parameters, we use a normal
distribution for constant mean and use a Gamma distribution for each of the
remaining parameters. A normal distribution is parameterized by $\mu$ (mean)
and $\sigma$ (standard deviation), while a Gamma distribution is parameterized
by $\alpha$ (shape) and $\beta$ (rate).

The specific distribution parameters used to generate the Synthetic Super-
dataset are as follows: constant mean is sampled from
$\mathrm{Normal}(\mu=1,\sigma=1)$, each dimension of length-scale is sampled
from $\mathrm{Gamma}(\alpha=10,\beta=30)$, signal variance is sampled from
$\mathrm{Gamma}(\alpha=1,\beta=1)$, and noise variance is sampled from
$\mathrm{Gamma}(\alpha=10,\beta=100000)$.

For the non-informative hierarchical GP, the specific parameters are as
follows: constant mean is sampled from $\mathrm{Uniform}(low=-100,high=100)$,
each dimension of length-scale is sampled from
$\mathrm{Uniform}(low=0.001,high=10)$, signal variance is sampled from
$\mathrm{Uniform}(low=0.000001,high=100)$, and noise variance is sampled from
$\mathrm{Uniform}(low=0.00000001,high=100)$.

For the hand-specified hierarchical GP, the specific parameters are as
follows: constant mean is sampled from $\mathrm{Normal}(\mu=0,\sigma=1)$, each
dimension of length-scale is sampled from $\mathrm{Gamma}(\alpha=1,\beta=10)$,
signal variance is sampled from $\mathrm{Gamma}(\alpha=1,\beta=5)$, and noise
variance is sampled from $\mathrm{Gamma}(\alpha=10,\beta=100)$.

We fit the GP parameters in each search space by minimizing Eq. 3 using L-BFGS
for the Synthetic Super-dataset, and using the Adam optimizer for HPO-B. For
L-BFGS on the Synthetic Super-dataset, the number of iterations is 500, and
all observations from a sub-dataset is used at each iteration of L-BFGS. For
Adam optimizer on HPO-B Super-dataset, the number of iterations is 10000, and
each sub-dataset is sub-sampled to 50 observations at each iteration of Adam
optimizer. The learning rate of Adam optimizer is 0.001. The reason for using
Adam in HPO-B experiments is that the sub-datasets in HPO-B can be too large
to be handled by L-BFGS. Furthermore, L-BFGS was recommended by Wang et al.,
(2022) as the standard objective optimization method for HyperBO, while Adam
was recommended by Wistuba and Grabocka, (2021) that applied the FSBO method
on HPO-B.

When testing the BO performances of compared methods, the BO budget is 50, and
5 random observations are made prior to BO for initialization. For
hierarchical GP-based methods, 100 samples of GP parameters are used to
compute the acquisition function value according to Eq. 5 at each BO
iteration. All methods are tested on 5 random seeds, while the set of initial
observations are the same for all methods given the same random seed. Each
sub-dataset of HPO-B are sub-sampled to 1000 observations for BO testing,
while no sub-sampling is needed for the Synthetic Super-dataset since each
sub-dataset contains 300 observations.

For evaluation of NLL values, each sub-dataset is sub-sampled to 100
observations due to the time complexity of computing NLL, and we report the
average NLL across 10 random seeds. The sub-sampling is the same for all
compared methods given the same random seed. For hierarchical GP-based
methods, 500 samples of GP parameters are used to compute the NLL value
according to Eq. LABEL:equ:estimate-nll.

##  Appendix D Additional experiment results

In this section, we present additional empirical evidences on the asymptotic
behaviors of our pre-training method, as well as analyses on the performance
on Bayesian optimization.

###  D.1 Asymptotic behavior of fitting a single GP

(a) Length-scale

(b) Constant Mean

(c) Signal Variance

(d) Noise Variance

Figure 11: Estimated length-scale, constant mean, signal variance and noise
variance parameters of a 1-dimensional GP as the number of sub-datasets
increases. Each sub-dataset contains 25 observations. The violin plot shows
the distribution of estimates for each number of sub-datasets over 50 random
seeds. The variance of estimated GP parameters decreases as the number of sub-
datasets increases.

Fig. 11 demonstrates the asymptotic behavior empirically for fitting the
parameters of a single GP using an increasing number of sub-datasets with
simulations on synthetic data generated with a 1-dimensional GP. The ground-
truth GP parameters are shown in the figure. The variance of estimated GP
parameters decreases as the number of sub-datasets increases.

Figure 12: NLL evaluations of training/testing data on HyperBO+ and HyperBO in
Setup B on the Synthetic Super-dataset.

###  D.2 NLL evaluations per search space in Setup B

Fig. 9 shows the NLL evaluations of training/testing data on HyperBO+ and
HyperBO for the HPO-B Super-dataset. HyperBO obtained high test NLLs despite
having much lower train NLLs in some search spaces (such as 5527, 5891, and
6794). This observation supports our conjecture in §6.3 that HyperBO overfits
in some search spaces of HPO-B. Fig. 12 shows the NLL evaluations of
training/testing data on HyperBO+ and HyperBO in Setup B on the Synthetic
Super-dataset. We can see that the overfitting does not occure on the
Synthetic Super-dataset as HyperBO+ has slightly higher NLL values than
HyperBO. The reason for this difference across the two super-datasets might be
that the Synthetic Super-datasets are generated so that all sub-datasets in a
search space are sampled from the same GP, while the HPO-B Super-dataset does
not guarantee this property as its data are collected from real-world
hyperparameter tuning tasks.

(a) Length-scale - $\alpha$ of Gamma Distribution

(b) Length-scale - Gamma Distribution

(c) Constant Mean - $\mu$ of Normal Distribution

(d) Constant Mean - Normal Distribution

(e) Signal Variance - $\alpha$ of Gamma Distribution

(f) Signal Variance - Gamma Distribution

(g) Noise Variance - $\alpha$ of Gamma Distribution

(h) Noise Variance - Gamma Distribution

Figure 13: Left: Estimated shape parameter $\hat{\alpha}$ (y-axis) of Gamma
distribution prior for the length-scale, signal variance and noise variance
parameters of a kernel function and estimated mean parameter $\hat{\mu}$
(y-axis) of the normal distribution prior that models the constant mean GP
parameter w.r.t. the number of training datasets (x-axis). We show the mean
and violin plots for $\hat{\alpha}$ over 5 random seeds. Right: The pre-
trained prior distributions with 2 training datasets (green) and 16 training
datasets (blue) are shown along with the ground-truth prior distribution (red)
for each GP parameter.

###  D.3 Asymptotic behavior of the two-step pre-training method

Fig. 13 shows on the left the estimated prior distribution parameter (mean
$\mu$ for normal distribution and shape $\alpha$ for Gamma distribution) for
GP parameters including length-scale, constant mean, signal variance and noise
variance w.r.t. the number of training datasets. This is demonstrated only on
the Synthetic Super-dataset since there is no known ground-truth prior
distributions for HPO-B Super-dataset.

In addition, the estimated prior distributions for each GP parameter with 2
training datasets and 16 training datasets along with the corresponding
ground-truth prior distribution are shown on the right of Fig. 13. For most of
the GP parameters (length-scale, constant mean and signal variance), the
estimated prior distribution parameter shows decreasing variance and gets
closer to the ground-truth value as the number of training datasets increases.

The estimated prior distribution also gets much closer to the ground-truth
prior distribution when 16 training datasets are used compared to only using 2
training datasets. For the noise variance parameter, the $\alpha$ parameter of
Gamma distribution is estimated to be a different value than the ground-truth
value, and the estimated prior distribution shows a different shape than the
ground-truth prior distribution even when using 16 training datasets. Note
that each dataset of the Synthetic Super-dataset contains only 10 sub-
datasets. This small number of sub-datasets in each search space is likely
part of the reason for the difference between the estimated noise variance
prior and the ground-truth prior. The variance of the estimated $\alpha$
parameter of Gamma distribution for noise variance still decreases as the
number of training datasets increases. The visualization shows that the
estimated prior distributions for noise variance also put most of the
probability density over values between [0, 0.0002] as the ground-truth prior
distribution does despite the difference in their shapes.

Figure 14: Average normalized simple regret of each method during BO in both
Setup A and Setup B for the Synthetic Super-dataset. The averages are taken
over 5 random seeds, and the highlighted areas show
$\mathrm{mean}\pm\mathrm{std}$ for each method. Results for acquisition
functions PI, EI, and UCB are shown here. Figure 15: Average normalized simple
regret of each method during BO for the HPO-B Super-dataset. The averages are
taken over 5 random seeds, and the highlighted areas show
$\mathrm{mean}\pm\mathrm{std}$ for each method. Results for acquisition
functions PI, EI, and UCB are shown here.

###  D.4 BO performances with different acquisition functions

Fig. 14 shows the average normalized simple regret of compared methods duing
BO iterations with different acquisition functions in both Setup A and Setup B
for the Synthetic Super-dataset and Fig. 15 shows the results with different
acquisition functions in HPO-B Super-dataset. The acquisition functions
included are Probability of Improvement (PI), Expected Improvement (EI), and
GP-UCB. The PI acquisition function has parameter $\zeta=0.1$, which means its
target value is the maximum observation plus $\zeta$. The EI acquisition
function uses the maximum observation as the target value. The GP-UCB
acquisition function has parameter $\beta=3$. It can be observed that PI
generally provides better BO performance for GP-based methods across the two
super-datasets than EI and GP-UCB. This is in line with the finding by Wang
and Jegelka, (2017) that PI can obtain high BO performance by setting good
target values. Therefore we use PI as the acquisition function for BO results
in the main paper.

We can see that given the same acquisition function, HyperBO+ outperforms the
non-informative hierarchical GP and the hand-specified hierarchical GP
baselines and achieves a performance that is either similar to or better than
that of HyperBO in most of the three cases (HPO-B Super-dataset and two setups
for the Synthetic Super-dataset).

Notice that HyperBO+ aims to recover the ground-truth prior distribution, and
as what we have shown here, we need a good acquisition function to unleash the
best performance of BO given our pre-trained prior. Using misspecified priors,
some sub-optimal acquisition functions may actually lead to better performance
than an acquisition function that achieves better results with the ground-
truth prior.