Title: OPTIMIZING FOR GENERALIZATION IN MACHINE LEARNING WITH CROSS-VALIDATION GRADIENTS Anonymous

Abstract: Cross-validation is the workhorse of modern applied statistics and machine learning, as it provides a principled framework for selecting the model that maximizes generalization performance. In this paper, we show that the cross-validation risk is differentiable with respect to the hyperparameters and training data for many common machine learning algorithms, including logistic regression, elastic-net regression, and support vector machines. Leveraging this property of differentiability, we propose a cross-validation gradient method (CVGM) for hyperparameter optimization. Our method enables efficient optimization in high-dimensional hyperparameter spaces of the cross-validation risk, the best surrogate of the true generalization ability of our learning algorithm.

1. INTRODUCTION:

The ultimate aim of a supervised learning method is generalization, that is, achieving good prediction ability on unseen test data given only a finite set of training data. The generalization capability of learning algorithms should be the primary criterion for model selection, yet an algorithm's generalization capability is a somewhat elusive quantity that is challenging to optimize for. In this paper we introduce a method to optimize directly for the closest available proxy to generalization performance: cross-validation loss.

We begin with a formal description of the overall goal in predictive learning, which also serves as an introduction to notation used throughout the paper. The task of predictive learning involves deriving a prediction function from a finite set of training data. More formally, suppose that (x, y) ∈ X × Y have some joint probability distribution. We have access to a finite dataset of N training examples z i ∈ Z = X × Y drawn i.i.d. from the joint distribution, denoted S = {z 1 , z 2 , . . . , z N }.

We are given (or specify ourselves) a cost function c(ŷ, y) : Y × Y → R + that quantifies the displeasure incurred when ŷ is predicted instead of y. Denoting the function space from input to outputs as F = Y X , we define the loss of a function on a training example z = (x, y) as l(f, z) = c(f (x), y). Then, given a prediction function f ∈ F, we define the population risk as

E z [l(f, z)],

and the target function f * ∈ F as the function that minimizes the population risk. The population risk represents how much loss we incur, on average, on the full joint distribution, and is the quantity we would like as small as possible. In this paper, we consider parametric prediction functions, that is, f is parameterized by a vector θ, denoted f (x; θ)1 . For example, in linear regression, X = R n , Y = R, c(ŷ, y) = (ŷ -y) 2 , and F is the set of all affine functions parameterized as f (x; θ) = θ T x + θ 0 .

We are then tasked with designing a learning algorithm A : Z N → F, which is a function that maps a dataset S to a prediction function. Without substantial knowledge of the actual joint distribution, or assumptions about the target function f * , it is extremely unlikely that A will ever reproduce the exact target function. However, our goal is to minimize the population risk of the learning algorithm

R(A, S) = E z [l(A(S), z)]

(1) which is a random variable that depends on S, our dataset. To make this problem of searching for learning algorithms tractable, we similarly parameterize our learning algorithm A by a vector α ∈ R d , denoted A α . These are known as the "hyperparameters" or "meta-parameters" of the learning algorithm, and can play many important roles: they can perform regularization, enforce sparsity, or even guide feature selection (Hastie et al., 2001) . The quantity we would then like to optimize is the expected population risk, or

EQUATION

It is impossible to exactly calculate (2) with a finite dataset S, as there are two expectations that both involve an unknown probability distribution. What we can do is construct a Monte Carlo estimate of the an algorithm's expected population risk using a technique known as cross-validation. We first partition S into K partitions T j , V j , j = 1, . . . , K (that is,

T j ∩ V j = ∅ and T j ∪ V j = [n]

). Then our cross-validation risk, as a function of α, is

EQUATION

and is readily calculated. We first apply the algorithm to each training set and then average the loss on each corresponding validation set. The first sum in (3) corresponds to the expectation in (2), and the second sum corresponds to the expectation in (1). Setting K = 1 reduces to simple outof-sample validation and an arbitrary K reduces to the common K-fold cross-validation estimate (provided T j form a partition of {1, . . . , N } and |T j | = N -N K ). Thus, this formulation can be viewed as a generalization of cross-validation. (See Kohavi & John (1995) for a longer discussion about this general framework.) In cases where the class of models to be used for learning are known, we have reduced the predictive learning problem to the problem of selecting of a hyperparameter vector α to minimize the cross-validation loss. Even in the simplest cases, however, the objective in (3) is nonconvex in α, and in many cases not even continuous (e.g., the 0 -1 classification loss), which can make optimization of this quantity tricky.

1.1. SUMMARY OF RESULTS:

Our first result is to demonstrate that we can find ∇ α L cv (α) for many common convex machine learning algorithms (e.g., logistic regression, elastic-net regression, support vector machines), provided the cross-validation loss function is differentiable (Section 3). In those algorithms, α often plays the role of regularizer or defines a feature map (in the case of SVM kernels). In the case where α is low-dimensional, (3) can be optimized by exhaustive search without incurring too much cost. However, if we want to design our machine learning algorithms with more expressive regularizations or feature maps, exhaustive search over our hyperparameter space becomes prohibitive.

Our second contribution is to propose the cross-validation gradient method (CVGM), which makes it possible to optimize cross-validation loss over high-dimensional hyperparameter spaces via gradient descent techniques (Section 4). We test the CVGM on an elastic-net regression problem to optimize two hyperparameters and on a more ambitious synthetic classification problem to optimize an entire neural network that serves as a kernel function. In this case, the parameters of the neural network are the hyperparameters (Section 5).

2. RELATED WORK:

There have been many proposed approaches for the problem of hyperparameter optimization, which roughly fall into two camps based on whether or not they use gradients.

2.1. GRADIENT-FREE METHODS:

Exhaustive Search Exhaustive search, also known as grid search, restricts the possible set of α to a (finite) set {α 1 , . . . , α n }, usually by discretizing the parameter search space into a regular grid.

Then one exhaustively computes (3) for each α i and chooses the argmin. The main disadvantage of exhaustive search is that its complexity (to find an approximate minimum) scales exponentially with the dimension d, making it prohibitive for practitioners to successfully apply exhaustive search to d greater than 5 or 6.

Random Search Random search for hyperparameters involves repeatedly specifying a probability distribution over R d , sampling from it, and evaluating (3). Quite unintuitively, random search can be more efficient than exhaustive search, even with a simple probability distribution. This is because, in practice, only a few of the hyperparameter dimensions matter (Bergstra & Bengio, 2012) .

Bayesian Optimization Bayesian regression allows us to predict a distribution over L cv (α), further allowing us to query α that maximize a surrogate function, e.g., probability of improvement or expected improvement (Močkus, 1975; Snoek et al., 2012) . The regression is usually carried out with Gaussian processes (GPs) (Rasmussen, 2004) . However, random search still remains a fierce competitor to the (substantially more complicated) Bayesian optimization approach.

2.2. GRADIENT-BASED METHODS:

Implicit Differentiation Most learning algorithms A α are solving some parameterized optimization problem, that is, optimizing some objective function. Under certain conditions, one can apply the well-known implicit function theorem (Dontchev & Rockafellar, 2014) to the optimality conditions of the objective function, and calculate the gradients of the loss function. Larsen et al. (1998) were the first to propose this, in the context of neural networks, when the objective function includes a regularization term that is linear in the regularization parameters. Bengio (2000) further derived the gradients for a general (unconstrained and differentiable) training criterion along with an efficient way of calculating the gradient for a quadratic training criterion, and applied the algorithm to weight decays for linear regression. These results were then extended to support vector machines (SVMs) (Chapelle et al., 2002; Keerthi et al., 2007) and applied to log-linear models (Foo et al., 2008) and ridge regression (Pedregosa, 2016) . This paper seeks to generalize these methods and provide exact conditions under which A α is actually differentiable. In short, when A α is a convex optimization problem parameterized by α, under certain conditions that are satisfied by many common learning algorithms, we can find exact cross-validation gradients.

Iterative Differentiation In addition to approaches based on implicit differentiation, there are also approaches based on iterative differentiation, i.e., they unroll the optimization procedure in A to calculate gradients. Many large-scale machine learning algorithms perform a variation of gradient descent, and since gradient descent is a sequence of analytic updates to the parameters, A α can be unrolled (or "reverse-mode" differentiated) with respect to α by recursively applying the chain rule to the updates in backwards order. Domke (2012) was the first to propose this, deriving backpropagation rules for the heavy-ball method and LBFGS. Since most large-scale machine learning problems in practice are (approximately) solved using variations of the stochastic subgradient method (also known as SGD), the "learning rate" parameter has a large impact on the convergence and training speed of nonconvex models, e.g., neural networks. Maclaurin et al. (2015) extended the results of Domke to the case of stochastic gradient methods, and as a result, the authors were able to update the learning rate throughout the learning process. The advantage of these methods are that they can be applied to any large-scale machine learning problem that uses a stochastic subgradient method. The main limitations of these methods, however, are that the use of finite precision arithmetic when recursively applying the chain rule can lead to inaccuracies in the gradient calculation and that one can encounter exploding or vanishing gradients from repeated application of the chain rule.

3. EXACT DIFFERENTIABILITY OF LEARNING ALGORITHMS:

Recent work by Barratt (2018) provided necessary and sufficient conditions for a parameterized convex optimization problem to be differentiable. We review the results here, and refer the reader to the paper for more details. The setting is a parameterized convex optimization problem minimize f 0 (x, α) subject to f (x, α) 0, h(x, α) = 0.

(4) where x ∈ R n is the optimization variable, the functions f 0 and f are convex for fixed α and h is affine for fixed α. Let s(α) = (x, λ, ν) T denote the optimal x, λ, ν for a given α in (4), where λ and ν are Lagrange multipliers, i.e., that satisfy the Karush-Kuhn-Tucker (KKT) conditions. Then define the vector-valued function

EQUATION

)

where L is the Lagrangian. The main result of the paper is that, for an optimal z = (x, λ, ν),

EQUATION

)

under the assumption that both f i and g are twice differentiable in x and α, strong duality holds, and

∇ α g(x, λ, ν, α) ∈ R(∇ x g(x, λ, ν, α)).

In other words, we can get the derivative of (x, λ, ν) with respect to α. We will focus on the derivative of x with respect to α in this paper, however, it would be interesting to consider the derivative with respect to the dual variables λ and ν.

Since most parametric machine learning procedures can be expressed as parameterized convex programs that satisfy these conditions2 , we can conclude that, in many cases, A α is in fact differentiable.

In fact, many machine learning procedures can even be expressed as quadratic programs (QPs)quadratic objectives with affine inequality and equality constraints -and satisfy the conditions for differentiability, as shown in Amos & Kolter (2017) (assuming a positive definite quadratic). Further, if l is differentiable with respect to the parameters of f , we can use the chain rule to find the gradient of the cross-validation loss with respect to the hyperparameters. We now present several examples of predictive learning algorithms that are in fact differentiable with respect to their hyperparameters. (Additional examples, including the support vector machine, can be found in the Supplementary Materials.)

Example 3.1 (Logistic regression). In logistic regression, X = R n and Y = {-1, 1}. As is standard in classification, we model Prob(x = 1). This probability is represented by the "sigmoid" function p(x; θ) = 1 1+exp(-θ T x) and we minimize a loss function that is proportional to the likelihood of the dataset under this model plus a regularization term

L(θ, C) = 1 2 θ T θ + C N i=1 log(exp(-y i x T i θ) + 1).

This (convex) optimization problem is unconstrained, so the derivative of the optimal solution with respect to the hyperparameter C is just

∇ C A = -∇ 2 θ L(θ, C) -1 ∇ C [∇ θ L(θ, C)]. Letting π i = 1 1+exp(-yiθ T xi) , the gradient is ∇ θ L(θ, C) = θ + C N i=1 (π i -1)y i x i and the Hessian is ∇ 2 θ L(θ, C) = I + CX T DX where D is a diagonal matrix with D ii = π i (1 -π i )

and the rows of X are x i , and is guaranteed to be positive definite. The righthand side is just Calculate the gradient

∇ C ∇ θ f (θ, C) = N i=1 (π i -1)y i x i .

g ← ∇ α   1 K K j=1 1 |V j | i∈Vj l(A α k (T j ), z i )   5:

Update α k+1 using a gradient method with the gradient g 6:

Project α k+1 onto the constraint set 7: end while

We can also take the derivative with respect to the training examples x i (or y i with a similar derivation) using the fact that

∇ xi ∇ θ f (θ, C) = C(π i -1)y i I.

Example 3.2 (Elastic-net regression). In regression, X = R n , Y = R, and c(ŷ, y) = (ŷ -y) 2 . The function f (x; θ) = θ T x, where the intercept term is omitted for illustration. Let the ith row of the data matrix X be equal to x i and the ith entry of the vector y be equal to y i . Elastic-net regression generalizes ridge and LASSO regression and optimizes the squared penalty with a weighted combination of 1 and 2 regularizers (Zou & Hastie, 2005) , or solves the optimization problem

EQUATION

)

The objective is convex, but not differentiable. To transform this into a differentiable parameterized convex optimization problem, we introduce two variables to represent the positive and negative parts of θ, denoted θ p and θ n . Then, letting v = [θ p θ n ] T , elastic-net can be expressed as the following quadratic program (QP) with 2n variables minimize

EQUATION

Since we can differentiate the solution to positive-definite QPs, we can find the gradient of the optimal solution θ = θ p -θ n with respect to the hyperparameters λ 1 and λ 2 , provided λ 2 > 0.

To the best knowledge of the authors, this is the first derivation of the gradients of the elastic-net solution with respect to elastic-net's hyperparameters.

4. CROSS-VALIDATION GRADIENT METHOD (CVGM):

Building off our findings that the solution to many parametric machine learning procedures are differentiable with respect to their hyperparameters, we can now design an algorithm to minimize (3).

The algorithm is summarized in Algorithm 1. The algorithm essentially performs projected gradient descent on (3), restricting α to a pre-defined constraint set. It runs the learning algorithm on the training part of each cross-validation split (line 3), then calculates the loss on the held-out part of each cross-validation split, and then uses the chain rule to calculate their gradients, which are then averaged (line 4). This averaged gradient is then used to update α in a first-order gradient method (line 5), and then α is projected back onto the constraint set (line 6). Once we run the CVGM to find α * , we then run the learning algorithm on the full dataset to find the final prediction function A α * (S).

There are several advantages to this method. First, it directly optimizes the quantity of interest using a gradient-based method, which can be much faster than exhaustive search. Second, if one is smart with their implementation, computing the gradient in line 4 of the algorithm costs little on top of evaluating the function L cv (α) itself. (See Barratt (2018) and Amos & Kolter (2017) for a discussion of this.) Third, our method plays well with parallel computation. The majority of computation time is spent finding the gradient in line 4 of the algorithm. Since the gradient operation is linear, we can split up K runs of the learning algorithm on the K datasets over K processors or compute nodes and then average the resulting gradients. Also, the algorithm can be run in parallel with different random initializations of α to find multiple hyperparameter settings.

There are several immediate improvements that can be made to the CVGM as stated. One improvement would be to make the sampling of cross-validation splits uniform, that is, each index appears an equal number of times in all of the V j and V j . This ensures that each data point shows up an equal number of times in the cross-validation loss. Another improvement would be to use a more sophisticated optimization method, e.g., accelerated or adaptive methods, but in our experiments we just use a gradient method with constant step size and found that it works quite well.

Our method requires two parameters: the number of partitions K (the batch size), and the fraction of samples to include in the training set partition p. We expect that a value of K between 16 and 128 and p > 1 2 should work well in almost all scenarios. A larger K leads to reduced variance, and a larger p leads to a reduced number of examples held out for validation.

5. NUMERICAL EXPERIMENTS:

We evaluate our method on synthetic regression and classification data, noting that further in-depth comparison on real datasets is needed in future work. One benefit of small synthetic experiments is that the true population risk is readily calculated, and also we can evaluate the performance of the method in the low-data regime. All of the code to run our experiments is freely available online3 .

5.1. SYNTHETIC REGRESSION DATA:

First, we evaluate our method on synthetic regression data. There are N = 30 observations and n = 10 features. However, only 8 of the features have non-zero coefficients. We generate data via the following scikit-learn (Pedregosa et al., 2011) command: X, y, coef = make_regression(N, n, n_informative=8, noise=100., \ tail_strength=0., coef=True) We also generate a test set of 1000 examples with the same command for evaluation. As is standard practice in machine learning, we normalize the features-that is, we normalize each feature to mean 0 and standard deviation 1 across the training set and then this same normalization is applied to the validation/test set before prediction.

We run an elastic-net regression (see Example 3.2) to learn a linear prediction function. For simplicity, the (unpenalized) intercept is learned using standard linear regression and then subtracted from y. For our projection step (line 6 of the algorithm), we require λ 2 ≥ , where = 1 × 10 -7 , and λ 1 ≥ 0. We use K = 128, p = .95, and a gradient descent step size of 2 × 10 -4 . The method is implemented in PyTorch using the qpth library, which is a fast, batched, and differentiable QP library, making the algorithm efficient and scalable (Amos & Kolter, 2017) . The authors note, however, that one could create a much faster implementation by making the solver specialized for elastic-net regression.

We compared the CVGM with exhaustive and random search, noting, however, that the hyperparameter optimization problem is in two dimensions and exhaustive/random search are likely to be quite competitive. We ran CVGM with an initial λ 1 = 1 × 10 -2 and λ 2 = 1 × 10 -4 for 100 steps.

For exhaustive search, we did a grid search over a log scale for

λ 1 ∈ [1 × 10 -4 , 1 × 10 -1 ] and λ 2 ∈ [1 × 10 -4 , 1 × 10 -1 ],

and kept the hyperparameters that achieved the lowest cross-validation loss. For random search, we sampled uniformly at random in a log scale from the same variable ranges as exhaustive search. The resulting final test losses (at iteration 100) of this experiment are in Table 1 and we also included a plot of the (test loss) progress of the algorithms in Figure 3 in the Supplementary Materials. CVGM ultimately achieves a test loss of 1.080, lower than the other two methods, and the test loss (which CVGM has no access to but we do compute during training) is for the most part monotonically decreasing throughout the procedure.

5.2. SYNTHETIC CLASSIFICATION DATA:

Next, we experiment with CVGM's ability to learn kernels from scratch on two dimensional synthetic classification data. We first generate a two dimensional dataset of N examples in polar coordinates from two classes that form rings of different radii and have significant overlap. One class has the distribution r ∼ N (1, .4) and θ ∼ Unif(-π, π), and the other class has the distribution r ∼ N (2, .4) and θ ∼ Unif(-π, π). The data is then transformed into Cartesian coordinates using the transformation (x, y) = (r cos θ, r sin θ). A training dataset of size N = 60 is displayed in the top part of Figure 1 . Clearly, the Bayes decision rule for this dataset is to separate the classes at x 2 = 1.5, and the best a linear classifier can do is 50 % test error. But for the sake of illustration of our method, we seek to learn a (parameterized) kernel φ that transforms x into a space where the data is linearly separable, or at least to a space where we can achieve low misclassification loss by learning a linear classifier with logistic regression. We will use a one-layer neural network kernel φ α : R 2 → R 2 , or

φ α (x) = W 2 σ(W 1 x + b 1 ) + b 2 with parameters α = [W 1 , b 1 , W 2 , b 2 ].

In our experiments, we use σ(x) = max(0, x), where the maximum is taken element-wise, and W 1 ∈ R 64×2 and W 2 ∈ R 2×64 . We will transform the data into the new two dimensional space using the neural network, and then fit a linear classifier there using logistic regression (see Example 3.1). In other words, given φ α , we minimize the following objective

L(θ) = 1 2 θ T θ + C N i=1 log(exp(-y i v T i θ) + 1).

where v i = φ α (x i ). We can then find the Jacobian of the optimal solution θ of this objective ∂θ vi using arguments in Example 3.1, and then using the chain rule to find derivatives of the optimal solution with respect to the neural network's parameters. We use the (differentiable) soft-margin loss for the cross-validation loss function l, thereby allowing us to find the derivatives of the crossvalidation loss with respect to the neural network parameters. A nice interpretation is that we are learning a two-layer neural network that first is fed to φ α and then to the logistic regression layer, but the first part of the neural network is trained using the CVGM, and the second is learned through (standard) logistic regression. The (differentiable) logistic regression layer is implemented as a modular PyTorch Function, and is in the source code. In our experiments, we fix C = 10, K = 256, p = .95, and use a gradient descent step size of 1 × 10 -1 . For the rest of the details of the experiments, we refer the reader to the source code.

The kernel manifold at select iterations of the CVGM is displayed in the bottom part of Figure 1 . After about 10 iterations, the method is able to learn a manifold under which the data is (approximately) linearly separable and achieves a test accuracy of 86.5 %, in comparison to the Bayes-optimal accuracy on that test set of 89.2 %.

In a separate experiment, we compared three separate methods: CVGM method with the neural network kernel as described above, a two-layer neural network with the same architecture as the CVGM method, and logistic regression. Training was done on dataset sizes from 8 to 200 over 25 random seeds. The CVGM model was trained using gradient descent on the binary cross entropy loss with a step size of 1 × 10 -2 for 100 steps, as we found that optimizing to convergence led to severe overfitting-this overfitting is more pronounced when there is less data and is likely a consequence of the low-data regime of this experiment. The mean test accuracies of the various learning algorithms over the random seeds, as well as the Bayes accuracy, are displayed in Figure 2 . CVGM outperforms the two other methods, especially in the low-data regime One important observation is that CVGM is able to stably learn 322 hyperparameters with as few as 8 data points.

6. CONCLUSION:

By showing that we can in fact differentiate the optimal solution to most convex machine learning algorithms, we have made the cross-validation loss, which is commonly viewed as a black-box function, a differentiable objective function. This opens up the possibility of optimizing over large hyperparameter spaces, as demonstrated by our second experiment, where we optimized 322 hyperparameters with as few as 20 training examples.

Practitioners know that of the most important parts of machine learning pipelines is feature engineering, which involves applying some function to raw data before feeding it to a machine learning algorithm. Typically, the optimal features are problem-dependent, requiring experts to spend time constructing and experimenting with hand-crafted functions. However, with the CVGM, practitioners can design differentiable parameterized feature engineering functions for their class of problems, and optimize the feature engineering pipeline directly for generalization capability using gradient descent. Hence, we believe that the CVGM we present is a step towards robust automatic feature learning, a prized goal of machine learning research.

Roughly, this corresponds to the difference in the expected loss between A not having access to z i and the algorithm having access to z i . The main theorem of the paper relates the empirical risk to population risk E[R -R emp ] = ∆. Because our goal is to optimize E[R], we can achieve this by optimizing the quantity ∆ + E[R emp ], or ∆ + E[R emp ] = E l(A(S), z i ) -l(A(S i ), z i ) + l(A(S), z i ) . The last two terms cancel, leaving us with E[l(A(S), z i )],

which we approximate with (3). Hence, our CVGM algorithm is implicitly choosing the hyperparameters to optimize the stability of the learning algorithm. 

D SUPPLEMENTARY FIGURES FOR NUMERICAL EXPERIMENTS:



E XOR EXPERIMENT:

We also experimented with learning a two-dimensional XOR function. The data (N = 100) comes from two classes. One class comes from (x, y) ∼ U[-3, 0.6], U[-0.6, 3] or (x, y) ∼ U[-0.6, 3], U[-3, 0.6] with equal probability. The other class comes from (x, y) ∼ U[-3, 0.6], U[-3, 0.6] or (x, y) ∼ U[-0.6, 3], U[-0.6, 3] with equal probability. The details are similar to our classification experiment, but instead we use a two-layer neural network, with 64 hidden units in the first layer, and 64 hidden units in the second. (This corresponds to 4482 hyperparameters.) The results of this experiment are displayed in Figure 4 . 

:

Nonparametric learning algorithms exist, e.g., k-nearest neighbor, but are challenging to analyze with our method.

Two notable exceptions to this are neural networks and decision trees, which both have nonconvex training criterions and thus one cannot guarantee finding a global minimum.

www.github.com/anonymous/crossval

Figures and Tables:

Algorithm 1 CVGM Require: α 0 : Initial hyperparameter vector Require: K: Number of partitions Require: p: Fraction of samples in training set 1: Sample K partitions T j , V j , j = 1, . . . , K uniformly at random, where |T j | = pN 2: while α k not converged do 3:Run the algorithm A α k (T j ) separately on each training set 4:

Figure 1: Top: Synthetic classification training data in two dimensions, along with the learned decision boundary (black). Bottom: Visualization of kernel applied to 1000 (unseen) test data points in iterations 1, 10, 20, and 100 from left to right, along with the linear classifier in that space (black). Note that the kernel quickly learns a manifold where the data is (approximately) linearly separable. The CVGM achieves 89.8 % test accuracy on the test set with only 60 training examples, and the Bayes (optimal) accuracy is 89.4 %.

Figure 2: Test loss of the CVGM method, a neural network, and vanilla logistic regression for various dataset sizes.

Figure 3: Synthetic regression data. Test loss of random search, CVGM, and exhaustive search, each run for 100 iterations.

Figure 4: XOR classification experiment. Top: Synthetic training data in two dimensions, along with the learned decision boundary (black). Bottom: Visualization of kernel applied to (unseen) test data points in iterations 1, 2, 5, and 50 from left to right, along with the linear classifier in that space (black). The classifier achieves 72 % test accuracy, where the Bayes accuracy is 85.1 %.

Synthetic Regression Results.

Table content: <table><tr><td>Name</td><td colspan="2">Test Loss Cross-Validation Steps</td></tr><tr><td colspan="2">CVGM Exhaustive search 1.207 1.080 Random search 1.084</td><td>100 100 100</td></tr></table>

Peer review (iclr2019):
 The paper considers the problem of automatic tuning of hyperparameters in machine learning models. To address this problem the authors propose to use the so called cross-validation gradients, which optimize a validation objective with respect to the hyperparameters of a model. This approach and the investigated setting falls into a class of optimization problems known as bilevel optimization. The main characteristic of this class of optimization problems is the nested structure, with an outer and inner optimization objectives/problems. The outer problem corresponds to the validation objective and it is defined via an optimal solution to the inner problem which corresponds to a training objective. The paper, however, fails to make a reference to a rather rich literature on bilevel optimization (e.g., see [3-5] and references therein). The approach, presented as Algorithm 1, does not seem different from [1] and [2] where hyperparameter optimization was considered for (kernel) support vector machines and (kernel) ridge regression. The data is initially split into k-folds (not necessarily of identical size) and each fold is used exactly once to define a validation objective whereas the complementary folds act as training data. The validation gradient is obtained by averaging the gradients of the k validation folds. Essentially, the same algorithm with k-fold cross-validation was considered in [1]. Thus, for me there does not seem to be any novelty in this approach and the paper itself. The experiments involve synthetic regression and classification datasets but there are no novel insights that advance what is already known about the hyperparameter optimization (e.g., see [3]). For example, there is no intuition on the geometry of the optimization problem and the optimality of the outer optimization problem which is non-convex (e.g., see [7]), or dependence of the outer solution on the accuracy of the inner solution. References: [1] S. Keerthi, V. Sindhwani, and O. Chapelle. An Efficient Method for Gradient-Based Adaptation of Hyperparameters. NIPS 2007. [2] O. Chapelle, V. Vapnik, O. Bousquet, and S. Mukherjee. Choosing Multiple Parameters for Support Vector Machines. Machine Learning, 2002. [-3] L. Franceschi, P. Frasconi, S. Salzo, R. Grazzi, and M. Pontil. Bilevel Programming for Hyperparameter Optimization and Meta-Learning. ICML 2018. [-4] G. Kunapuli, K.P. Bennet, J. Hu, and J-S. Pang. Bilevel Model Selection for SVMs. American Mathematical Society, 2008. [-5] E.S.H. Neto and A.R. de Pierro. On Perturbed Steepest Descent Methods with Inexact Line Search for Bilevel Convex Optimization. Journal of Mathematical Programming and Operations Research, 2011. [-6] B. Colson, P. Marcotte, and G. Savard. A Trust-Region Method for Nonlinear Bilevel Programming: Algorithms and Computational Experience. Computational Optimization and Applications, 2005. [-7] M. Janzamin, H. Sedghi, and A. Anandkumar. Beating the Perils of Non-Convexity: Guaranteed Training of Neural Networks using Tensor Method. arXiv preprint arXiv:1506.08473v3, 2016.