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, fullyconnected, 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 autotuning, or finding a good setting of these parameters, is now referred to as hyperparameter optimization (HPO), or simply automatic machine learning (autoML). 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
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:

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.

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

is structured.
The third point is very important since clearly HPO is informationtheoretically 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. Multiarmed 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: sampleefficient 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 stateoftheart solvers from Bayesian optimization, Multiarmed bandit techniques, and Random Search. Projecting to even higher numbers of hyperparameters, we perform simulations that show several ordersofmagnitude of speedup versus Bayesian optimization techniques.
1.2 Previous work
The literature on discretedomain HPO can be roughly divided into two: probabilistic approaches and decisiontheoretic 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. Gradientbased 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).
Decisiontheoretic 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 multiarmed 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 lowdegree polynomials) used the celebrated “lowdegree” algorithm of Linial, Mansour, and Nisan [LMN93]
. Their algorithm uses random sampling to estimate each lowdegree Fourier coefficient to high accuracy.
We make use of the approach of Stobbe and Krause [SK12], who showed how to learn lowdegree, 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 lowdegree” (in the sense that most of the mass of the Fourier spectrum resides on monomials of lowdegree). 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 informationtheoretic 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, realvalued 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 ifDefinition 1.
[Rau10] We say a family of functions ( maps to ) is a Random Orthonormal Family with respect to if
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
where
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 lowdegree, 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 equalwhere is assumed to be zeromean, 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 ):
(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 nonzero 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 :
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.
(2) 
Theorem 6 (Noiseless recovery).
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 lowdegree:
Lemma 7 (Noisy recovery).
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 :
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
(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 varianceWith 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
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
squarederror at most via rescaling 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 learningtheoretic 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, lowdegree 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 nonzero 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 Quasipolynomial Time and Polynomial Sample Complexity
Here we observe that our results imply new bounds for decisiontree learning. For example, we obtain the first quasipolynomial 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 squarederror 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 “LowDegree” Algorithm. Prior work for learning decision trees (more generally Boolean functions that are approximated by lowdegree polynomials) used the celebrated “lowdegree” algorithm of Linial, Mansour, and Nisan [LMN93]. Their algorithm uses random sampling to estimate each lowdegree Fourier coefficient to high accuracy. In contrast, the compressedsensing approach of Stobbe and Krause [SK12] takes advantage of the incoherence of the design matrix and gives results that seem unattainable from the “lowdegree” algorithm.
For learning noiseless, Boolean decision trees, the lowdegree 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 realvalued 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 lowdegree 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 nonconvex 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 tiebreaking 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).
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 Harmonica^{2}^{2}2A python implementation of Harmonica can be found at https://github.com/callowbird/Harmonica with Spearmint^{3}^{3}3https://github.com/HIPS/Spearmint.git [SLA12], Hyperband, SH^{4}^{4}4We implemented a parallel version of Hyperband and SH in Lua. and Random Search. Both Spearmint and Hyperband are stateoftheart algorithms, and it is observed that Random Search 2x (Random Search with doubled function evaluation resources) is a very competitive benchmark that beats many algorithms^{5}^{5}5E.g., see [Rec16a, Rec16b]. .
Our first experiment is over training residual network on Cifar10 dataset^{6}^{6}6https://github.com/facebook/fb.resnet.torch. 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 handtuning 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 stagebystage. More specifically, during the feature selection stages, we run Harmonica for tuning an
layer neural network withtraining 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 epochs^{7}^{7}7Other 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 (Harmonica1), 2 stages (Harmonica2) and 3 stages (Harmonica3). 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:
Harmonica1 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]). Harmonica2 uses slightly more time, but is able to find better results, which are comparable with Spearmint with running time.
Improving upon humantuned parameters:
Harmonica3 obtains a better test error () as compared to the best handtuning rate reported in [HZRS16]^{8}^{8}8 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 Cifar10.. Harmonica3 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 Harmonica3. In Section 5.3 we show that by running Harmonica3 for longer time, one can obtain a few other solutions better than hand tuning.
Performance of provable methods:
Harmonica1 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 (Harmonica1, Harmonica2, Harmonica3), with 1 stage and 2 stages using Random Search as the base algorithm (Harmonica1RandomSearch, Harmonica2RandomSearch), and with 2 stages and 3 stages running SH as the base for longer time (Harmonica2Long, Harmonica3Long). 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.
Parameter  Stage 1  Stage 2  Stage 3 

#Samples 
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 hyperparameter 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 hyperparameter optimization. In J. ShaweTaylor, 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. Gradientbased optimization of hyperparameters. Neural Computation, 12(8):1889–1900, 2000.
 [BKW03] Avrim Blum, Adam Kalai, and Hal Wasserman. Noisetolerant 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 TatSeng Chua. Drmad: Distilling reversemode 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 1315, 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, 2126 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 listdecoding size of fouriersparse boolean functions. In David Zuckerman, editor, 30th Conference on Computational Complexity, CCC 2015, June 1719, 2015, Portland, Oregon, USA, volume 33 of LIPIcs, pages 58–71. Schloss Dagstuhl  LeibnizZentrum fuer Informatik, 2015.
 [HR16] Ishay Haviv and Oded Regev. The restricted isometry property of subsampled fourier matrices. In Proceedings of the Twentyseventh Annual ACMSIAM 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 humanlevel performance on imagenet classification.
In2015 IEEE International Conference on Computer Vision, ICCV 2015, Santiago, Chile, December 713, 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 2730, 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, 611 July 2015, pages 448–456, 2015.
 [JT16] Kevin G. Jamieson and Ameet Talwalkar. Nonstochastic best arm identification and hyperparameter optimization. In Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, AISTATS 2016, Cadiz, Spain, May 911, 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 813 2014, Montreal, Quebec, Canada, pages 3122–3130, 2014.
 [LBGR15] Jelena Luketina, Mathias Berglund, Klaus Greff, and Tapani Raiko. Scalable gradientbased 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 BanditBased Approach to Hyperparameter Optimization. ArXiv eprints, 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. Gradientbased 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. http://www.argmin.net/2016/06/23/hyperband/, 2016.
 [Rec16b] Benjamin Recht. The news on autotuning. http://www.argmin.net/2016/06/20/hypertuning/, 2016.
 [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 2123, 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 36, 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, 1621 June 2013, pages 1139–1147, 2013.
 [SSA13] Kevin Swersky, Jasper Snoek, and Ryan Prescott Adams. Multitask 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 58, 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 nonstationary functions. In Proceedings of the 31th International Conference on Machine Learning, ICML 2014, Beijing, China, 2126 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 510, 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 39, 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
Option Name  Description 
01. Weight initialization  Use standard initializations or other initializations? 
02. Weight initialization (Detail 1)  Xavier Glorot [GB10], Kaiming [HZRS15], , or ? 
03. Optimization method  SGD or ADAM? [KB14] 
04. Initial learning rate  or ? 
05. Initial learning rate (Detail 1)  , , , or ? 
06. Initial learning rate (Detail 2)  0.3, 0.1, 0.03, 0.01, 0.003, 0.001, 0.0003, or 0.0001? 
07. Learning rate drop  Do we need to decrease learning rate as we train? Yes or No? 
08. Learning rate first drop time  If drop learning rate, when is the first time to drop by ? Epoch 40 or Epoch 60? 
09. Learning rate second drop time  If drop learning rate, when is the second time to drop by ? Epoch 80 or Epoch 100? 
10. Use momentum [SMDH13]  Yes or No? 
11. Momentum rate  If use momentum, rate is 0.9 or 0.99? 
12. Initial residual link weight  What is the initial residual link weight? All constant 1 or a random number in ? 
13. Tune residual link weight  Do we want to use back propagation to tune the weight of residual links? Yes or No? 
14. Tune time of residual link weight  When do we start to tune residual link weight? At the first epoch or epoch 10? 
15. Resblock first activation  Do we want to add activation layer after the first convolution? Yes or No? 
16. Resblock second activation  Do we want to add activation layer after the second convolution? Yes or No? 
17. Resblock third activation  Do we want to add activation layer after adding the residual link? Yes or No? 
18. Convolution bias  Do we want to have bias term in convolutional layers? Yes or No? 
19. Activation  What kind of activations do we use? ReLU or others? 
20. Activation (Detail 1)  ReLU, ReLU, Sigmoid, or Tanh? 
21. Use dropout [SHK14]  Yes or No? 
22. Dropout rate  If use dropout, rate is high or low? 
23. Dropout rate (Detail 1)  If use dropout, the rate is 0.3, 0.2, 0.1, or 0.05? 
24. Batch norm [IS15]  Do we use batch norm? Yes or No? 
25. Batch norm tuning  If we use batch norm, do we tune the parameters in the batch norm layers? Yes or No? 
26. Resnet shortcut type  What kind of resnet shortcut type do we use? Identity or others? 
27. Resnet shortcut type (Detail 1)  Identity, Identity, Type B or Type C? 
28. Weight decay  Do we use weight decay during the training? Yes or No? 
29. Weight decay parameter  If use weight decay, what is the parameter? or ? 
30. Batch Size  What is the batch size we should use? Big or Small? 
31. Batch Size (Detail 1)  256, 128, 64, or 32? 
32. Optnet  An option specific to the code^{9}^{9}9https://github.com/facebook/fb.resnet.torch. Yes or No? 
33. Share gradInput  An option specific to the code. Yes or No? 
34. Backend  What kind of backend shall we use? cudnn or cunn? 
35. cudnn running state  If use cudnn, shall we use fastest of other states? 
36. cudnn running state (Detail 1)  Fastest, Fastest, default, deterministic 
37. nthreads  How many threads shall we use? Many or few? 
38. nthreads (Detail 1)  8, 4, 2, or 1? 
3960. Dummy variables  Just dummy variables, no effect at all. 
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.
Stage  Feature Name  Weights 

11  24. Batch norm  8.05 
12  19. Activation  3.47 
13  04. Initial learning rate * 05. Initial learning rate (Detail 1)  3.12 
14  19. Activation * 24. Batch norm  2.55 
15  04. Initial learning rate  2.34 
16  28. Weight decay  1.90 
17  24. Batch norm * 28. Weight decay  1.79 
18  34. Optnet * 35. Share gradInput * 52. Dummy ^{10}^{10}10This is an interesting feature. In the code repository that we use, optnet, shared gradInput are two special options of the code and cannot be set true at the same time, otherwise the training becomes unpredictable.  1.54 
21  03. Optimization method  4.22 
22  03. Optimization method * 10. Use momentum  3.02 
23  15. Resblock first activation  2.80 
24  10. Use momentum  2.19 
25  15. Resblock first activation * 17. Resblock third activation  1.68 
26  01. Good initialization  1.26 
27  01. Good initialization * 10. Use momentum  1.12 
28  01. Good initialization * 03. Optimization method  0.67 
31  29. Weight decay parameter  0.49 
32  28. Weight decay  0.26 
33  06. Initial learning rate (Detail 3) * 28. Weight decay  0.23 
34  25. Batch norm tuning  0.21 
35  28. Weight decay * 29. Weight decay parameter  0.20 
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.
Randomsearchbased 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 Cifar10 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.