OPTIMIZATION WITH A CASE STUDY ON TUNING DEEP NEURAL NETWORK OPTIMIZERS
Year: 2021

Abstract:
The performance of deep neural networks can be highly sensitive to the choice of a variety of meta-parameters, such as optimizer parameters and model hyperparameters. Tuning these well, however, often requires extensive and costly experimentation. Bayesian optimization (BO) is a principled approach to solve such expensive hyperparameter tuning problems efficiently. Key to the performance of BO is specifying and refining a distribution over functions, which is used to reason about the optima of the underlying function being optimized. In this work, we consider the scenario where we have data from similar functions that allows us to specify a tighter distribution a priori. Specifically, we focus on the common but potentially costly task of tuning optimizer parameters for training neural networks. Building on the meta BO method from Wang et al. (2018b), we develop practical improvements that (a) boost its performance by leveraging tuning results on multiple tasks without requiring observations for the same meta-parameter points across all tasks, and (b) retain its regret bound for a special case of our method. As a result, we provide a coherent BO solution for iterative optimization of continuous optimizer parameters. To verify our approach in realistic model training setups, we collected a large multi-task hyperparameter tuning dataset by training tens of thousands of configurations of near-state-of-the-art models on popular image and text datasets, as well as a protein sequence dataset. Our results show that on average, our method is able to locate good hyperparameters at least 3 times more efficiently than the best competing methods.

text: The performance of deep neural networks can be highly sensitive to the choice of a variety of meta-parameters, such as optimizer parameters and model hyperparameters. Tuning these well, however, often requires extensive and costly experimentation. Bayesian optimization (BO) is a principled approach to solve such expensive hyperparameter tuning problems efficiently. Key to the performance of BO is specifying and refining a distribution over functions, which is used to reason about the optima of the underlying function being optimized. In this work, we consider the scenario where we have data from similar functions that allows us to specify a tighter distribution a priori. Specifically, we focus on the common but potentially costly task of tuning optimizer parameters for training neural networks. Building on the meta BO method from Wang et al. (2018b), we develop practical improvements that (a) boost its performance by leveraging tuning results on multiple tasks without requiring observations for the same meta-parameter points across all tasks, and (b) retain its regret bound for a special case of our method. As a result, we provide a coherent BO solution for iterative optimization of continuous optimizer parameters. To verify our approach in realistic model training setups, we collected a large multi-task hyperparameter tuning dataset by training tens of thousands of configurations of near-state-of-the-art models on popular image and text datasets, as well as a protein sequence dataset. Our results show that on average, our method is able to locate good hyperparameters at least 3 times more efficiently than the best competing methods.

heading: 1 INTRODUCTION
text: The careful tuning of a variety of meta-parameters, such as optimizer parameters and model hyperparameters has become a basic necessity for deep learning (Bergstra et al., 2011; Feurer et al., 2015). Such tuning requires extensive experimentation, retraining models repeatedly with different configurations, and can be challenging at realistic budgets because the tuning landscape is typically non-stationary, noisy and ill-behaved. Tuning has become sufficiently costly that finding more efficient and effective tuning procedures has the potential to save a substantial amount of resources or, alternatively, improve the accuracy of the final models at a given budget.
Some hyperparameters might be common across a large number of tuning problems, such as those pertaining to an optimization algorithm. For example, Adam is used across many deep learning applications and has four parameters that require careful tuning (Nado et al., 2021). Thus, if we have access to the performance of different optimizer-specific hyperparameters on different model training tasks, we may be able to transfer the knowledge among those tasks. This kind of meta-level learning is common among practitioners themselves: when faced with a new tuning problem, one might first try reusing hyperparameter settings that worked well on another problem. The underlying assumption is that hyperparameters should perform similarly across tasks. In this work, we aim to formalize this assumption and automate optimizer hyperparameter tuning by leveraging knowledge from previous experiments. Although our experiments consider optimizer parameter tuning as a practically important sub-problem of hyperparameter tuning for deep neural networks, our method applies to any hyperparameters that are common across multiple tasks.
Bayesian optimization (BO) has become a popular methodology for optimizing the hyperparameters of machine learning models (Snoek et al., 2012; Bergstra et al., 2011) and represents the state-ofthe-art (Turner et al., 2021). BO involves specifying a probabilistic model over the function to be optimized and using this to reason about the location of the optimum. The optimization proceeds by iteratively updating the model with new data and using the posterior distribution to reason about where to next evaluate, trading off exploration and exploitation. The model is typically specified with only a-priori assumptions of smoothness, for example using a Gaussian process (GP) with a smooth covariance function. Even if the model is well-specified, BO can be slow to converge due to the generality of the prior assumptions. This would seem wasteful for problems that are repeated often or share considerable structure with previous experiments.
One natural option is to cast our problem as meta Bayesian optimization, where the goal is to learn to optimize a function by generalizing from past experience with other similar functions. Indeed, several BO methods exist in the literature, but they are unsuitable for our scenario where we envision potentially thousands of related tasks within e.g. the context of a hyperparameter tuning service. Existing meta BO methods either scale cubically in the number of evaluations and tasks (Swersky et al., 2013; Bardenet et al., 2013) (See §4.3 for more details), impose a restrictive set of assumptions on the available data (Wang et al., 2018b; Swersky et al., 2013) to obtain efficient solutions, or make assumptions on the availability of GP parameters (Volpp et al., 2020) or descriptive task-level features (Brazdil et al., 1994; Bardenet et al., 2013; Yogatama & Mann, 2014).
To address these issues, we introduce HyperBO: a meta BO method that builds upon Wang et al. (2018b) with a relatively simple assumption: all the related functions being optimized are samples from the same GP prior distribution over functions. Concretely, HyperBO assumes the functions are conditionally independent given the hyperparameters, mean and covariance function of the GP. Compared to Wang et al. (2018b), HyperBO does not impose any strict conditions on data or model structures, and a special case of HyperBO retains strong regret bounds. From a computational perspective, HyperBO scales linearly in the number of tasks during training, and does not depend on the number of tasks when deployed. By not imposing assumptions about the data collection conditions, it can be used with large offline datasets or a few related optimization trajectories.
To evaluate HyperBO, we collected a large multi-task hyperparameter tuning dataset by training tens of thousands of configurations of near-state-of-the-art models on popular image and text datasets, as well as on a protein sequence dataset. We compare HyperBO to several hyperparameter tuning baselines in the sequential BO setting. Our results show that optimizers that use hyperparameters suggested by our method are able to obtain better performing models requiring at least 3 times fewer function evaluations than other baselines.
Our main contributions are two-fold: (1) a practical meta BO approach that makes minimal assumptions; and (2) a large multi-task hyperparameter tuning dataset that not only benefits our method but also serves as an ideal benchmark to test future multi-task or meta-learning BO methods.1
Related work There is a rich literature of innovative methodologies to improve the efficiency of BO given related tasks or additional context. Here we discuss the most closely related work and explain why these don’t solve the specific scenario which we envision. Specifically, our goal is a methodology that is scalable enough to share information across thousands of tasks, each with potentially hundreds of observations, such as in the context of a large BO service or library.
In this work we use the term meta-BO to refer to the class of BO methods that use data from existing tasks to optimize a new task. Multi-task BO (Swersky et al., 2013; Poloczek et al., 2017; Yogatama & Mann, 2014) and transfer learning BO using contextual GPs (Krause & Ong, 2011; Bardenet et al., 2013; Poloczek et al., 2016) are both meta BO approaches. Some meta BO methods have also been studied for hyperparameter tuning tasks in machine learning (Feurer et al., 2015).
HyperBO assumes all tasks are independent (after conditioning on the GP), whereas both multi-task and contextual BO rely heavily on the assumption that tasks are related. Thus the latter approaches typically scale cubically in both the number of tasks and observations in each task, meaning that they cannot gracefully scale across both without heavy approximations. When assuming that all inputs are equal across tasks, multi-task BO can be sped up using a Kronecker decomposition of the kernel
1We are working on open-sourcing the code base and dataset. The dataset is collected based on an opensourced code base (Gilmer et al., 2021).
to a task kernel and an input kernel which can be inverted separately; a similar assumption is made by Wang et al. (2018b). In comparison, HyperBO scales linearly in the number of tasks (see §B).
End-to-end learning (Chen et al., 2017; Volpp et al., 2020) is another popular meta BO approach for hyperparameter tuning that learns a strategy to suggest new query points based on past history of BO. One limitation of such approaches is that the total number of BO iterations must be determined a-priori. Furthermore, by nature of using a highly parameterized model to train the strategy, we lose the interpretability of intermediate steps that GPs and acquisition functions provide.
HyperBO builds upon Wang et al. (2018b) and Kim et al. (2017; 2019). We resolve their issues with optimizing over a continuous space rather than a discrete set and limitations on using the same set of inputs across tasks. Kim et al. (2017; 2019) estimated a multivariate Gaussian that models values for search strategies in robot manipulation tasks, and thus only considered discrete inputs. Wang et al. (2018b) provided regret bounds for Kim et al. (2017; 2019), which was identified as meta BO without the knowledge of the GP prior. For both finite discrete search spaces and continuous ones, Wang et al. (2018b) requires observations on the same set of inputs across tasks, which is an assumption that is not required for HyperBO; HyperBO still inherits the same regret bound as Wang et al. (2018b) for the special case where the same-inputs assumption is satisfied. Similar ideas also appeared in Perrone et al. (2018), which can be viewed as a special case of HyperBO or Wang et al. (2018b).

heading: 2 PROBLEM FORMULATION
text: We consider the standard black-box function optimization scenario: given a real-valued function f defined over a compact, hyper-rectangular space X ⊂ Rd and given observations of similar functions f1, · · · , fN , we seek an x ∈ X optimizing f . We inherit our problem formulation from Wang et al. (2018b), but we relax impractical assumptions on data availability (we do not require all observations to be made on the same inputs across tasks) and model restrictions. Assumptions and the goal. Concretely, we assume that there exists a Gaussian process GP(µ, k) with unknown mean function µ : X→ R and kernel k : X× X→ R. Let N be the number of tasks and let Mi be the number of observations we have for the ith task. Conditioned on independent function samples fi ∼ GP(µ, k) and inputs x(i)j ∈ X, i ∈ [N ], j ∈ [Mi], we observe evaluations y (i) j ∼ N (fi(x (i) j ), σ
2) perturbed by i.i.d. additive Gaussian noise N (0, σ2). Taken together, the collection of sub-datasets Dfi = {(x (i) j , y (i) j )} Mi j=1 define a dataset DN = {Dfi}Ni=1. Finally, our goal is to maximize a new function independently sampled from the same GP, f ∼ GP(µ, k); that is, solve arg maxx∈X f(x) given dataset DN but unknown functions µ, k and unknown parameter σ
2. An example. In our optimizer hyperparameter tuning application, a task corresponds to finding the best optimizer hyperparameters to train a given model on a particular dataset,2 e.g. training a ResNet (He et al., 2016) on ImageNet (Russakovsky et al., 2015). Notice that we do not assume that the mean function µ, kernel k and noise variance σ2 are given. This is consistent with the reality of solving real-world black-box optimization problems including hyperparameter tuning. We must learn those unknown functions and parameters from data. However, in practice, searching in functional spaces to find the right mean µ or kernel k is a daunting task. Hence for practical concerns, a well defined search space for functions is required. More details on this can be found at §3.1. Metrics. For simplicity, throughout this paper, we focus on the setting where the target function f can only be optimized by iteratively choosing where to evaluate, and defer batch evaluation setups to Sec. F. As we run BO on the target function f for T iterations, we accumulate a set of observations Df = {(xt, yt)}Tt=1, yt ∼ N (f(xt), σ2). We evaluate the quality of the optimization using the simple regret metric: RT = maxx∈X f(x)− f(x̂), where x̂ is the final recommendation at the end of the optimization process. There are various ways of setting x̂ based on the observations Df ; we use the input that achieved the best evaluation: x̂ = xτ ; τ = arg maxt∈[T ] yt. Bayesian viewpoint. As mentioned above, the observed functions f1, · · · , fN and the evaluation target f are assumed to be independent draws from the same GP. This assumption is consistent with a hierarchical Bayes interpretation (Fig. 1), where all observed functions are independent conditioned on the GP. Notice that for BO, each selected input x(i)j depends on all previous observations. But we only describe the generative model of a hierarchical GP for simplicity.
2Technically, we also consider different batch sizes to be different tasks.
More specifically, we assume that underlying functions of hyperparameter optimization tasks are defined by a parameter θ ∼ p(θ;α); mean and kernel functions µ and k are defined by deterministic functions parameterized by θ. The independent function samples {fi}i∈[N ] are draws from GP(µ, k). We then learn function f from observations DN on all other conditionally i.i.d. function samples f1, · · · , fN . We forgo a fully Bayesian approach that samples from the posterior over θ at every BO iteration. Our method can be viewed as a type-II maximum likelihood approximation of such a Bayesian solution. Notations. Let [n] denote {1, · · · , n},∀n ∈ Z+. For conciseness, we write the evaluation of a function f : X→ R on matrix x = [xi] n
i=1 ∈ Rn×d as µ(x) := [µ(xi)]ni=1. Similarly, for two matrices (i.e., the concatenation of input vectors) x,x′, we write the corresponding kernel matrix as k(x,x′) := [k(xi, x′j)]i∈[n],j∈[n′], and shorten k(x) := k(x,x). We denote a (multivariate) Gaussian distribution with mean u and variance Σ by N (u,Σ), and a Gaussian process (GP) with mean function µ and covariance function k by GP(µ, k). Let σ2 be the noise variance in observations. Given a set of observations D = {(xt, yt)}Tt=1,yT = [yt]Tt=1 ∼ N (f(xT ), σ2I),xT = [xt]Tt=1 and f ∼ GP(µ, k), we denote the corresponding conditional GP distribution as GP(µ, k | D). Let ψ(x) = k(x,xT )(k(xT ) + σ2I)−1. Recall that the conditional distribution, GP(µ, k | D) = GP(µD, kD), is given for any x, x′ ∈ X as
µD(x) = µ(x) + ψ(x)(yT − µ(xT )), kD(x, x′) = k(x, x′)− ψ(x)k(xT , x′). (1)
3 OUR METHOD
Algorithm 1 HyperBO with acquisition function α(·). 1: function HYPERBO (f,DN ) 2: GP(µ̂, k̂)← TRAINGP(DN ) 3: Df ← ∅ 4: for t = 1, · · · , T do 5: xt ← arg max x∈X α ( x;GP(µ̂, k̂ | Df )
) 6: yt ← OBSERVE(f(xt)) 7: Df ← Df ∪ {(xt, yt)} 8: end for 9: return Df
10: end function As shown in Alg. 1, our approach trains the GP hyperparameters on a representative set of datasets and fixes them for the duration of the optimization procedure; we refer to this approach as HyperBO. HyperBO runs in two steps. First, we learn a GP model GP(µ̂, k̂) to approximate the ground-truth (unknown) GP that generated the dataset DN . Then, we do standard BO to optimize a new function f with the learned GP GP(µ̂, k̂). The initial learning process (line 2) is the critical difference between HyperBO and standard BO algorithms, as well as the key contribution of this paper.
Based on the Bayesian graphical model interpretation (Fig. 1), our goal is to obtain a point estimate θ̂ for the parameter θ. Given this estimate, we can then estimate the mean function µ̂ and the kernel k̂, which defines our learned model GP(µ̂, k̂). During the BO iterations (Alg. 1, lines 4-8), we update the conditional GP, but do not re-estimate the GP mean and kernel. By separating the data for the conditional GP update and GP parameter training, we minimize the computational cost while still maintaining good performance both theoretically and empirically. Moreover, we avoid the BO chicken-and-egg dilemma (Wang et al., 2018b) where the search strategy is trained on data collected in the BO process and the data points are selected by the search strategy.
Next, we introduce our GP training strategy based on two types of objectives: marginal data likelihood (§ 3.1) and distance between estimates and model predictions (§ 3.2). In the appendix, §B reveals the complexity of HyperBO that is linear in the number of tasks and § 3.3 shows that a special case of HyperBO retains strong regret bounds (Wang et al., 2018b).

heading: 3.1 MARGINAL LIKELIHOOD
text: A straightforward way to train a GP is by optimizing the log marginal likelihood over the GP’s hyperparameters. This is also known as type II maximum likelihood approximation (Rasmussen
& Williams, 2006). In our case, we derive the data likelihood for the observations from multiple functions that are assumed to be given, which is a key difference to regular GP or BO setups. The log marginal likelihood for our method is
log p(DN | µ, k, σ2) = ∑N
i=1
( −1
2 ȳ>(i)K
−1ȳ(i) − 1 2 log |K| − Mi 2 log 2π
) , (2)
where ȳ(i) = y(i) − µ ( x(i) ) , K = k ( x(i) ) + σ2I , x(i) = [x(i)j ] Mi j=1 and y (i) = [y (i) j ] Mi j=1.
Our solution to the choice of mean function, kernel function and noise variance then becomes
µ̂, k̂, σ̂2 = arg maxµ,k,σ2 log p(DN | µ, k, σ2). (3)
For mean function µ and kernel k, this optimization is done in functional space. While methods exist to search for functional structures (Kemp & Tenenbaum, 2008; Malkomes & Garnett, 2018), one may opt for a simple search strategy within a group of functional structures (e.g. mean µ ∈ {linear, constant} and kernel k ∈ {exponentiated quadratic,Matérn}). For all combinations of mean/kernel structures or functional classes, we then optimize the parameterization of them and noise variance σ2 to eventually solve Eq. 3. Details of how we defined the search space can be found in §4.

heading: 3.2 DISTANCE BETWEEN ESTIMATES AND MODEL PREDICTIONS
text: Although the marginal likelihood is a straightforward objective to optimize, it may not be straightforward to interpret how high of a likelihood is high enough for us to stop our search for a decent model. Nevertheless, we may be able to directly estimate the sample mean and covariance, and the distance between those estimates and model predictions could be a good indicator of how good the model is. We will show in §3.3 that a distance objective may lead us to better theoretical properties.
Here we consider a special case of dataset DN where part of it has matching inputs across some sampled functions. More formally, suppose we have a matching dataset D′N = {(xj ,yj)}Mj=1 where M is a positive integer, xj ∈ X,yj = [y(i)j ]Ni=1 ∈ RN and y (i) j ∼ N (f(xj), σ2). Empirically, dataset D′N can be constructed by querying a set of functions f1, · · · , fN at the same set of input locations x = [xj ] M j=1 ∈ RM×d to obtain an observation matrix y = [yj ]Mj=1 ∈ RM×N .
By definition of a GP, the vector of all function queries f(x) is distributed according to a multivariate Gaussian distribution N (µ(x), k(x)). With our observation model, we get the distribution for observations y ∼ N (µ(x), k(x) + Iσ2) for some unknown mean function µ and kernel k. However, given that we have access to all observations y, we can estimate the mean on inputs x as µ̃ = 1N y1N ∈ R M and estimated covariance as K̃ = 1N (y − µ̃1 > N )(y − µ̃1>N )> ∈ RM×M ; here 1N is a column vector of size N filled with 1s. We use a biased estimate of covariance to be consistent with the corresponding maximum likelihood estimator in Eq. 3. 3 Notice that the estimated covariance includes in diagonal terms the variance of the observation noise.
For any distance function between the estimateN (µ̃, K̃) and model predictionN (µ(x), k(x)+Iσ2), we obtain an objective to minimize, D ( N (µ̃, K̃),N (µ(x), k(x) + Iσ2) ) . While there are different
measures of distributional discrepancy, we adopt the KL divergence. Let µ = µ(x) and K = k(x) + Iσ2. The KL divergence is defined as
DKL ( N (µ̃, K̃),N (µ,K) ) = 1
2
( tr(K−1K̃) + (µ− µ̃)>K−1(µ− µ̃) + ln |K|
|K̃| −M
) , (4)
and we can estimate the mean, kernel and noise variance by minimizing DKL. While it is difficult to gauge how much a probability density is enough to obtain a good model, Eq. 4 is a “distance” that goes to 0 as the difference between two distributions reduces. One may choose to do early stopping or model selection based on how close Eq. 4 is to 0. Through information theory, we also know that the KL divergence in Eq. 4 describes the number of extra bits (or nats) to encode the multivariate normalN (µ̃, K̃). Overall we found the KL divergence in Eq. 4 relatively more interpretable than the marginal likelihood in Eq. 3.
3One may choose to re-scale learned kernel by N N−1 to be unbiased.
The KL divergence in Eq. 4 introduces a different optimization landscape than the marginal likelihood in Eq. 3. The KL divergence also makes use of the matching dataset D′N in a way that the marginal likelihood cannot. In fact, all matching inputs in the marginal likelihood in Eq. 3 are implicit: all inputs are passed in to mean/kernel functions, and so there is no way that Eq. 3 can be informed that some inputs are the same across tasks. As shown in §4, the KL divergence in Eq. 4 interestingly led to better results in our experiments.

heading: 3.3 THEORETICAL ANALYSES
text: While it is nontrivial to prove regret bounds for general scenarios without strict assumptions, it is straightforward to show a regret bound for our method with objective DKL of Eq. 4 in the matching dataset case where BO is running on a finite set of inputs (Wang et al., 2018b). Theorem 1. Let N ≥ 4 log 6δ + T + 2. With probability at least 1− δ, simple regret in T iterations of Alg. 1 with special cases of either GP-UCB or PI satisfies
RT < O
(√ 1
N − T +
√ log 1
δ
) O ( 1
2T max A⊂X,|A|=T log |I + σ−2k(A)|+ σ
) . (5)
More details can be found at §D. Theorem 1 shows that the regret bound has a linear dependency on the observation noise σ. This is expected because in practice, we select the best observation rather than best function value (before observing a noisy version of it) to compute the simple regret. Another reason is that we learn the noise parameter σ jointly with the kernel, as shown by Eq. 4. Hence when computing acquisition functions, the noise σ is always included in the predicted variance.
Intuitively, the more sub-datasets we have in the dataset, the larger N is, the better we are able to estimate the GP model, and the closer the regret bound is to the case where the GP model is assumed known. Interestingly, the number of BO iterations T makes the regret smaller in the second term but larger in the first term in Eq. 5. Usually as we get more observations, we get more information about the maximizer, and we are able to optimize the function better. However, as we get more observations on the new function, GP conditional predictions have more freedom to deviate from the ground truth (see Lemma 1 of Wang et al. (2018b)). As a result, we get less and less confident about our predictions, which is eventually reflected in a looser regret upper bound.
It is tempting to prove similar bounds for more general settings where inputs are not the same across all sub-datasets and BO happens in continuous space. Though the only prerequisite is to show that the difference between the learned mean/kernel and the ground truth mean/kernel is small, this prerequisite is as difficult as showing we can find a model that has bounded generalization error across the entire continuous input space of an arbitrary function. Instead of making unrealistic assumptions just to satisfy such prerequisite, we leave the regret bound for general settings as an open question.

heading: 4 EXPERIMENTS
text: Our goal in this paper is to provide a practical approach for hyperparameter optimization when we are given data on a range of tasks over the same search space. To analyze the effectiveness of our proposal, we take the optimizer hyperparameter tuning problem in deep learning as a case study. Our implementation of HyperBO is based on JAX (Bradbury et al., 2018).4
To reduce ambiguity, we distinguish between datasets that individual neural networks are trained on and the dataset we collected that includes optimizer hyperparameter points with their validation errors (and other metrics). We will call the former (e.g. MNIST, CIFAR10) task datasets and call the latter the tuning dataset. The tuning dataset is what we described as dataset DN in §2.

heading: 4.1 HYPERPARAMETER TUNING DATASET
text: In order to collect our hyperparameter tuning dataset, the PD1 Neural Net Tuning Dataset, we defined a set of 24 neural network tuning tasks5 and a single, broad search space for Nesterov
4We are working on open-sourcing our code as well as trained GP models. 5The batch size 1024 ResNet50 ImageNet task only has 100 hyperparameter points because we abandoned it
when scaling up data collection in order to save compute resources. It is used in training, but not evaluation.
momentum. Each task is defined by a task dataset (e.g. ImageNet), a specific neural network model (e.g. ResNet50), and a batch size. Tab. 1 shows all the tasks that we consider in the tuning dataset. We used an existing code base (Gilmer et al., 2021) for neural network model training. The dataset used roughly 12,000 machine-days of computation for approximately 50,000 hyperparameter evaluations.
For each task, we trained the model on the task dataset repeatedly using Nesterov momentum (Nesterov, 1983; Sutskever et al., 2013), with the task’s minibatch size, with different hyperparameter settings drawn from the 4-dimensional search space detailed in Tab. 2. We tuned the base learning rate, η, on a log scale, the momentum, β, with 1− β on a log scale, and the polynomial learning rate decay schedule power p and decay steps fraction λ. We used a polynomial decay schedule with the
following form: ητ = η1000 + ( η − η1000 ) ( 1− min(τ,λT )λT )p , where τ is the training step and T is the total number of training steps for the task.
Table 1: Tasks
Task Dataset Model Batch Sizes CIFAR10 Wide ResNet {256, 2048} CIFAR100 Wide ResNet {256, 2048} Fashion MNIST Max pool CNN ReLU {256, 2048} Fashion MNIST Max pool CNN tanh {256, 2048} Fashion MNIST Simple CNN {256, 2048} ImageNet ResNet50 {512, 1024, 2048} LM1B Transformer {2048} MNIST Max pool CNN relu {256, 2048} MNIST Max pool CNN tanh {256, 2048} MNIST Simple CNN {256, 2048} SVHN (no extra) Wide ResNet {256, 1024} WMT15 German-English xformer {64} uniref50 Transformer {128}
Table 2: 4-dimensional input search space (see text for more details)
Hyperparameter Range Scaling η [10−5, 10] Log p [0.1, 2.0] Linear
1− β [10−3, 1.0] Log λ [0.01, 0.99] Linear
We collected two types of data: matched and unmatched data. Matched data used the same set of uniformly-sampled hyperparameter points across all tasks and unmatched data sampled new points for each task. All other training pipeline hyperparameters were fixed to hand-selected, task-specific default values. All of our tasks are classification problems, so they all used the same training loss, although occasionally task-specific regularization terms were added. For each trial (training run for a single hyperparameter point), we recorded validation error (both cross entropy error and misclassification rate). In many cases, poor optimizer hyperparameter choices can cause training to diverge. We detected divergent training when the training cost became NaN and then marked the trial but did not discard it. Please see the Appendix, supplementary material, and code (Onomous, 2021) for additional details about the tasks and training procedure. The different tuning tasks vary in difficulty and numbers of data points, but generally there are roughly 500 matched datapoints and 1500 unmatched datapoints per tuning task. For unmatched data only, we attempted to generate roughly similar numbers of non-divergent points across tasks, so tasks with a higher probability of sampling a hyperparameter point that causes training to diverge will tend to have more trials.

heading: 4.2 DESCRIPTION OF ALL COMPARED METHODS
text: Our method HyperBO has several variants including using different acquisition functions and different objectives. In §4, unless otherwise mentioned, we used a thresholded probability of improvement (PI) as the acquisition function. We set PI in line 5 of Alg. 1 as µ̂Df (x)−maxt(yt+0.1)
σ̂Df (x) . We empirically
evaluated a variety of acquisition functions, but found PI thresholded at 0.1 to be surprisingly effective. Because we model the observations as log error rate, this actually trades off exploration and exploitation - i.e. with larger error rates this seeks relatively more substantial improvements than with small error rates. The list of 5 different acquisition functions we tested is as follows: PI with 0.1 threshold, expected improvement and UCB with 2, 3, 4 coefficients. See their comparisons at §E.5.
We use H* NLL to denote HyperBO with negative log marginal likelihood as the objective and H* KL to denote HyperBO with KL divergence on matching datapoints as objective. Both objectives are optimized via L-BFGS (Nocedal, 1980) with µ(x) = θ>0 tanh(θ > 1 x), k(x, x
′) = Matérn32(tanh(θ>1 x), tanh(θ > 1 x ′)), θ1 ∈ R4×8. These two settings of HyperBO are relatively representative of the performance of variants of HyperBO. See comparisons over objectives at §E.4.
Our baselines include (a) Rand: Random search in the corresponding scaled space in Tab. 2. (b) STBO: Single-task BO where in every BO iteration, STBO optimizes the GP hyperparameters via marginal likelihood on data of the test task. This implementation corresponds to the basic off-the-shelf BO setups. (c) STBOH: Single-task GP-UCB with a hand-tuned prior on hyper-parameters including
UCB coefficient. (d) MIMO: Multi-task BO with GP bases as an ensemble of feedforward neural networks with shared subnetworks (Havasi et al., 2020). (e) RFGP: Multi-task BO with GP bases as random features. More details can be found at §E.1.

heading: 4.3 RESULTS ON OFFLINE OPTIMIZER HYPERPARAMETER TUNING TASKS
text: Many tasks in §4.1 can use up a lot of compute resources and time, which makes it infeasible to perform a wide variety of experiments to analyze the characteristics of BO methods. Hence we adopt an offline approximation, which runs BO only on the finite set of points that each tuning sub-dataset contains. More details and analyses are available in Appendix E.2. Holding out relevant tasks. Fig. 2 (a) shows the performance profiles, the fraction of all test tasks that each method from §4.2 is able to outperform a baseline criterion at each BO iteration. We can see that MIMO is able to outperform other methods in the beginning 20 BO iterations, but its leading position soon gets surpassed by HyperBO (H* NLL and H* KL). Fig. 2 (b,c) illustrates the BO convergence curves of all competing methods, together with the vertical slice at the 100th iterations. RFGP and STBO are both falling much behind Rand. STBO trains the GP on the data that the GP suggests to query. Optimizing the marginal data likelihood on at most 100 datapoints in fact may not lead to a better model than random initialization (see Tab. 5 in §F). Surprisingly, the contextual information learned by RFGP did not generalize to a new task. On the other hand, MIMO is able to obtain a slightly better error rate than STBOH. Overall, learning the GP prior through data as with HyperBO outperforms other meta BO methods, and is a more principled and effective way to obtain the GP prior when compared with hand-tuning.
Effect of number of training tasks. We now investigate the impact of the number of training tasks on the performance of meta BO methods. In Fig 3 we show the BO simple regrets on tasks from Table 1 (except ImageNet ResNet50 2048) that use meta BO models trained on different number of training tasks. To analyze the performance of all methods on less-related tasks, we first remove training tasks that have the same task dataset as our current tuning task for testing, and then remove randomly selected training datasets from the rest.
HyperBO variants reduced the simple regret as more training tasks are given. Interestingly, H* NLL and H* KL are already slightly better than Rand and STBOH when they started off with only 3 training tasks. There are reasonable fluctuations in the results but overall the trend of regret is going down as the number of training tasks increases. MIMO also reduced regret when the number of tasks increased from 8 to 18. RFGP, however, fails to learn from training tasks possibly because it did not learn good task embeddings for GP regression models.
Figure 4: Medians and 20/80 percentiles of simple regrets for methods that uses models trained on 0.2% to 90% of data in each task.
Effect of number of data points in training tasks. One remaining question is, how does Mi in §2, the number of data points in each training tasks, affect the performance of meta BO methods. We analyze the impact of Mi by removing a portion of all data that we have access to for each task. In particular, we set the percentage of remaining data to be 0.2%, 0.5%, 1%, 3%, 5%, 10%, 30%, 50%, 70%, 90%. Remaining datapoints are selected uniformly randomly, which breaks the structure of matching data. Hence we do not include H* KL in this comparison, as H* KL only makes use of matching data.
Fig. 4 shows how the simple regret changes as the fraction of training data grows. Below 10% training data, we observe clear trend that more data lead to lower regret for both H* NLL and MIMO, and relatively no change for RFGP. We also found that the performance of HyperBO (H* NLL) does not change much as the fraction of training data increases from 5% to 90%. However, MIMO and RFGP suffers significantly from more data as the fraction of training data increases from 5% to 50%. It is not entirely clear why MIMO and RFGP have such behaviors. One conjecture is that neural network based Bayesian linear regression models may get too confident once the amount of data reaches a certain threshold. This means much less exploration if those models are used for BO.

heading: 4.4 RESULTS ON ONLINE OPTIMIZER HYPERPARAMETER TUNING TASKS
text: We now evaluate HyperBO methods, H* NLL and H* NLLKL (NLL plus KL divergence on matching datapoints as the objective), in the online setting, where we optimize over the full hypercube and some hyperparameters may be infeasible to evaluate. See full results at §E.3. Fig. 5 shows that it can be difficult for STBO or MIMO to recover from a “bad” datapoint, but HyperBO methods are robust and performed the best among all methods being compared.

heading: 5 DISCUSSION AND CONCLUSION
text: While we focused on obtaining a better prior in BO in this work, the following directions are orthogonal to what we studied: different search spaces across tasks, batch evaluation, high-dimensional or large scale data, etc. However, it should be straightforward to combine their solutions with HyperBO. Please find more discussions at §F together with implications of our assumptions.
HyperBO is a novel meta BO approach that supports practical applications that involve continuous inputs queried at possibly non-aligned locations across tasks. HyperBO uses a simple yet effective idea that is easy to implement and efficient to run. We evaluated HyperBO on real-world big model optimizer tuning tasks, and the results demonstrated its superior performance over state-of-the-art competing methods.

heading: A HIERARCHICAL GP, A BAYESIAN PERSPECTIVE
text: In Bayesian viewpoint of §2, we assume that the overall setting of the hyperparameter optimization task is defined by a parameter θ ∼ p(θ;α); mean and kernel functions µ and k are drawn from p(µ, k | θ). The independent function samples {fi}i∈[N ] are themselves draws from GP(µ, k). The generative story is as follows:
• Draw GP parameter θ from p(θ;α) and observation noise parameter σ from p(σ;α). • Draw mean function µ and kernel function k from p(µ, k | θ). • For each task i from 1 to N ,
– Draw a function fi from GP(µ, k). – For each data point j from 1 to Mi, ∗ Given input x(i)j , we draw the observation y (i) j ∼ N (fi(x (i) j ), σ 2).
We simplify this hierarchical setting by defining p(µ, k | θ) to be a sum of Dirac delta functions: both mean function µ and kernel k are deterministic functions parameterized by θ. Thus, we can infer GP parameter θ and noise σ from their posterior p(θ, σ | DN ∪Df ;α) and obtain an informed prediction for the target function
p(f | DN ∪Df ) = ∫ θ p(f | θ)p(θ | DN ∪Df ;α)
= ∫ θ p(f | θ) ∫ σ p(θ, σ | DN ∪Df ;α)
In other words, we learn function f from observations on all other conditionally i.i.d. function samples f1, · · · , fN . We forgo a fully Bayesian approach that samples from the posterior over θ at every BO iteration, although our method, HyperBO, can be viewed as a type-II maximum likelihood approximation of such a Bayesian solution.

heading: B COMPUTATIONAL COMPLEXITY
text: The marginal likelihood in Eq. 2 naturally decomposes into a sum of GP data likelihood terms on each sub-dataset Dfi . The time complexity to compute Eq. 2 is O(M3N), where N is the number of sub-datasets and M = maxNi=1Mi is the maximum number of data points for these sub-datasets. Notice that our method scales linearly in the number of tasks, N , in contrast to the cubic O(M3N3) scaling of multi-task or contextual BO methods (Swersky et al., 2013; Bardenet et al., 2013; Poloczek et al., 2016; Yogatama & Mann, 2014). The only cubic cost of HyperBO is on the number of data points in sub-datasets.
To train a GP with K optimization steps on Eq. 2, the time complexity is O(M3NK). The distance regularizers introduced in §3.2 requires estimating mean and covariance, which takes O(M2N) for matrix multiplication. The KL divergence in Eq. 4 has complexity O(M3) to compute and O(M3K) to optimize.
If there is any better probabilistic model than a GP to fit the data with less compute time, we can easily swap it in and reduce the O(M3) complexity that the GP contributed to the O(M3N) complexity of Eq. 2. For example, if we approximate a GP with a linear model on V random features (Rahimi et al., 2007), the complexity of Eq. 2 becomes O(V 3N). Another example is to train Eq. 2 with stochastic optimization methods, where the complexity of Eq. 2 on the full dataset can be reduced to O(B2MN), where B is the mini-batch size. Running stochastic optimization will then take O(B2MNK), where K is the number of optimization epochs.

heading: C OBJECTIVE FUNCTIONS
text: In §3, we presented NLL and KL divergence as objectives. Below we derive the KL divergence between a regular multivariate Gaussian and a degenerate multivariate Gaussian, which is the case for most of our matching data settings in §4.1: the number of matching data points is greater than the number of training tasks. In the end of this section, we introduce a new objective function, combining NLL and KL, interpreted as MAP with a data-dependent prior.
KL divergence for a degenerate multivariate Gaussian Eq. 4 of §3.2 gives the KL divergence between two Gaussians in the non-degenerate case. In practice, when we minimize Eq. 4, we can simply remove the constants and do the following
µ̂, k̂, σ̂2 = arg minµ,k,σ2DKL ( N (µ̃, K̃),N (µ,K) ) = arg minµ,k,σ2 1
2
( tr(K−1K̃) + (µ− µ̃)>K−1(µ− µ̃) + ln |K|
|K̃| −M ) = arg minµ,k,σ2 tr(K
−1K̃) + (µ− µ̃)>K−1(µ− µ̃) + ln |K|. (6) Here the variables we care about, µ, k, σ2, only appear in mean vector µ and covariance matrix K over the matching data. Even if the sample mean and covariance estimate N (µ̃, K̃) is degenerate, the optimization objective stays the same as reflected by the derivation below.
If N (µ̃, K̃) is degenerate, its base measure is at most N -dimensional rather than M -dimensional, given that there exists a full rank matrix A ∈ RM×R such that K̃ = AA> (R ≤ N ). Note that M is the number of matching data points, N the number of training tasks and R is the rank of matrix A and K̃. The KL divergence DKL ( N (µ̃, K̃),N (µ,K) ) is not well-defined because the base measure of
N (µ̃, K̃) is different from the base measure of N (µ,K), given K is full-rank. However, it is still possible to derive a pseudo KL divergence D∗KL ( N (µ̃, K̃),N (µ,K) ) as below.
Let the degenerate Gaussian be p(x) = N (µ̃, K̃) = |2πK̃|− 1 2 ∗ exp ( − 12 (x− µ̃)K̃ +(x− µ̃)> ) and
the non-degenerate one be q(x) = N (µ,K), where | · |∗ is the pseudo-determinant and K̃+ the peudo-inverse of K̃. We define the support of distribution p as S(p) = {µ̃+ K̃ 12 v | v ∈ RM}. The pseudo KL divergence between p(x) and q(x) now becomes
D∗KL ( N (µ̃, K̃),N (µ, K) ) =
∫ S(p) p(x) (ln p(x)− ln q(x))
= − 1
2 ∫ S(p) p(x) ( ln |2πK̃|∗ − ln |2πK|+ (x− µ̃)>K̃+(x− µ̃)− (x− µ)>K−1(x− µ) ) = 1
2
( (M − R) ln 2π + ln
|K| |A>A|
− Ep[tr(K̃+(x− µ̃)(x− µ̃)>) + tr(K−1(x− µ)(x− µ)>)] )
= 1
2
( (M − R) ln 2π + ln
|K| |A>A|
− tr(K̃+K̃) + Ep[tr(K−1(x− µ)(x− µ)>)] )
= 1
2
( (M − R) ln 2π + ln
|K| |A>A|
− tr(AA+) + Ep[tr(K−1(x− µ̃)(x− µ̃)> +K−1(2xµ̃> − 2xµ> − µ̃µ̃> + µµ>))] )
= 1
2
( (M − R) ln 2π + ln
|K| |A>A|
− R + tr(K−1K̃ +K−1(µ̃− µ)(µ̃− µ)>) )
= 1
2
( tr(K−1K̃) + (µ− µ̃)>K−1(µ− µ̃) + ln |K| − ln |A>A| − R + (M − R) ln 2π ) . (7)
If the covariance matrix K̃ is in fact full rank, i.e. R = M , pseudo KL in Eq. 7 then recovers the KL divergence DKL ( N (µ̃, K̃),N (µ,K) ) in Eq. 4 for non-degenerate Gaussians. If we were to
minimize this pseudo KL divergenceD∗KL ( N (µ̃, K̃),N (µ,K) ) , we still would get the minimization task in Eq. 6. However, pseudo KL divergence D∗KL does not satisfy properties including nonnegativity. In practice, we may also choose to add small epsilon values to the diagonal terms of both K and K̃ to make K̃ “less degenerate” and enable the existence of DKL ( N (µ̃, K̃),N (µ,K) ) .
Combining marginal likelihood and distance metrics It is also possible to additively combine the KL divergence with the negative log marginal likelihood objective, and treat this distance as a regularizer. In the case of KL divergence, it is equivalent to adding a data-dependent prior on the GP itself: µ,K ∼ exp(−λDKL(N (µ̃,K̃),N (µ,K)))/Z for some normalization constant Z, and the posterior is
p(µ, k, σ2 | DN ; µ̃, K̃) ∝ p(DN | µ, k, σ2) exp(−λDKL ( N (µ̃, K̃),N (µ,K) ) ). (8)
We can then obtain an MAP estimation for the unknown functions and variables µ, k, σ2.

heading: D DETAILS OF REGRET BOUNDS
text: Theorem 1 is a direct result of Theorem 16 in Wang et al. (2018b). The only subtle difference is that we adopted a biased estimate of the covariance matrix, a factor of NN−1 different from the unbiased estimate. But this can be corrected in the acquisition functions. We provide more details below. Proposition 2. For any M,d,N ∈ Z+,x ∈ RM×d,µ ∈ RM , V ∈ RN×M and K = V >V , there exists a Gaussian process GP(µ̂, k̂) such that DKL ( N (µ,K),N (µ̂(x), k̂(x)) ) ≡ 0.
Proposition 2 is easy to show. We can train a simple memory based model for mean function µ̂ and kernel k̂. The model stores each element of vector µ and matrix K at the corresponding locations in input x. When making a prediction on a new input x′ ∈ Rd, the model simply retrieves the values of the closest element in x. Given Proposition 2, a regret bound follows (Wang et al., 2018b).
By Proposition 2, we are able to obtain a GP(µ, k) where µ(x) and k(x) are equal to the sample mean and covariance on a matching dataset (x,y), following the notations in §3.2. Hence, our problem setup is consistent with the case with discrete input space in Wang et al. (2018b). The following theorem is a rewrite of Theorem 16 in Wang et al. (2018b), taking into account our biased estimators. Theorem 3. Assume there exist constant c ≥ maxx∈X k(x) and a training dataset is available whose size is N ≥ 4 log 6δ + T + 2. Define
ιt−1 = √√√√6(N − 3 + t+ 2√t log 6δ + 2 log 6δ) δN(N − t− 1) , bt−1 = 1 N − t log 6 δ , for any t ∈ [T ],
and ρT = max A∈X,|A|=T 1 2 log |I + σ −2k(A)|. Then, with probability at least 1 − δ, the best-sample simple regret in T iterations of meta BO with GP-UCB that uses
ζt =
( 6N(N − 3 + t+ 2 √ t log 6δ + 2 log 6 δ )/(δN(N − t− 1)) ) 1 2
+ (2N log( 3δ )) 1 2(
(N − 1)(1− 2( 1N−t log 6 δ )
1 2 ) ) 1 2
(9)
as its hyperparameter in αGP-UCBt−1 (x) = µ̂t−1(x) + ζtk̂t−1(x) 1 2 satisfies rGP-UCBT ≤ ηGP-UCB √
2cρT T log(1 + cσ−2)
+ σ2 − (2 log(3δ )) 1 2σ2
√ c+ σ2 ,
where ηGP-UCB = ( ιT−1+(2 log( 3 δ )) 1 2√
1−2 √ bT−1
√ 1 + 2 √ bT−1 + 2bT−1 + ιT−1 + (2 log( 3 δ )) 1 2 ).
With probability at least 1− δ, the best-sample simple regret in T iterations of meta BO with PI that uses f̂∗ ≥ maxx∈X f(x) as its target value satisfies
rPIT < η PI
√ 2cρT
T log(1 + cσ−2) + σ2 −
(2 log( 32δ )) 1 2σ2
2 √ c+ σ2 ,
where ηPI = ( f̂ ∗−µτ−1(x∗)√ kτ−1(x∗)+σ2 + ιτ−1)
√ 1+2b 1 2 τ−1+2bτ−1
1−2b 1 2 τ−1
+ ιτ−1 + (2 log( 3 2δ )) 1 2 , τ =
arg mint∈[T ] kt−1(xt).
The proof can be found in Wang et al. (2018b). Theorem 1 is a condensed version of Theorem 3.
While Theorem 3 provides us with some understanding of HyperBO in a specific setting, in practice, we need to query in a continuous input space that goes beyond the finite set of points present in the training dataset. It may or may not be possible to obtain data on a wide range of tasks to ensure N ≥ 4 log 6δ + T + 2. In fact, in all of our experiment, this criterion on number of tasks is not satisfied. However, we still obtained good performance.

heading: E EXPERIMENT DETAILS AND MORE RESULTS
text: In this section, we provide more detailed setups and empirical results on the impact of objective functions and acquisition functions in HyperBO. All experiment setups are the same as §E.2.1: offline and holding out related tasks.

heading: E.1 BASELINES
text: Our first set of baselines include those that do not use information from training tasks:
• Rand: Random search in the corresponding scaled space in Tab. 2. • STBO: Single task BO with a constant mean function, Matern32 kernel and PI acquisition
function (same as above). Every BO iteration, STBO optimizes the GP hyperparameters via marginal likelihood on data of the test task. This implementation corresponds to the basic off-the-shelf BO setups. • STBOH: Single task GP-UCB (coefficient=1.8) with constant mean, Matern52 kernel and
hand-tuned prior on hyper-parameters including UCB coefficient. Specifically, log amplitude follows Normal(-1, 1), log length scale (one per input parameter) follows Normal(0,1), and log observation noise variance follows Normal(-6, 3). The hyperparameters are postprocessed by tensorflow-probability’s SoftClip bijector to constrain the values between 1-st and 99-th quantiles. These prior distributions were manually tuned to obtain reasonable convergence rates on 24 analytic functions in COCO (Hansen et al., 2021). The GP parameters are then optimized via maximum marginal likelihood every BO iteration.
For multi-task BO baselines, we included scalable methods that replace the GP with a regression model that can be trained using SGD and thus scales linearly in the number of observations. Following the multi-task setup of Springenberg et al. (2016), we jointly trained a 5-dimensional embedding of each task, which was then added to the input of the following two models.
• MIMO: We trained an ensemble of feedforward neural networks with shared subnetworks (Havasi et al., 2020). We used 1 shared dense layer of size 10 and 2 unshared layers of size 10. We used tanh activation based on (Snoek et al., 2015, Figure 2). The network has one output unit with linear activation and another with softmax(10−4, 1) activation, corresponding respectively to the mean and standard deviation parameters of a normal distribution. We trained for 1000 epochs using the Adam optimizer with learning rate 10−4 and batch size 64. • RFGP: We used the open-source implementation of approximate GP by Liu et al. (2020).
We trained for 1000 epochs using the Adam optimizer with learning rate 10−3 and batch size 64.
All methods share the same input and output warping. The input warping is done according to the scaling function in Tab. 2: η ← log η, 1− β ← log(1− β). The output warping is done for the best validation error rate r ← − log(r + 10−10).

heading: E.2 OFFLINE BO EXPERIMENTS
text: In all offline experiments, we ran offline BO on the data from the test task starting from zero initial data from this task. Each method was repeated 5 times with different random seeds to initialize its model. We ran all methods without de-duplication to best simulate online BO. We evaluate on regret on error rate which denotes the simple regret on the finite set of data points in each tuning sub-dataset.

heading: E.2.1 HOLDING OUT RELEVANT TASKS
text: We first conducted experiments in a setting where a new task dataset is presented, and a BO method is trying to tune the optimizer hyperparameters for a selected model on that task dataset. A training dataset for meta BO is composed of at most 18 tuning sub-datasets on training tasks that do not involve the same task dataset as the test task. All methods then proceed to solve the test task on the new task dataset. Fig. 6 shows performance profiles of the BO methods described in §4.2. The performance profiles show the fraction of all test tasks that each method is able to outperform a baseline criterion at each BO iteration.6 We chose the criterion to be the median of best error rates achieved by all methods at 3 different BO iterations: 25th, 50th or 100th. The larger the fraction of tasks at each BO iteration, the better the method is. From all 3 criteria, we can see that MIMO is able to outperform other methods in the beginning 10 to 20 BO iterations, but its leading position soon gets surpassed by HyperBO (H* NLL and H* KL) and STBOH. HyperBO methods are gaining a similar if not larger fraction than the best alternative, STBOH, throughout BO iterations. Fig. 6 (c) has the most stringent performance criterion, and it shows that HyperBO with the KL objective outperforms HyperBO with the NLL objective in this set of experiments with a small margin. And both methods in HyperBO are doing considerably better than others.
Fig. 7 illustrates the BO convergence curves of all competing methods, together with the vertical slices at the 50th and 100th iterations. RFGP and STBO are both falling much behind Random search. STBO trains the GP on the data that the GP suggests to query, which creates a loop that could be harmful for data acquisition. Optimizing the marginal data likelihood on at most 100 datapoints in fact may not lead to a better model than random initialization (see Tab. 5 in §F). Surprisingly, RFGP, though equipped with the tuning dataset and initially reached some good values, performed similarly to STBO in the end. Clearly, the contextual information learned by RFGP did not generalize to a new task. On the other hand, MIMO is able to obtain a slightly better error rate than STBOH.
Fig. 6 and Fig. 7 both show that learning the GP prior through data like what HyperBO does is performing much better than other meta BO methods, and it is a more principled and effective
6We show performance relative to a baseline because of varying scales across tasks.
approach to obtain the GP prior than hand-tuning. As a reference, we include Tab. 3 which shows the task-wise best validation error rates obtained by the top 5 methods in 100 BO iterations.
To more precisely quantify HyperBO’s advantage, we also computed how much faster HyperBO can get a better error rate than best alternatives, which can be different from task to task. We found that on average, on over 50% tasks, H* NLL is at least 2.86 times faster than best non-HyperBO alternatives; while on over 57% tasks, H* KL is at least 3.26 times faster than best non-HyperBO alternatives. Moreover, on over 73% tasks, H* NLL is at least 7.74 times faster than random search; and on over 75% tasks, H* KL is at least 6.07 times faster than random search.

heading: E.2.2 EFFECT OF NUMBER OF TRAINING TASKS
text: We now investigate the impact of number of training tasks on the performance of meta BO methods. In Fig 8 we show the BO simple regrets on tasks from Table 1 (except ImageNet ResNet50 2048) that use meta BO models trained on different number of training tasks. To analyze the performance of all methods on less-related tasks, we first remove training tasks that have the same task dataset as our current tuning task for testing, and then remove randomly selected training datasets from the rest.
HyperBO variants were able to reduce the simple regret as more training tasks are given. Interestingly, both H* NLL and H* KL are already slightly better than Rand and STBOH when they started off with only 3 training tasks. There are reasonable fluctuations in the results but overall the trend of regret is going down as the number of training tasks increases. MIMO also reduced regret when the number of tasks increased from 8 to 18. RFGP, however, fails to learn from training tasks possibly because it did not learn good task embeddings for GP regression models.

heading: E.2.3 EFFECT OF NUMBER OF DATA POINTS IN TRAINING TASKS
text: One remaining question is, how does Mi in §2, the number of data points in each training tasks, affect the performance of meta BO methods. We analyze the impact of Mi by removing a portion of all data that we have access to for each task. In particular, we set the percentage of remaining data to be 0.2%, 0.5%, 1%, 3%, 5%, 10%, 30%, 50%, 70%, 90%. Remaining datapoints are selected uniformly randomly, which breaks the structure of matching data. Hence we do not include H* KL in this comparison, as H* KL only makes use of matching data.
Fig. 9 shows how the simple regret changes as the fraction of training data grows. Below 10% training data, we observe clear trend that more data lead to lower regret for both H* NLL and MIMO, and relatively no change for RFGP. We also found that the performance of HyperBO (H* NLL) does not change much as the fraction of training data increases from 5% to 90%. However, MIMO and RFGP
suffers significantly from more data as the fraction of training data increases from 5% to 50%. It is not entirely clear why MIMO and RFGP have such behaviors. One conjecture is that neural network based Bayesian linear regression models may get too confident once the amount of data reaches a certain threshold. This means much less exploration if those models are used for BO.

heading: E.2.4 TRAINING ON ALL BUT ONE TASK
text: We also studied the case where meta BO approaches have access to both training tasks that do not use the same task dataset and training tasks that use the same task dataset but different model configurations. This is especially common when we do architecture search: we aim to find the best model and we are tuning the optimizer hyperparameters for a new machine learning model given tuning data of the same task dataset on some other models.
For this section only, we added a new baseline, MAF: We refer to the meta BO method from Volpp et al. (2020) as MAF (Meta Acquisition Function) to avoid confusion. MAF used reinforcement learning to learn an acquisition function modeled by a neural network over a set of transfer learning tasks. All MAF results were generated using the code from Volpp et al. (2020). See App. E.6 for experimental details. As MAF takes significantly longer to run than HyperBO and other methods, we only include its results for this section.
We carried out a series of leave-one-out experiments, where we picked one task as the BO test function and let meta BO methods train on the remaining tasks. In Fig. 10, we aggregated results from all 23 tasks to show the trend of how each method performs.
The conclusions are similar to those from §E.2.1. As expected, STBO without any tricks to avoid pitfalls of vanilla BO did not show very good results. We inspected its learned GP which mimicked a
Dirac function that is flat almost everywhere except some locations, and hence it got very confident that it landed at a good spot and lost its ability to explore.
STBOH, on the other hand, achieved very competitive results. This is because it used hand-tuned priors on all of its GP parameters, although they were tuned on somewhat different problems than the ones we consider. As part of the goals of meta learning, we would like to show that it is possible for meta BO methods to exceed or at least match STBOH.
Both HyperBO variants obtained better results than the hand-tuned STBOH. Especially in the beginning few BO iterations, it is able to locate much better hyperparameters than all other methods.
Tab. 4 presents mean and standard error of the best validation error rates achieved in 100 BO iterations on the 23 tasks. HyperBO and its variants were able to achieve the best performance on 20 out of 23 tasks. In Fig. 11, we show the optimization curves of 4 individual tasks that are considered most
difficult because few similar task datasets are present in their training data. On all of these 4 difficult tasks, HyperBO identified good hyperparameters much sooner than its competitors.

heading: E.3 RESULTS ON ONLINE OPTIMIZER HYPERPARAMETER TUNING TASKS
text: Finally, we look into the online BO setting where we optimize over the full hypercube. In the online setting, some combinations of hyperparameters may be infeasible to evaluate. For example, an overly big learning rate may lead to divergence in gradients, in which case we do not obtain a valid model. To address this, we pre-process the function values to [−2, 2) such that infeasible evaluations map to −2, while bad evaluations approach asymptotically to −2. More precisely, for each subdataset Dfi , we applied for each successful y ∈ {yj(i)}Mij=1 the following mapping:
y ← softplus(y − y) softplus(ymax − y) ∗ 4− 2
where y is the median of {yj(i)}Mij=1.
Table 4: The mean and standard error of best validation error rates (%) for each test task in the offline leave-one-out experiments. We show results of the top 6 methods, and we highlight the lowest error rates in bold.
Rand STBOH MIMO MAF H* NLL H* KL WMT XFormer 64 34.27 ± 0.16 34.15 ± 0.15 34.40 ± 0.13 34.09 ± 0.09 33.91 ± 0.01 33.97 ± 0.02 Uniref50 Transformer 128 79.06 ± 0.04 78.92 ± 0.12 79.17 ± 0.13 79.34 ± 0.27 78.71 ± 0.06 78.64 ± 0.00 LM1B Transformer 2048 61.96 ± 0.03 61.95 ± 0.04 61.96 ± 0.05 62.02 ± 0.10 61.81 ± 0.01 61.81 ± 0.01 SVHN WRN 1024 3.99 ± 0.04 4.05 ± 0.10 3.83 ± 0.04 4.10 ± 0.09 4.10 ± 0.02 4.08 ± 0.01 SVHN WRN 256 3.71 ± 0.01 3.72 ± 0.02 3.65 ± 0.01 3.69 ± 0.03 3.78 ± 0.01 3.72 ± 0.03 ImageNet ResNet50 256 23.03 ± 0.07 22.66 ± 0.07 22.73 ± 0.07 26.44 ± 1.98 22.53 ± 0.00 22.58 ± 0.04 ImageNet ResNet50 512 23.02 ± 0.11 22.74 ± 0.05 23.01 ± 0.05 25.46 ± 1.41 22.65 ± 0.02 22.79 ± 0.03 MNIST CNNPoolTanh 2048 0.55 ± 0.01 0.53 ± 0.01 0.53 ± 0.01 0.52 ± 0.01 0.59 ± 0.02 0.54 ± 0.00 MNIST CNNPoolTanh 256 0.51 ± 0.01 0.48 ± 0.01 0.47 ± 0.00 0.47 ± 0.01 0.46 ± 0.01 0.47 ± 0.01 MNIST CNNPoolReLU 2048 0.69 ± 0.01 0.73 ± 0.02 0.67 ± 0.02 0.68 ± 0.01 0.64 ± 0.00 0.70 ± 0.03 MNIST CNNPoolReLU 256 0.51 ± 0.01 0.55 ± 0.03 0.50 ± 0.01 0.51 ± 0.01 0.49 ± 0.00 0.49 ± 0.00 MNIST CNNReLU 2048 1.14 ± 0.03 1.20 ± 0.09 1.10 ± 0.01 1.17 ± 0.02 1.06 ± 0.00 1.11 ± 0.02 MNIST CNNReLU 256 1.09 ± 0.02 1.06 ± 0.01 1.08 ± 0.02 1.07 ± 0.02 1.03 ± 0.00 1.04 ± 0.01 Fashion CNNPoolTanh 2048 7.14 ± 0.06 7.10 ± 0.05 7.01 ± 0.04 7.12 ± 0.04 7.00 ± 0.04 7.02 ± 0.07 Fashion CNNPoolTanh 256 6.51 ± 0.03 6.67 ± 0.18 6.40 ± 0.05 6.47 ± 0.03 6.40 ± 0.04 6.34 ± 0.04 Fashion CNNPoolReLU 2048 7.47 ± 0.02 7.48 ± 0.04 7.54 ± 0.06 7.63 ± 0.04 7.47 ± 0.03 7.47 ± 0.02 Fashion CNNPoolReLU 256 6.78 ± 0.04 6.74 ± 0.01 7.03 ± 0.07 6.84 ± 0.05 6.74 ± 0.03 6.81 ± 0.05 Fashion CNNReLU 2048 7.70 ± 0.03 7.47 ± 0.09 7.60 ± 0.04 40.40 ± 17.80 7.54 ± 0.01 7.57 ± 0.02 Fashion CNNReLU 256 7.70 ± 0.04 7.46 ± 0.11 7.84 ± 0.06 24.13 ± 14.54 7.29 ± 0.05 7.25 ± 0.05 CIFAR100 WRN 2048 21.28 ± 0.27 20.78 ± 0.19 21.75 ± 0.15 50.70 ± 15.44 21.22 ± 0.23 20.82 ± 0.19 CIFAR100 WRN 256 19.17 ± 0.19 19.02 ± 0.03 19.12 ± 0.04 19.84 ± 0.13 19.00 ± 0.00 19.04 ± 0.05 CIFAR10 WRN 2048 3.73 ± 0.05 3.43 ± 0.07 3.46 ± 0.05 3.40 ± 0.06 3.55 ± 0.10 3.43 ± 0.05 CIFAR10 WRN 256 2.84 ± 0.04 2.88 ± 0.06 2.89 ± 0.06 3.04 ± 0.05 2.82 ± 0.03 2.74 ± 0.01
In this section, we set HyperBO variants and STBO to share exactly the same GP-UCB acquisition function as STBOH, MIMO and RFGP. The UCB coefficient for all methods is 1.8. The variants of HyperBO are as follows:
• H* NLL: HyperBO with UCB as the acquisition function and negative log marginal likelihood (NLL) as the objective function.
• H* NLLKL: HyperBO with UCB as the acquisition function and NLL plus 10 times KL divergence on matching datapoints as the objective function. See §C for more details.
In Fig. 12, we include the online tuning results for selected tasks due to limited compute resources. We noticed that for some methods, e.g. STBO and MIMO, it is very difficult for them to recover from a “bad” datapoint. This is partly because predictions from these models are significantly tied to the initial observations. For example, STBO may overfit to the initial bad value and believe there are bad values in the entire search space. Nevertheless, in 7 out of 9 tasks, HyperBO methods performed the best among all methods being compared.
E.4 IMPACT OF OBJECTIVE FUNCTIONS
Here we investigate how different objective functions in HyperBO can impact its performance. Besides NLL and KL, which are already described in details in §3, we also include NLL+KL, which corresponds to Eq. 8 with λ = 10. λ = 10 is an arbitrary choice, and one may find other better options to set λ in Eq. 8.
Figure 13 shows the performance profiles and BO simple regret curves of NLL, KL and NLL+KL when HyperBO uses different acquisition functions. As comparisons, we also included the better performing baselines: Rand, STBOH and MIMO. While the ranks of performance by different objectives do vary depending on which acquisition function we used, it is clear that all HyperBO variants outperform the baselines. For EI and PI, the KL objective gives better performance, but NLL or NLL+KL might be preferred for UCB with coefficient 2, 3 or 4.
E.5 IMPACT OF ACQUISITION FUNCTIONS
As explained briefly in §4, we used 5 acquisition functions in our experiments: the vannila EI method, PI with coefficient 0.1, αPI ( x;GP(µ̂, k̂ | Df ) ) =
µ̂Df (x)−maxt(yt+0.1) σ̂Df (x) , and UCB with coefficient ζ = 2, 3 and 4 in αUCB ( x;GP(µ̂, k̂ | Df ) ) = µ̂Df (x) + ζσ̂Df (x).
Our goal is to verify that HyperBO maintains good performance across different choices of acquisition functions. To do so, we investigated in the performance of HyperBO variants under different objectives. We avoid over cluttering the figures by only including STBOH as baseline, since it is roughly the best baseline according to the main results in Fig. 6.
As shown in Figure 14, HyperBO with EI and PI as acquisition functions perform relatively better than HyperBO with UCB variants. However, HyperBO with UCB3 can still be very competitive when it is coupled with NLL objective. Overall, HyperBO with all of the 5 acquisition function options outperforms the best performing baselines.

heading: E.6 MAF IMPLEMENTATION DETAILS
text: We compared to (Volpp et al., 2020) using the code and default hyperparameters provided by the authors.7 This code assumes that each task is additionally accompanied by the optimal set of hyperparameters for the GP used to model the task (including the task used for evaluation). Following the MAF approach, we learned these hyperparameters using the GPY library, and provided them to
7https://github.com/boschresearch/MetaBO
the MAF algorithm. Given that MAF takes significantly longer to run than HyperBO, each subdataset was evaluated using only one random seed.
Each neural acquisition function was trained for a total of 1000 iterations. As was done in (Volpp et al., 2020), we selected the optimal training iteration for the neural acquisition function by cross-validation on the transfer learning tasks; in this case, we randomly sampled 3 transfer learning task, and chose the training iteration with the lowest average simple regret.
Finally, to reuse the MAF code, we also had to ensure that (a) each subtask had the same number of evaluation points, and (b) that there were no duplicated hyperparameters. For this reason, we first removed all duplicate hyperparameters within each subdataset, then capped each subdataset to the first 1559 points (the size of the smallest sub-dataset).

heading: F DISCUSSION
text: In this work, we focused on the question of how to make use of multi-task data to enable better Bayesian optimization. For the convenience of our investigation, we made simplifications such as sequential evaluations and a shared search space across tasks, although these are mostly unnecessary (see below). Our method also relies on an important assumption: functions of all tasks are i.i.d. samples from the same GP. In this section, we explore how reasonable the i.i.d. assumption is and discuss extensions to our work that would enable even more flexible uses.
Assumption on i.i.d. GP samples. To get a better idea on how much our assumptions helped on training the GP, we compare NLLs associated with 23 tasks in §4.1 with models obtained via 3 scenarios:
(a) No training: a randomly initialized model with no training;
(b) Single task: models trained on 100 randomly selected data points of the test task;
(c) H*: models trained on 18 irrelevant tasks selected in §E.2.2.
Here case (c) corresponds to the method HyperBO used for training a GP and case (b) corresponds to the model STBO can obtain with 100 initial observations. In Tab. 5, we show NLLs of these 3 methods on all tasks8 and NLLs on the test task. Note that the held-out tasks for some test tasks are the same because of the held-out rules in §E.2.1.
Comparing NLLs of the test tasks using models without training and trained via marginal likelihood like STBO, it is perhaps surprising to see that training on a subset of data points of the sub-dataset of the test task not only did not contribute to lowering NLL on the entire sub-dataset, but it even made it worse in 20 out of 23 test tasks. The training process by optimizing the NLL on a part of a sub-dataset leads to severe over-fitting. We can observe the same results of NLLs on all tasks. Without any training, our NLL is 148211.2. Yet single-task training leads to higher NLLs for all models trained on different sub-datasets.
Our method H*, on the other end, consistently achieves lower NLLs on both the test task and all tasks. Although it is not entirely clear what the relation is between a better NLL of the GP and better BO results, achieving lower NLLs typically means that the model has a better fit to the dataset. Hence, by the assumption of typical BO methods, the test function should look like a sample from our model, and so lower NLLs of model will contribute to matching the assumption of typical BO methods. By enhancing the assumption with ours on i.i.d. GP samples, Tab. 5 shows we then will be able to obtain models with a much better fit to the data.
We also computed the (pseudo) KL divergence across all matching datasets in the last columns of Tab. 5. See Appendix C for a comprehensive analysis on pseudo KL divergence for degenerate multivariate Gaussians. Note that pseudo KL divergence can be negative. Here we use pseudo KL divergence if required by the matching dataset. Again, single-task training leads to unstable (pseudo) KL values, sometimes even higher than without training (2177.2). On the contrary, training with H* leads to much more stable and lower values for KL. This indicates that the model learned to predict
8All tasks include ImageNet ResNet50 2048. But it is excluded in the test tasks in Tab. 5 because it has much fewer data points than the others.
similarly to the sample mean/covariance estimate, which is known to help better selection of BO query points by Theorem 1.
Batch evaluation. For simplicity of this paper, we did not consider batch evaluation but rather only focused on the prior selection dimension of the challenges in BO. However, it is straightforward to adopt any batch BO methods in conjunction with HyperBO to support obtaining observations in parallel. For example, we can directly use batch methods in Snoek et al. (2012), Kathuria et al. (2016), or Wang et al. (2017) to replace line 5 of Alg. 1.
High-dimensional and large scale data. Similar to batch BO, our method can also be naturally combined with most high-dimensional and large scale BO methods to offer more capabilities. For these cases, typically a probabilistic model different from vanilla GPs may be adopted. In line 2 of Alg. 1, it is straightforward to adapt our method to optimize the cumulative marginal likelihood in Eq. 3 instead for the new model. Our meta-learning idea in this paper in fact also brings benefit to high-dimensional and large scale BO methods so that they can better identify their critical special structures, e.g. low-dimensional embedding Wang et al. (2016), cylindrical kernels Oh et al. (2018) or additive Mondrian kernels Wang et al. (2018a).
Different search spaces. Roughly speaking, there could be two circumstances for difference search spaces. Case I is that tasks share the same search variables, but the search ranges for some variables are different. For example, we may have each function fi : Xi → R, i ∈ [N ] and Xi = ∏d j=1[lij , hij ] ⊂ Rd. In this case, our solution still applies by simply setting a union search
space as X = ⋃N i=1 Xi for learning and use the designated search space of new tasks for optimization.
Case II is more complicated: the search space for each function fi is Xi ⊂ Rdi and each dimension of Xi may have a different meaning than another search space Xj (i 6= j). This paper does not have a solution for this scenario. Further research will be needed to reduce Case II to Case I which can be then immediately combined with HyperBO.

Review summaries:
This paper is concerned with speeding up Bayesian optimization by using evaluation data from previous, related tasks defined over the same configuration space. The authors propose to model the data from each experiment (or "task") by independent Gaussian processes, which all share the same mean and covariance function. This surrogate model can be learned from past data.  The paper also presents experiments on a fairly simple search space of 4 optimizer parameters. This is done for a bunch of datasets and NN models. And there is a pretty simple extension of theoretical results from (Wang, 2018b).

This paper suggests a meta Bayesian optimization strategy that optimizes free parameters of GP including a prior function and noise variance, where multiple sets of historical observations are given. In particular, the proposed method chooses a free parameters using one of three approaches: (i) optimizing a marginal likelihood, (ii) measuring KL divergence, (iii) considering both marginal likelihood and KL divergence. The authors finally show the theoretical analyses on regret bounds and the numerical results on hyperparameter optimization.

 HyperBO assumes the tasks are independent given the hyperparameters, unlike typical metalearning approaches which assume tasks are related. This allows for an efficient Kronecker decomposition of the kernel and thus linear, rather than cubic scaling, across tasks.   Using this model, HyperBO performs BO as usual; maximize the acquisition function to obtain the next point to evaluate. HyperBO also makes the critical assumption of an offline pre-training of hyperparameters on a representative set of completed tasks; during optimization itself the hyperparameters are fixed.  

This paper presents a Bayesian optimization method based on meta-BO. The motivation is tasks can share the same parameter structure and this shared information, e.g. correlation between tasks, can be transferred to new and similar tasks. An example is to optimize the the hyper-parameters of a same optimizer across different architectures and different datasets. This problem is a very important one in the community of Bayesian optimization and a reasonable method can lead to a potentially dramatic decrease in the required computation, especially when the objective function is very expensive. This work tries to overcome limitations of existing methods. For example, the method proposed in this work does not need to evaluate all objective functions associated with all tasks on the same parameters.