Year: 2018

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.

heading: 1 INTRODUCTION
text: 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 zi ∈ Z = X × Y drawn i.i.d. from the joint distribution, denoted
S = {z1, z2, . . . , zN}.
We are given (or specify ourselves) a cost function c(ŷ, y) : Y × Y → R+ that quantifies the displeasure incurred when ŷ is predicted instead of y. Denoting the function space from input to outputs as F = YX , 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 = Rn, Y = R, c(ŷ, y) = (ŷ − y)2, and F is the set of all affine functions parameterized as f(x; θ) = θTx+ θ0.
We are then tasked with designing a learning algorithmA : ZN → 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
1Nonparametric learning algorithms exist, e.g., k-nearest neighbor, but are challenging to analyze with our method.
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 α ∈ Rd, 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
L(α) = E S
[R(Aα, S)] . (2)
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 Tj , Vj , j = 1, . . . ,K (that is, Tj ∩ Vj = ∅ and Tj ∪ Vj = [n]). Then our cross-validation risk, as a function of α, is
Lcv(α) = 1
K K∑ j=1 1 |Vj | ∑ i∈Vj l(Aα(Tj), zi) (3)
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 Tj form a partition of {1, . . . , N} and |Tj | = N − NK ). 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.

heading: 1.1 SUMMARY OF RESULTS
text: Our first result is to demonstrate that we can find ∇αLcv(α) 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).

heading: 2 RELATED WORK
text: 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.

heading: 2.1 GRADIENT-FREE METHODS
text: 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 Rd, 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 Lcv(α), further allowing us to query α that maximize a surrogate function, e.g., probability of improvement or expected improvement (Moč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.

heading: 2.2 GRADIENT-BASED METHODS
text: 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.

heading: 3 EXACT DIFFERENTIABILITY OF LEARNING ALGORITHMS
text: 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 f0(x, α) subject to f(x, α) 0,
h(x, α) = 0.
(4)
where x ∈ Rn is the optimization variable, the functions f0 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
g(x̃, λ̃, ν̃, α) = ∇xL(x̃, λ̃, ν̃, α)diag(λ̃)f(x̃, α) h(x̃, α)  . (5) where L is the Lagrangian. The main result of the paper is that, for an optimal z = (x̃, λ̃, ν̃),
∇αs(α) = −∇zg(x̃, λ̃, ν̃, α)−1∇αg(x̃, λ̃, ν̃, α). (6) under the assumption that both fi and g are twice differentiable in x and α, strong duality holds, and ∇αg(x̃, λ̃, ν̃, α) ∈ R(∇xg(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 = Rn andY = {−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(−yixTi θ) + 1).
This (convex) optimization problem is unconstrained, so the derivative of the optimal solution with respect to the hyperparameter C is just
∇CA = −∇2θL(θ, C)−1∇C [∇θL(θ, C)]. Letting πi = 11+exp(−yiθT xi) , the gradient is
∇θL(θ, C) = θ + C N∑ i=1 (πi − 1)yixi
and the Hessian is ∇2θL(θ, C) = I + CXTDX where D is a diagonal matrix with Dii = πi(1− πi) and the rows of X are xi, and is guaranteed to be positive definite. The righthand side is just
∇C∇θf(θ, C) = N∑ i=1 (πi − 1)yixi.
2Two 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.
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 Tj ,Vj , j = 1, . . . ,K uniformly at random, where |Tj | = bpNc 2: while αk not converged do 3: Run the algorithm Aαk(Tj) separately on each training set 4: Calculate the gradient
g ← ∇α  1 K K∑ j=1 1 |Vj | ∑ i∈Vj l(Aαk(Tj), zi)  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 xi (or yi with a similar derivation) using the fact that
∇xi∇θf(θ, C) = C(πi − 1)yiI.
Example 3.2 (Elastic-net regression). In regression, X = Rn, Y = R, and c(ŷ, y) = (ŷ − y)2. The function f(x; θ) = θTx, where the intercept term is omitted for illustration. Let the ith row of the data matrix X be equal to xi and the ith entry of the vector y be equal to yi. 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
minimize 1
2N ‖Xθ − y‖22 + λ1‖θ‖1 +
1 2 λ2‖θ‖22. (7)
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 1 2 vT [ 1 NX TX + λ2I − 1NX TX
− 1NX TX 1NX TX + λ2I
] v + [ − 1NX
T y + λ1I 1 NX T y + λ1I
]T v
subject to v 0. (8)
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.

heading: 4 CROSS-VALIDATION GRADIENT METHOD (CVGM)
text: 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 Lcv(α) 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 Vj and Vj . 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 > 12 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.

heading: 5 NUMERICAL EXPERIMENTS
text: 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.

heading: 5.1 SYNTHETIC REGRESSION DATA
text: 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:

heading: X, y, coef = make_regression(N, n, n_informative=8, noise=100., \
text: 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
3www.github.com/anonymous/crossval
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.

heading: 5.2 SYNTHETIC CLASSIFICATION DATA
text: 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 φα : R2 → R2, or φα(x) = W2σ(W1x+ b1) + b2
with parameters α = [W1, b1,W2, b2]. In our experiments, we use σ(x) = max(0, x), where the maximum is taken element-wise, and W1 ∈ R64×2 and W2 ∈ R2×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(−yivTi θ) + 1).
where vi = φα(xi). 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.

heading: 6 CONCLUSION
text: 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.

heading: A SUPPORT VECTOR MACHINE EXAMPLE
text: Support vector machines perform classification, where X = Rn and Y = {0, 1}. The function class is again linear, or f(x; θ) = θTx. The loss function used in SVMs is the hinge loss, or c(ŷ, y) = max(0, 1− yŷ). In the `1 and `2-regularized SVM, we optimize
L(θ) = 1
N N∑ i=1 max(0, 1− yiθTxi) + λ1‖θ‖1 + 1 2 λ2‖θ‖22.
Introducing the vectors θp, θn (where θ = θp − θn), variables ti, i = 1, . . . , N , we can rewrite the problem as
minimize 1
N N∑ i=1 ti + λ11 T (θp + θn) + 1 2 λ2‖θp − θn‖22
subject to θp 0, θn 0 ti ≥ 0, i = 1, . . . , N 1− yi(θp − θn)Txi ≤ ti, i = 1, . . . , N ,
(9)
which is a QP with considerable structure. Thus the solution θ is differentiable with respect to λ1 and λ2, again provided λ2 > 0. A similar method, i.e., replacing xi with φ(xi) can be used to show that kernel-based SVMs are also differentiable with respect to the kernel parameters.sectionB: Learning Loss Functions

heading: B LEARNING LOSS FUNCTIONS
text: Much of classification can be viewed as optimizing a convex surrogate of the 0 − 1 loss function Bartlett et al. (2006). Thus, a reasonable loss function is a convex combination of such convex surrogates. Four common loss functions are:
• hinge: lh(ŷ, y) = max(0, 1− yŷ). • exponential: le(ŷ, y) = exp(−yŷ). • truncated quadratic: lt(ŷ, y) = max{1− yŷ, 0}2. • logistic: ll(ŷ, y) = ln(1 + exp(−2yŷ)).
Except for the hinge loss, all of these loss functions are differentiable. Thus we can define the optimization problem
minimize 1
N N∑ i=1 α1ti + α2le(θ Txi, yi) + α3lt(θ Txi, yi) + α4ll(θ Txi, yi) +R(θ, α)
subject to 1− yiθTxi ≤ ti, ti ≥ 0, i = 1, . . . , N
(10)
that is convex and differentiable in α. We can then run a projected gradient method over α in the probability simplex.

heading: C CONNECTIONS TO STABILITY
text: Bousquet and Elisseeff Bousquet & Elisseeff (2002) introduced several mathematically precise notions of the “stability” of a learning algorithm. They consider a modified dataset, constructed by replacing one element: Si = {z1, . . . , zi−1, z′i, zi+1, . . . , zN .} The stability of a learning algorithm A is then defined as
∆ = E S,z′i
[ l(A(S), z′i)− l(A(Si), z′i) ] .
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−Remp] = ∆. Because our goal is to optimize E[R], we can achieve this by optimizing the quantity ∆ + E[Remp], or
∆ + E[Remp] = E [ l(A(S), z′i)− l(A(Si), 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.

heading: D SUPPLEMENTARY FIGURES FOR NUMERICAL EXPERIMENTS
text: 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.

Review summaries:
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).