# Hyperparameter Optimization: A Spectral Approach

We give a simple, fast algorithm for hyperparameter optimization inspired by techniques from the analysis of Boolean functions. We focus on the high-dimensional regime where the canonical example is training a neural network with a large number of hyperparameters. The algorithm --- an iterative application of compressed sensing techniques for orthogonal polynomials --- requires only uniform sampling of the hyperparameters and is thus easily parallelizable. Experiments for training deep neural networks on Cifar-10 show that compared to state-of-the-art tools (e.g., Hyperband and Spearmint), our algorithm finds significantly improved solutions, in some cases better than what is attainable by hand-tuning. In terms of overall running time (i.e., time required to sample various settings of hyperparameters plus additional computation time), we are at least an order of magnitude faster than Hyperband and Bayesian Optimization. We also outperform Random Search 8x. Additionally, our method comes with provable guarantees and yields the first improvements on the sample complexity of learning decision trees in over two decades. In particular, we obtain the first quasi-polynomial time algorithm for learning noisy decision trees with polynomial sample complexity.

• 55 publications
• 13 publications
• 17 publications
04/25/2016

### CMA-ES for Hyperparameter Optimization of Deep Neural Networks

Hyperparameters of deep neural networks are often optimized by grid sear...
02/23/2020

### Weighting Is Worth the Wait: Bayesian Optimization with Importance Sampling

Many contemporary machine learning models require extensive tuning of hy...
09/20/2019

### Bayesian Optimization for Iterative Learning

The success of deep (reinforcement) learning systems crucially depends o...
12/05/2018

### An empirical study on hyperparameter tuning of decision trees

Machine learning algorithms often contain many hyperparameters whose val...
06/29/2022

### Open Problem: Properly learning decision trees in polynomial time?

The authors recently gave an n^O(loglog n) time membership query algorit...
09/28/2021

### Faster Improvement Rate Population Based Training

The successful training of neural networks typically involves careful an...
03/29/2021

### One Network Fits All? Modular versus Monolithic Task Formulations in Neural Networks

Can deep learning solve multiple tasks simultaneously, even when they ar...

## 1 Introduction

Large scale machine learning and optimization systems usually involve a large number of free parameters for the user to fix according to their application. A timely example is the training of deep neural networks for a signal processing application: the ML specialist needs to decide on an architecture, depth of the network, choice of connectivity per layer (convolutional, fully-connected, etc.), choice of optimization algorithm and recursively choice of parameters inside the optimization library itself (learning rate, momentum, etc.).

Given a set of hyperparameters and their potential assignments, the naive practice is to search through the entire grid of parameter assignments and pick the one that performed the best, a.k.a. “grid search”. As the number of hyperparameters increases, the number of possible assignments increases exponentially and a grid search becomes quickly infeasible. It is thus crucial to find a method for automatic tuning of these parameters.

This auto-tuning, or finding a good setting of these parameters, is now referred to as hyperparameter optimization (HPO), or simply automatic machine learning (auto-ML). For continuous hyperparameters, gradient descent is usually the method of choice [MDA15, LBGR15, FLF16]. Discrete parameters, however, such as choice of architecture, number of layers, connectivity and so forth are significantly more challenging. More formally, let

 f:{−1,1}n↦[0,1]

be a function mapping hyperparameter choices to test error of our model. That is, each dimension corresponds to a certain hyperparameter (number of layers, connectivity, etc.), and for simplicity of illustration we encode the choices for each parameter as binary numbers . The goal of HPO is to approximate the minimizer in the following setting:

1. Oracle model: evaluation of for a given choice of hyperparameters is assumed to be very expensive. Such is the case of training a given architecture of a huge dataset.

2. Parallelism is crucial: testing several model hyperparameters in parallel is entirely possible in cloud architecture, and dramatically reduces overall optimization time.

3. is structured.

The third point is very important since clearly HPO is information-theoretically hard and evaluations of the function are necessary in the worst case. Different works have considered exploiting one or more of the properties above. The approach of Bayesian optimization [SLA12] addresses the structure of , and assumes that a useful prior distribution over the structure of is known in advance. Multi-armed bandit algorithms [LJD16], and Random Search [BB12], exploit computational parallelism very well, but do not exploit any particular structure of . These approaches are surveyed in more detail later.

### 1.1 Our contribution

In this paper we introduce a new spectral approach to hyperparameter optimization. Our main idea is to make assumptions on the structure of in the Fourier domain. Specifically we assume that can be approximated by a sparse and low degree polynomial in the Fourier basis. This means intuitively that it can be approximated well by a decision tree.

The implication of this assumption is that we can obtain a rigorous theoretical guarantee: approximate minimization of over the boolean hypercube with function evaluations only linear in sparsity that can be carried out in parallel

. We further give improved heuristics on this basic construction and show experiments showing our assumptions are validated in practice for HPO as applied to deep learning over image datasets.

Thus our contributions can be listed as:

• A new spectral method called Harmonica that has provable guarantees: sample-efficient recovery if the underlying hyperparameter objective is a sparse (noisy) polynomial and easy to implement on parallel architectures.

• Improved bounds on the sample complexity of learning noisy, size decision trees over

variables under the uniform distribution. We observe that the classical sample complexity bound of

due to Linial et al. [LMN93] to quadratic in the size of the tree while matching the best known quasipolynomial bound in running time.

• We demonstrate significant improvements in accuracy, sample complexity, and running time for deep neural net training experiments. We compare ourselves to state-of-the-art solvers from Bayesian optimization, Multi-armed bandit techniques, and Random Search. Projecting to even higher numbers of hyperparameters, we perform simulations that show several orders-of-magnitude of speedup versus Bayesian optimization techniques.

### 1.2 Previous work

The literature on discrete-domain HPO can be roughly divided into two: probabilistic approaches and decision-theoretic methods. In critical applications, researchers usually use a grid search over all parameter space, but that becomes quickly prohibitive as the number of hyperparameter grows. Gradient-based methods such as [MDA15, LBGR15, FLF16, Ben00] are applicable only to continuous hyperparameters which we do not consider.

#### Probabilistic methods and Bayesian optimization.

Bayesian optimization (BO) algorithms [BBBK11, SLA12, SSA13, SSZA14, GKX14, WZH13, IAFS17]

tune hyperparameters by assuming a prior distribution of the loss function, and then keep updating this prior distribution based on the new observations. Each new observation is selected according to an acquisition function, which balances exploration and exploitation such that the new observation gives us a better result, or helps gain more information about the loss function. The BO approach is inherently serial and difficult to parallelize, and its theoretical guarantees have thus far been limited to statistical consistency (convergence in the limit).

#### Decision-theoretic methods.

Perhaps the simplest approach to HPO is random sampling of different choices of parameters and picking the best amongst the chosen evaluations [BB12]. It is naturally very easy to implement and parallelize. Upon this simple technique, researchers have tried to allocate different budgets to the different evaluations, depending on their early performance. Using adaptive resource allocation techniques found in the multi-armed bandit literature, Successive Halving (SH) algorithm was introduced [JT16]. Hyperband further improves SH by automatically tuning the hyperparameters in SH [LJD16].

#### Learning decision trees.

Prior work for learning decision trees (more generally Boolean functions that are approximated by low-degree polynomials) used the celebrated “low-degree” algorithm of Linial, Mansour, and Nisan [LMN93]

. Their algorithm uses random sampling to estimate each low-degree Fourier coefficient to high accuracy.

We make use of the approach of Stobbe and Krause [SK12], who showed how to learn low-degree, sparse Boolean functions using tools from compressed sensing (similar approaches were taken by Kocaoglu et al. [KSDK14] and Negahban and Shah [NS12]). We observe that their approach can be extended to learn functions that are both “approximately sparse” (in the sense that the norm of the coefficients is bounded) and “approximately low-degree” (in the sense that most of the mass of the Fourier spectrum resides on monomials of low-degree). This implies the first decision tree learning algorithm with polynomial sample complexity that handles adversarial noise. In addition, we obtain the optimal dependence on the error parameter .

For the problem of learning exactly -sparse Boolean functions over variables, Haviv and Regev [HR15] have recently shown that uniformly random samples suffice. Their result is not algorithmic but does provide an upper bound on the information-theoretic problem of how many samples are required to learn. The best algorithm in terms of running time for learning -sparse Boolean functions is due to [FGKP09], and requires time . It is based on the Blum, Kalai, and Wasserman algorithm for learning parities with noise [BKW03].

#### Techniques.

Our methods are heavily based on known results from the analysis of boolean functions as well as compressed sensing. The relevant material and literature are given in the next section.

## 2 Setup and definitions

The problem of hyperparameter optimization is that of minimizing a discrete, real-valued function, which we denote by (we can handle arbitrary inputs, binary is chosen for simplicity of presentation).

In the context of hyperparameter optimization, function evaluation is very expensive, although parallelizable, as it corresponds to training a deep neural net. In contrast, any computation that does not involve function evaluation is considered less expensive, such as computations that require time for “somewhat large” or are subexponential (we still consider runtimes that are exponential in to be costly).

### 2.1 Basics of Fourier analysis

The reader is referred to [O’D14] for an in depth treatment of Fourier analysis of Boolean functions. Let be a function over domain . Let

. We write and say that are -close if

 Ex∼D[(f(x)−g(x))2] ≤ ε.
###### Definition 1.

[Rau10] We say a family of functions ( maps to ) is a Random Orthonormal Family with respect to if

 ED[ψi(X)⋅ψj(X)]=δij={1 if i=j0 otherwise.

The expectation is taken with respect to probability distribution . We say that the family is -bounded if for every . Henceforth we assume .

An important example of a random orthonormal family is the class of parity functions with respect to the uniform distribution on :

###### Definition 2.

A parity function on some subset of variables is the function where .

It is easy to see that the set of all parity functions , one for each , form a random orthonormal family with respect to the uniform distribution on .

This random orthonormal family is often referred to as the Fourier basis, as it is a complete orthonormal basis for the class of Boolean functions with respect to the uniform distribution on . More generally, for any , can be uniquely represented in this basis as

 f(x)=∑S⊆[n]^fSχS(x)

where

 ^fS=⟨f,χS⟩=Ex∈{−1,1}n[f(x)χS(x)]

is the Fourier coefficient corresponding to where is drawn uniformly from . We also have Parseval’s identity: .

In this paper, we will work exclusively with the above parity basis. Our results apply more generally, however, to any orthogonal family of polynomials (and corresponding product measure on ). For example, if we wished to work with continuous hyperparameters, we could work with families of Hermite orthogonal polynomials with respect to multivariate spherical Gaussian distributions.

We conclude with a definition of low-degree, approximately sparse (bounded norm) functions:

###### Definition 3 (Approximately sparse function).

Let be the parity basis, and let be a class of functions mapping to . Thus for , . We say that:

• A function is -sparse if , ie., f has at most nonzero entries in its polynomial expansion.

• is -concentrated if .

• is -bounded if for every , is -concentrated and in addition has norm bounded by , that is, for every we have .

It is easy to see that the class of functions with bounded norm is more general than sparse functions. For example, the Boolean function has norm bounded by but is not sparse.

We also have the following simple fact:

###### Fact 4.

[Man94] Let be such that . Then there exists such that is sparse and . The function is constructed by taking all coefficients of magnitude or larger in ’s expansion as a polynomial.

### 2.2 Compressed sensing and sparse recovery

In the problem of sparse recovery

, a learner attempts to recover a sparse vector

which is sparse, i.e. , from an observation vector that is assumed to equal

 y=Ax+e,

where is assumed to be zero-mean, usually Gaussian, noise. The seminal work of [CRT06, Don06] shows how can be recovered exactly under various conditions on the observation matrix and the noise. The usual method for recovering the signal proceeds by solving a convex optimization problem consisting of minimization as follows (for some parameter ):

 minx∈Rn{∥x∥1+λ∥Ax−y∥22}. (1)

The above formulation comes in many equivalent forms (e.g., Lasso), where one of the objective parts may appear as a hard constraint.

For our work, the most relevant extension of traditional sparse recovery is due to Rauhut [Rau10], who considers the problem of sparse recovery when the measurements are evaluated according to a random orthonormal family. More concretely, fix with non-zero entries. For -bounded random orthonormal family }, and independent draws from corresponding distribution define the matrix such that . Rauhut gives the following result for recovering sparse vectors :

###### Theorem 5 (Sparse Recovery for Random Orthonormal Families, [Rau10] Theorem 4.4).

Given as input matrix and vector with for some vector with , mathematical program (1) finds a vector such that (for constants and )

 ∥x−x∗∥2 ≤ c1σs(x)1√s+c2η

with probability as long as, for sufficiently large constant C,

 m ≥ CK2logK⋅slog3s⋅log2N⋅log(1/δ).

The term is equal to . Recent work [Bou14, HR16] has improved the dependence on the polylog factors in the lower bound for .

## 3 Basic Algorithm and Main Theoretical Results

The main component of our spectral algorithm for hyperparameter optimization is given in Algorithm 1. It is essentially an extension of sparse recovery (basis pursuit or Lasso) to the orthogonal basis of polynomials in addition to an optimization step. See Figure 1 for an illustration. We prove Harmonica’s theoretical guarantee, and show how it gives rise to new theoretical results in learning from the uniform distribution.

In the next section we describe extensions of this basic algorithm to a more practical algorithm with various heuristics to improve its performance.

###### Theorem 6 (Noiseless recovery).

Let be a -bounded orthonormal polynomial basis for distribution . Let be a -bounded function as per definition 3 with respect to the basis . Then Algorithm 1, in time and sample complexity , returns such that

 x⋆∈argminf(x)

This theorem, and indeed most of the theoretical results of this paper, follow from the main recovery properties of Procedure 2. Our main technical lemma follows the same outline of the compressed sensing result due to Stobbe and Krause [SK12] but with a generalization to functions that are approximately sparse and low-degree:

###### Lemma 7 (Noisy recovery).

Let be a -bounded orthonormal polynomial basis for distribution . Let be a -bounded as per definition 3 with respect to the basis . Then Procedure 2 finds a function in time and sample complexity .

In the rest of this section we proceed to prove the main lemma and derive the theorem. Recall the Chebyshev inequality:

###### Fact 8 (Multidimensional Chebyshev inequality).

Let be an dimensional random vector, with expected value , and covariance matrix .

If is a positive definite matrix, for any real number :

 P(√(X−μ)TV−1(X−μ)>δ) ≤ mδ2
###### Proof of Lemma 7.

For ease of notation we assume . Let be an -bounded function written in the orthonormal basis as . We can equivalently write as , where is a degree polynomial that only includes coefficients of magnitude at least and the constant term of the polynomial expansion of .

Since , by Fact 4 we have that is sparse. The function is thus the sum of the remaining terms not included in .

Draw (to be chosen later) random labeled examples and enumerate all basis functions for as . Form matrix such that and consider the problem of recovering sparse given where is the vector of coefficients of , the th entry of equals , and .

We will prove that with constant probability over the choice random examples, . Applying Theorem 5 by setting and observing that , we will recover such that for some constant . As such, for the function we will have by Parseval’s identity. Note, however, that we may rescale by constant factor to obtain error and only incur an additional constant (multiplicative) factor in the sample complexity bound.

By the definition of , we have

 ∥g∥2=⎛⎝∑S,|S|>d^f(S)2+∑R^f(R)2⎞⎠ (3)

where each is of magnitude at most . By Fact 4 and Parseval’s identity we have . Since is -concentrated we have . Thus, is at most . Therefore, by triangle inequality .

It remains to bound . Note that since the examples are chosen independently, the entries

are independent random variables. Since

is a linear combination of orthonormal monomials (not including the constant term), we have

. Here we can apply linearity of variance (the covariance of

and is zero for all ) and calculate the variance

 Var(g(zi))=(∑S,|S|>d^f(S)2+∑R^f(R)2)

With the same calculation as (3), we know is at most .

Now consider the covariance matrix of the vector which equals (recall every entry of has mean ). Then is a diagonal matrix (covariance between two independent samples is zero), and every diagonal entry is at most . Applying Fact 8 we have

 P(∥e∥2>√ε2δ) ≤ mδ2.

Setting , we conclude that . Hence with probability at least , we have that . From Theorem 5, we may choose . This completes the proof. Note that the probability above can be boosted to any constant probability with a constant factor loss in sample complexity. ∎

Remark: Note that the above proof also holds in the adversarial or agnostic noise setting. That is, an adversary could add a noise vector to the labels received by the learner. In this case, the learner will see label vector . If , then we will recover a polynomial with squared-error at most via re-scaling by a constant factor and applying the triangle inequality to .

While this noisy recovery lemma is the basis for our enhanced algorithm in the next section as well as the learning-theoretic result on learning of decision trees detailed in the next subsection, it does not imply recovery of the global optimum. The reason is that noisy recovery guarantees that we output a hypothesis close to the underlying function, but even a single noisy point can completely change the optimum.

Nevertheless, we can use our techniques to prove recovery of optimality for functions that are computed exactly by a sparse, low-degree polynomial.

###### Proof of Theorem 6.

There are at most polynomials with . Let the enumeration of these polynomials be . Draw labeled examples independently from and construct an matrix with . Since can be written as an sparse linear combination of , there exists an -sparse vector such that where the th entry of is . Hence we can apply Theorem 5 to recover exactly. These are the non-zero coefficients of ’s expansion in terms of . Since is recovered exactly, its minimizer is found in the optimization step. ∎

### 3.1 Application: Learning Decision Trees in Quasi-polynomial Time and Polynomial Sample Complexity

Here we observe that our results imply new bounds for decision-tree learning. For example, we obtain the first quasi-polynomial time algorithm for learning decision trees with respect to the uniform distribution on with polynomial sample complexity and an optimal dependence on the error parameter :

###### Corollary 9.

Let and let be the class of all decision trees of size on variables. Then is learnable with respect to the uniform distribution in time and sample complexity . Further, if the labels are corrupted by arbitrary noise vector such that

, then the output classifier will have squared-error at most

.

###### Proof.

As mentioned earlier, the orthonormal polynomial basis for the class of Boolean functions with respect to the uniform distribution on is the class of parity functions for . Further, it is easy to show that for Boolean function , if then . The corollary now follows by applying Lemma 7 and two known structural facts about decision trees: 1) a tree of size is -concentrated and has norm bounded by (see e.g., Mansour [Man94]) and 2) by Fact 4, for any function with norm bounded by (i.e., a decision tree of size ), there exists an sparse function such that . The noise tolerance property follows immediately from the remark after the proof of Lemma 7. ∎

Comparison with the “Low-Degree” Algorithm. Prior work for learning decision trees (more generally Boolean functions that are approximated by low-degree polynomials) used the celebrated “low-degree” algorithm of Linial, Mansour, and Nisan [LMN93]. Their algorithm uses random sampling to estimate each low-degree Fourier coefficient to high accuracy. In contrast, the compressed-sensing approach of Stobbe and Krause [SK12] takes advantage of the incoherence of the design matrix and gives results that seem unattainable from the “low-degree” algorithm.

For learning noiseless, Boolean decision trees, the low-degree algorithm uses quasipolynomial time and sample complexity to learn to accuracy . It is not clear, however, how to obtain any noise tolerance from their approach.

For general real-valued decision trees where is an upper bound on the maximum value at any leaf of a size tree, our algorithm will succeed with sample complexity and be tolerant to noise while the low-degree algorithm will use (and will have no noise tolerance properties). Note the improvement in the dependence on (even in the noiseless setting), which is a consequence of the RIP property of the random orthonormal family.

## 4 Harmonica: The Full Algorithm

Rather than applying Algorithm 1 directly, we found that performance is greatly enhanced by iteratively using Procedure 2 to estimate the most influential hyperparameters and their optimal values.

In the rest of this section we describe this iterative heuristic, which essentially runs Algorithm 1 for multiple stages. More concretely, we continue to invoke the PSR subroutine until the search space becomes small enough for us to use a “base” hyperparameter optimizer (in our case either SH or Random Search).

The space of minimizing assignments to a multivariate polynomial is a highly non-convex set that may contain many distinct points. As such, we take an average of several of the best minimizers (of subsets of hyperparameters) during each stage.

In order to describe this formally we need the following definition of a restriction of function:

###### Definition 10 (restriction [O’d14]).

Let , , and be given. We call a restriction pair of function . We denote the function over variables given by setting the variables of to .

We can now describe our main algorithm (Algorithm 3). Here is the number of stages for which we apply the PSR subroutine, and the restriction size serves as a tie-breaking rule for the best minimizers (which can be set to ).

### 4.1 Algorithm attributes and heuristics

#### Scalability.

If the hidden function if -sparse, Harmonica can find such a sparse function using samples. If at every stage of Harmonica, the target function can be approximated by an sparse function, we only need samples where is the number of stages. For real world applications such as deep neural network hyperparameter tuning, it seems (empirically) reasonable to assume that the hidden function is indeed sparse at every stage (see Section 5).

For Hyperband [LJD16], SH [JT16] or Random Search, even if the function is -sparse, in order to cover the optimal configuration by random sampling, we need samples.

#### Optimization time.

Harmonica runs the Lasso [Tib96] algorithm after each stage to solve (2), which is a well studied convex optimization problem and has very fast implementations. Hyperband and SH are also efficient in terms of running time as a function of the number of function evaluations, and require sorting or other simple computations. The running time of Bayesian optimization is cubic in number of function evaluations, which limits applicability for large number of evaluations / high dimensionality, as we shall see in Section 5.4.

#### Parallelizability.

Harmonica, similar to Hyperband, SH, and Random Search, has straightforward parallel implementations. In every stage of those algorithms, we could simply evaluate the objective functions over randomly chosen points in parallel.

It is hard to run Bayesian optimization algorithm in parallel due to its inherent serial nature. Previous works explored variants in which multiple points are evaluated at the same time in parallel [WF16], though speed ups do not grow linearly in the number of machines, and the batch size is usually limited to a small number.

#### Feature Extraction.

Harmonica is able to extract important features with weights in each stages, which automatically sorts all the features according to their importance. See Section A.2.

## 5 Experiments with training deep networks

We compare Harmonica222A python implementation of Harmonica can be found at https://github.com/callowbird/Harmonica with Spearmint [SLA12], Hyperband, SH444We implemented a parallel version of Hyperband and SH in Lua. and Random Search. Both Spearmint and Hyperband are state-of-the-art algorithms, and it is observed that Random Search 2x (Random Search with doubled function evaluation resources) is a very competitive benchmark that beats many algorithms555E.g., see [Rec16a, Rec16b]. .

Our first experiment is over training residual network on Cifar-10 dataset. We included binary hyperparameters, including initialization, optimization method, learning rate schedule, momentum rate, etc. Table 2 (Section A.1) details the hyperparameters considered. We also include dummy variables to make the task more challenging. Notice that Hyperband, SH, and Random Search are agnostic to the dummy variables in the sense that they just set the value of dummy variables randomly, therefore select essentially the same set of configurations with or without the dummy variables. Only Harmonica and Spearmint are sensitive to the dummy variables as they try to learn the high dimensional function space. To make a fair comparison, we run Spearmint without any dummy variables.

As most hyperparameters have a consistent effect as the network becomes deeper, a common hand-tuning strategy is “tune on small network, then apply the knowledge to big network” (See discussion in Section A.3

). Harmonica can also exploit this strategy as it selects important features stage-by-stage. More specifically, during the feature selection stages, we run Harmonica for tuning an

layer neural network with

training epochs. At each stage, we take

samples to extract important features, and set restriction size (see Procedure 2). After that, we fix all the important features, and run the SH or Random Search as our base algorithm on the big layer neural network for training the whole epochs777Other algorithms like Spearmint, Hyperband, etc. can be used as the base algorithms as well. . To clarify, “stage” means the stages of the hyperparameter algorithms, while “epoch” means the epochs for training the neural network.

### 5.1 Performance

We tried three versions of Harmonica for this experiment, Harmonica with 1 stage (Harmonica-1), 2 stages (Harmonica-2) and 3 stages (Harmonica-3). All of them use SH as the base algorithm. The top test error results and running times of the different algorithms are depicted in Figure 2. SH based algorithms may return fewer than 10 results. For more runs of variants of Harmonica and its resulting test error, see Figure 3 (the results are similar to Figure 2).

#### Test error and scalability:

Harmonica-1 uses less than time of Hyperband and time compared with Random Search, but gets better results than the competing algorithms. It beats the Random Search x benchmark (stronger than Random Search 2x benchmark of [LJD16]). Harmonica-2 uses slightly more time, but is able to find better results, which are comparable with Spearmint with running time.

#### Improving upon human-tuned parameters:

Harmonica-3 obtains a better test error () as compared to the best hand-tuning rate reported in [HZRS16]888 is the rate obtained by residual network, and there are new network structures like wide residual network [ZK16] or densenet [HLW16] that achieve better rates for Cifar-10.. Harmonica-3 uses only 6.1 GPU days, which is less than half day in our environment, as we have 20 GPUs running in parallel. Notice that we did not cherry pick the results for Harmonica-3. In Section 5.3 we show that by running Harmonica-3 for longer time, one can obtain a few other solutions better than hand tuning.

#### Performance of provable methods:

Harmonica-1 has noiseless and noisy recovery guarantees (Lemma 7), which are validated experimentally.

### 5.2 Average Test Error For Each Stage

We computed the average test error among random samples for an layer network with epochs after each stage. See Figure 4. After selecting features in stage , the average test error drops from to , which indicates the top features are very important. As we proceed to stage , the improvement on test error becomes less significant as the selected features at stage have mild contributions.

### 5.3 Hyperparameters for Harmonica

To be clear, Harmonica itself has six hyperparameters that one needs to set including the number of stages, regularizer for Lasso, the number of features selected per stage, base algorithm, small network configuration, and the number of samples per stage. Note, however, that we have reduced the search space of general hyperparameter optimization down to a set of only six hyperparameters. Empirically, our algorithm is robust to different settings of these parameters, and we did not even attempt to tune some of them (e.g., small network configuration).

#### Base algorithm and #stages.

We tried different versions of Harmonica, including Harmonica with 1 stage, 2 stages and 3 stages using SH as the base algorithm (Harmonica-1, Harmonica-2, Harmonica-3), with 1 stage and 2 stages using Random Search as the base algorithm (Harmonica-1-Random-Search, Harmonica-2-Random-Search), and with 2 stages and 3 stages running SH as the base for longer time (Harmonica-2-Long, Harmonica-3-Long). As can be seen in Figure 3, most variants produce better results than SH and use less running time. Moreover, if we run Harmonica for longer time, we will obtain more stable solutions with less variance in test error.

#### Lasso parameters are stable.

See Table 1 for stable range for regularization term and the number of samples. Here stable range means as long as the parameters are set in this range, the top features and the signs of their weights (which are what we need for computing in Procedure 2) do not change. In other words, the feature selection outcome is not affected. When parameters are outside the stable ranges, usually the top features are still unchanged, and we miss only one or two out of the five features.

#### On the degree of features.

We set degree to be three because it does not find any important features with degree larger than this. Since Lasso can be solved efficiently (less than minutes in our experiments), the choice of degree can be decided automatically.

### 5.4 Experiments with Synthetic functions

Our second experiment considers a synthetic hierarchically bounded function . We run Harmonica with samples, features selected per stage, for stages, using degree features. See Figure 5 for optimization time comparison. We only plot the optimization time for Spearmint when , which takes more than one day for samples. Harmonica is several magnitudes faster than Spearmint. In Figure 6, we show that Harmonica is able to estimate the hidden function with error proportional to the noise level.

The synthetic function is defined as follows. has three stages, and in -th stage (), it has sparse vectors for . Each contains pairs of weight and feature for , where . and is a monomial on with degree at most . Therefore, for input , the sparse vector . Since , is binary. Therefore, contains binaries and represents a integer in , denoted as . Let , where is the noise uniformly sampled from ( is the noise level). In other words, in every stage we will get a sparse vector . Based on , we pick a the next sparse function and proceed to the next stage.

## 6 Acknowledgements

Thanks to Vitaly Feldman for pointing out the work of Stobbe and Krause [SK12]. We thank Sanjeev Arora for helpful discussions and encouragement. Elad Hazan is supported by NSF grant 1523815. This project is supported by a Microsoft Azure research award and Amazon AWS research award.

## References

• [BB12] James Bergstra and Yoshua Bengio. Random search for hyper-parameter optimization. J. Mach. Learn. Res., 13:281–305, February 2012.
• [BBBK11] James S. Bergstra, Rémi Bardenet, Yoshua Bengio, and Balázs Kégl. Algorithms for hyper-parameter optimization. In J. Shawe-Taylor, R. S. Zemel, P. L. Bartlett, F. Pereira, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 24, pages 2546–2554. Curran Associates, Inc., 2011.
• [Ben00] Yoshua Bengio. Gradient-based optimization of hyperparameters. Neural Computation, 12(8):1889–1900, 2000.
• [BKW03] Avrim Blum, Adam Kalai, and Hal Wasserman. Noise-tolerant learning, the parity problem, and the statistical query model. J. ACM, 50(4):506–519, July 2003.
• [Bou14] Jean Bourgain. An Improved Estimate in the Restricted Isometry Problem, pages 65–70. Springer International Publishing, Cham, 2014.
• [CRT06] E. J. Candes, J. Romberg, and T. Tao. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans. Inf. Theor., 52(2):489–509, February 2006.
• [Don06] D. L. Donoho. Compressed sensing. IEEE Trans. Inf. Theor., 52(4):1289–1306, April 2006.
• [FGKP09] V. Feldman, Parikshit Gopalan, Subhash Khot, and Ashok Kumar Ponnuswami. On agnostic learning parities, monomials,and halfspaces. SIAM Journal on Computing, 39(2):606–645, 2009.
• [FLF16] Jie Fu, Hongyin Luo, Jiashi Feng, Kian Hsiang Low, and Tat-Seng Chua. Drmad: Distilling reverse-mode automatic differentiation for optimizing hyperparameters of deep neural networks. CoRR, abs/1601.00917, 2016.
• [GB10] Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In

Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, AISTATS 2010, Chia Laguna Resort, Sardinia, Italy, May 13-15, 2010

, pages 249–256, 2010.
• [GKX14] Jacob R. Gardner, Matt J. Kusner, Zhixiang Eddie Xu, Kilian Q. Weinberger, and John P. Cunningham. Bayesian optimization with inequality constraints. In Proceedings of the 31th International Conference on Machine Learning, ICML 2014, Beijing, China, 21-26 June 2014, pages 937–945, 2014.
• [HLW16] Gao Huang, Zhuang Liu, and Kilian Q. Weinberger. Densely connected convolutional networks. CoRR, abs/1608.06993, 2016.
• [HR15] Ishay Haviv and Oded Regev. The list-decoding size of fourier-sparse boolean functions. In David Zuckerman, editor, 30th Conference on Computational Complexity, CCC 2015, June 17-19, 2015, Portland, Oregon, USA, volume 33 of LIPIcs, pages 58–71. Schloss Dagstuhl - Leibniz-Zentrum fuer Informatik, 2015.
• [HR16] Ishay Haviv and Oded Regev. The restricted isometry property of subsampled fourier matrices. In Proceedings of the Twenty-seventh Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’16, pages 288–297, Philadelphia, PA, USA, 2016. Society for Industrial and Applied Mathematics.
• [HZRS15] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.

Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.

In

2015 IEEE International Conference on Computer Vision, ICCV 2015, Santiago, Chile, December 7-13, 2015

, pages 1026–1034, 2015.
• [HZRS16] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In

2016 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2016, Las Vegas, NV, USA, June 27-30, 2016

, pages 770–778, 2016.
• [IAFS17] Ilija Ilievski, Taimoor Akhtar, Jiashi Feng, and Christine Shoemaker. Efficient hyperparameter optimization for deep learning algorithms using deterministic rbf surrogates. 2017.
• [IS15] Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In Proceedings of the 32nd International Conference on Machine Learning, ICML 2015, Lille, France, 6-11 July 2015, pages 448–456, 2015.
• [JT16] Kevin G. Jamieson and Ameet Talwalkar. Non-stochastic best arm identification and hyperparameter optimization. In Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, AISTATS 2016, Cadiz, Spain, May 9-11, 2016, pages 240–248, 2016.
• [KB14] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. CoRR, abs/1412.6980, 2014.
• [KSDK14] Murat Kocaoglu, Karthikeyan Shanmugam, Alexandros G. Dimakis, and Adam R. Klivans. Sparse polynomial learning and graph sketching. In Zoubin Ghahramani, Max Welling, Corinna Cortes, Neil D. Lawrence, and Kilian Q. Weinberger, editors, Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014, Montreal, Quebec, Canada, pages 3122–3130, 2014.
• [LBGR15] Jelena Luketina, Mathias Berglund, Klaus Greff, and Tapani Raiko. Scalable gradient-based tuning of continuous regularization hyperparameters. CoRR, abs/1511.06727, 2015.
• [LJD16] L. Li, K. Jamieson, G. DeSalvo, A. Rostamizadeh, and A. Talwalkar. Hyperband: A Novel Bandit-Based Approach to Hyperparameter Optimization. ArXiv e-prints, March 2016.
• [LMN93] Nathan Linial, Yishay Mansour, and Noam Nisan. Constant depth circuits, fourier transform, and learnability. J. ACM, 40(3):607–620, July 1993.
• [Man94] Yishay Mansour. Learning Boolean Functions via the Fourier Transform, pages 391–424. Springer US, Boston, MA, 1994.
• [MDA15] Dougal Maclaurin, David Duvenaud, and Ryan P. Adams. Gradient-based hyperparameter optimization through reversible learning. In Proceedings of the 32Nd International Conference on International Conference on Machine Learning - Volume 37, ICML’15, pages 2113–2122. JMLR.org, 2015.
• [NS12] Sahand Negahban and Devavrat Shah. Learning sparse boolean polynomials. In Allerton, pages 2032–2036. IEEE, 2012.
• [O’D14] Ryan O’Donnell. Analysis of Boolean Functions. Cambridge University Press, New York, NY, USA, 2014.
• [Rau10] Holger Rauhut. Compressive sensing and structured random matrices. Theoretical foundations and numerical methods for sparse recovery, 9:1–92, 2010.
• [Rec16a] Benjamin Recht. Embracing the random.
• [Rec16b] Benjamin Recht. The news on auto-tuning.
• [SHK14] Nitish Srivastava, Geoffrey E. Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. Dropout: a simple way to prevent neural networks from overfitting. Journal of Machine Learning Research, 15(1):1929–1958, 2014.
• [SK12] Peter Stobbe and Andreas Krause. Learning fourier sparse set functions. In Neil D. Lawrence and Mark A. Girolami, editors, Proceedings of the Fifteenth International Conference on Artificial Intelligence and Statistics, AISTATS 2012, La Palma, Canary Islands, April 21-23, 2012, volume 22 of JMLR Proceedings, pages 1125–1133. JMLR.org, 2012.
• [SLA12] Jasper Snoek, Hugo Larochelle, and Ryan P. Adams. Practical bayesian optimization of machine learning algorithms. In Advances in Neural Information Processing Systems 25: 26th Annual Conference on Neural Information Processing Systems 2012. Proceedings of a meeting held December 3-6, 2012, Lake Tahoe, Nevada, United States., pages 2960–2968, 2012.
• [SMDH13] Ilya Sutskever, James Martens, George E. Dahl, and Geoffrey E. Hinton. On the importance of initialization and momentum in deep learning. In Proceedings of the 30th International Conference on Machine Learning, ICML 2013, Atlanta, GA, USA, 16-21 June 2013, pages 1139–1147, 2013.
• [SSA13] Kevin Swersky, Jasper Snoek, and Ryan Prescott Adams. Multi-task bayesian optimization. In Advances in Neural Information Processing Systems 26: 27th Annual Conference on Neural Information Processing Systems 2013. Proceedings of a meeting held December 5-8, 2013, Lake Tahoe, Nevada, United States., pages 2004–2012, 2013.
• [SSZA14] Jasper Snoek, Kevin Swersky, Richard S. Zemel, and Ryan P. Adams. Input warping for bayesian optimization of non-stationary functions. In Proceedings of the 31th International Conference on Machine Learning, ICML 2014, Beijing, China, 21-26 June 2014, pages 1674–1682, 2014.
• [Tib96] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society (Series B), 58:267–288, 1996.
• [WF16] Jian Wu and Peter I. Frazier. The parallel knowledge gradient method for batch bayesian optimization. In Advances in Neural Information Processing Systems 29: Annual Conference on Neural Information Processing Systems 2016, December 5-10, 2016, Barcelona, Spain, pages 3126–3134, 2016.
• [WZH13] Ziyu Wang, Masrour Zoghi, Frank Hutter, David Matheson, and Nando de Freitas. Bayesian optimization in high dimensions via random embeddings. In IJCAI 2013, Proceedings of the 23rd International Joint Conference on Artificial Intelligence, Beijing, China, August 3-9, 2013, pages 1778–1784, 2013.
• [ZK16] Sergey Zagoruyko and Nikos Komodakis. Wide residual networks. CoRR, abs/1605.07146, 2016.

## Appendix A Experimental details

### a.1 Options

See Table 2 for the specific hyperparameter options that we use in Section 5. For those variables with options (), we use binary variables under the same name to represent them. For example, we have two variables (01, 02) and their binary representation to denote four kinds of possible initializations: Xavier Glorot [GB10], Kaiming [HZRS15], , or .

### a.2 Importance features

We show the selected important features and their weights during the first stages in Table 3, where each feature is a monomial of variables with degree at most . We do not include the 4th stage because in that stage there are no features with nonzero weights.

Smart choices on important options. Based on Table 3, Harmonica will fix the following variables (sorted according to their importance): Batch Norm (Yes), Activation (ReLU), Initial learning rate (), Optimization method (Adam), Use momentum (Yes), Resblock first activation (Yes), Resblcok third activation (No), Weight decay (No if initial learning rate is comparatively small and Yes otherwise), Batch norm tuning (Yes). Most of these choices match what people are doing in practice.

A metric for the importance of variables. The features that Harmonica finds can serve as a metric for measuring the importance of different variables. For example, Batch Norm turns out to be the most significant variable, and ReLU is second important. By contrast, Dropout, when Batch Norm is presented, does not have significant contributions. This actually matches with the observations in [IS15].

No dummy/irrelevant variables selected. Although there are dummy variables, we never select any of them. Moreover, the irrelevant variables like cudnn, backend, nthreads, which do not affect the test error, were not selected.

### a.3 Generalizing from small networks to big networks

In our experiments, Harmonica first runs on a small network to extract important features and then uses these features to do fine tuning on a big network. Since Harmonica finds significantly better solutions, it is natural to ask whether other algorithms can also exploit this strategy to improve performance.

Unfortunately, it seems that all the other algorithms do not naturally support feature extraction from a small network. For Bayesian Optimization techniques, small networks and large networks have different optimization spaces. Therefore without some modification, Spearmint cannot use information from the small network to update the prior distribution for the large network.

Random-search-based techniques are able to find configurations with low test error on the small network, which might be good candidates for the large network. However, based on our simulation, good configurations of hyperparameters from random search do not generalize from small networks to large networks. This is in contrast to important features in our (Fourier) space, which do seem to generalize.

To test the latter observation using Cifar-10 dataset, we first spent 7 GPU days on layer network to find top configurations among random selected configurations. Then we apply these configurations, as well as locally perturbed configurations (each of them is obtained by switching one random option from one top- configuration), so in total “promising” configurations, to the large layer network. This simulation takes GPU days, but the best test error we obtained is only , even worse than purely random search. Since Hyperband is essentially a fast version of Random Search, it also does not support feature extraction.

Hence, being able to extract important features from small networks seems empirically to be a unique feature of Harmonica.