Quentin Bertrand
email: quentin.bertrand@inria.fr

Quentin Klopfenstein
email: quentin.klopfenstein@u-bourgogne.fr

Mathurin Massias
email: mathurin.massias@gmail.com

Mathieu Blondel
email: mblondel@google.com

Samuel Vaiter
email: samuel.vaiter@math.cnrs.fr

Alexandre Gramfort
email: alexandre.gramfort@inria.fr

Joseph Salmon
email: joseph.salmon@umontpellier.fr

Massimiliano Pontil

Implicit Differentiation for Fast Hyperparameter Selection in Non-Smooth Convex Learning

Keywords: Convex optimization, hyperparameter optimization, hyperparameter selection, bilevel optimization, Lasso, generalized linear models

Finding the optimal hyperparameters of a model can be cast as a bilevel optimization problem, typically solved using zero-order techniques. In this work we study first-order methods when the inner optimization problem is convex but non-smooth. We show that the forward-mode differentiation of proximal gradient descent and proximal coordinate descent yield sequences of Jacobians converging toward the exact Jacobian. Using implicit differentiation, we show it is possible to leverage the non-smoothness of the inner problem to speed up the computation. Finally, we provide a bound on the error made on the hypergradient when the inner optimization problem is solved approximately. Results on regression and classification problems reveal computational benefits for hyperparameter optimization, especially when multiple hyperparameters are required.

Introduction

Almost all models in machine learning require at least one hyperparameter, the tuning of which drastically affects accuracy. This is the case for many popular estimators, where e λ 1 |β j | + 1 2 e λ 2 β 2 j 1 n X y ∞ sparse log. reg.

1 n n i=1 ln(1 + e -y i X i: β ) e λ |β j | 1 2n X y ∞ dual SVM 1 2 (y X) β 2 -p j=1 β j ι [0,e λ ] (β j ) -
the regularization hyperparameter controls the trade-off between a data fidelity term and a regularization term. Such estimators, including Ridge regression [START_REF] Hoerl | Ridge regression: Biased estimation for nonorthogonal problems[END_REF], Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], sparse logistic regression [START_REF] Koh | An interior-point method for large-scale l1-regularized logistic regression[END_REF], support-vector machine/SVM [START_REF] Boser | A training algorithm for optimal margin classifiers[END_REF][START_REF] Platt | Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods[END_REF] are often cast as an optimization problem (Table 1)

β(λ) ∈ arg min β∈R p Φ(β, λ) f (β) + p j=1 g j (β j , λ) g(β,λ) , (1)
with smooth f : R p → R (i.e., with Lipschitz gradient), proper closed convex (possibly nonsmooth) functions g j (•, λ), and a regularization hyperparameter λ ∈ R r . In the examples of Table 1, the computation of f involves a design matrix X ∈ R n×p ; and the cost of computing ∇f (β) is O(np). In the SVM example, since we consider the dual problem, we chose to reverse the roles of n and p to enforce β ∈ R p . We often drop the λ dependency and write β instead of β(λ) when it is clear from context. For a fixed λ, the issue of solving efficiently Problem (1) has been largely explored. If the functions g j are smooth, one can use solvers such as L-BFGS [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF], SVRG [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF][START_REF] Zhang | Linear convergence with condition number independent access of full gradients[END_REF], or SAGA [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF]. When the functions g j are non-smooth, Problem (1) can be tackled efficiently with stochastic algorithms [START_REF] Pedregosa | Breaking the nonsmooth barrier: A scalable parallel method for composite optimization[END_REF] or using working set methods [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF][START_REF] Tibshirani | Strong rules for discarding predictors in lasso-type problems[END_REF] combined with coordinate descent [START_REF] Tseng | Block-coordinate gradient descent method for linearly constrained nonsmooth separable optimization[END_REF], see overview by [START_REF] Massias | Dual extrapolation for sparse generalized linear models[END_REF]. The question of model selection, i.e., how to select the hyperparameter λ ∈ R r (potentially multidimensional), is more open, especially when the dimension r of the regularization hyperparameter λ is large.

For the Lasso, a broad literature has been devoted to parameter tuning. Under strong hypothesis on the design matrix X, it is possible to derive guidelines for the setting of the regularization parameter λ [START_REF] Lounici | Sup-norm convergence rate and sign concentration property of Lasso and Dantzig estimators[END_REF][START_REF] Bickel | Simultaneous analysis of Lasso and Dantzig selector[END_REF][START_REF] Belloni | Square-root Lasso: pivotal recovery of sparse signals via conic programming[END_REF]. Unfortunately, these guidelines rely on quantities which are typically unknown in practice, and Lasso users still have to resort to other techniques to select the hyperparameter λ.

A popular approach for hyperparameter selection is hyperparameter optimization [START_REF] Kohavi | Automatic parameter selection by minimizing estimated error[END_REF][START_REF] Hutter | Beyond manual tuning of hyperparameters[END_REF][START_REF] Feurer | Hyperparameter optimization[END_REF]: one selects the hyperpa- rameter λ such that the regression coefficients β(λ) minimize a given criterion C : R p → R.

Here C should ensure good generalization, or avoid overcomplex models. Common examples (see Table 2) include the hold-out loss [START_REF] Devroye | Distribution-free performance bounds for potential function rules[END_REF], the cross-validation loss (CV, [START_REF] Stone | Estimating WAIS IQ from Shipley Scale scores: Another cross-validation[END_REF], see Arlot and Celisse 2010 for a survey), the AIC [START_REF] Akaike | A new look at the statistical model identification[END_REF], BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF] or SURE [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] criteria. Formally, the hyperparameter optimization problem is a bilevel optimization problem [START_REF] Colson | An overview of bilevel optimization[END_REF] arg min (2)

Popular approaches to solve (the generally non-convex) Problem (2) include zero-order optimization (gradient-free) techniques such as grid-search, random-search [START_REF] Rastrigin | The convergence of the random search method in the extremal control of a many parameter system[END_REF][START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF][START_REF] Bergstra | Hyperopt: A python library for optimizing the hyperparameters of machine learning algorithms[END_REF] or Sequential Model-Based Global Optimization (SMBO), often referred to as Bayesian optimization [START_REF] Mockus | The bayesian approach to local optimization[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF][START_REF] Brochu | A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[END_REF][START_REF] Snoek | Practical bayesian optimization of machine learning algorithms[END_REF]. Grid-search is a naive discretization of Problem (2). It consists in evaluating the outer function L on a grid of hyperparameters, solving one inner optimization Problem (1) for each λ in the grid (see Figure 1). For each inner problem solution β(λ) , the criterion C( β(λ) ) is evaluated, and the model achieving the lowest value is selected. Random-search has a similar flavor, but one randomly selects where the criterion must be evaluated. Finally, SMBO models the objective function L via a function amenable to uncertainty estimates on its predictions such as a Gaussian process. Hyperparameter values are chosen iteratively to maximize a function such as the expected improvement as described, e.g., by [START_REF] Bergstra | Algorithms for hyper-parameter optimization[END_REF]. However, these zero-order methods share a common drawback: they scale exponentially with the dimension of the search space (Nesterov, 2004, Sec. 1.1.2).

When the hyperparameter space is continuous and the regularization path λ → β(λ) is well-defined and almost everywhere differentiable, first-order optimization methods are well suited to solve the bilevel optimization Problem (2). Using the chain rule, the gradient of 1. For a linear model y = Xβ + ε, with ε ∼ N (0, σ 2 ), the degree of freedom (dof, Efron 1986) is defined as dof(β) = n i=1 cov(yi, (Xβ)i)/σ 2 . 2. The smoothed Hinge loss is given by (

x) = 1 2 -x if x ≤ 0, 1 2 (1 -x) 2 if 0 ≤ x ≤
C(β (λ)
)

Figure 1: 5-fold cross-validation error C(β (λ) ): (top) Lasso CV error with respect to λ for multiple hyperparameter optimization methods on the real-sim data set, and (bottom) elastic net CV error with respect to λ 1 and λ 2 on the rcv1 data set. Crosses represent the 10 (top) or 25 (bottom) first error evaluations for each method.

L with respect to λ, also referred to as the hypergradient, evaluates to

∇ λ L(λ) = Ĵ (λ) ∇C( β(λ) ) , (3)
with Ĵ(λ) ∈ R p×r the Jacobian of the function λ → β(λ) ,

Ĵ(λ)      ∂ β(λ) 1 ∂λ 1 . . . ∂ β(λ) 1 ∂λr . . . . . . . . . ∂ β(λ) p ∂λ 1 . . . ∂ β(λ) p ∂λr      . ( 4
)
An important challenge of applying first-order methods to solve Problem (2) is evaluating the hypergradient in Equation (3). There are three main algorithms to compute the hypergradient ∇ λ L(λ): implicit differentiation [START_REF] Larsen | Design and regularization of neural networks: the optimal use of a validation set[END_REF][START_REF] Bengio | Gradient-based optimization of hyperparameters[END_REF] and automatic differentiation using the reverse-mode [START_REF] Linnainmaa | The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors[END_REF][START_REF] Lecun | Efficient backprop[END_REF] or the forward-mode [START_REF] Wengert | A simple automatic derivative evaluation program[END_REF][START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF][START_REF] Franceschi | Forward and reverse gradient-based hyperparameter optimization[END_REF]. As illustrated in Figure 1, once the hypergradient in Equation (3) has been computed, one can solve Problem (2) with first-order schemes, e.g., gradient descent.

Contributions. We are interested in tackling the bilevel optimization Problem (2), with a non-smooth inner optimization Problem (1). More precisely,

• We show that classical algorithms used to compute hypergradients for smooth inner problem have theoretically grounded non-smooth counterparts. We provide in Theorem 9 an implicit differentiation formula for non-smooth optimization problems. We obtain in Theorem 13, for the first time in the non-smooth case, error bounds with respect to the hypergradient when the inner problem and the linear system involved are only solved approximately. We obtain in Theorem 12 convergence rates on the hypergradient for iterative differentiation of non-smooth optimization problems.

• Based on the former contributions we propose an algorithm to tackle Problem (2). We develop an efficient implicit differentiation algorithm to compute the hypergradient in Equation ( 3), leveraging the sparsity of the Jacobian and enabling the use of stateof-the-art solvers (Algorithm 5). We combine in Algorithm 6 this fast hypergradient computation with a gradient descent scheme to solve Problem (2).

• We provide extensive experiments on diverse data sets and estimators (Section 4). We first show that implicit differentiation significantly outperforms other hypergradient methods (Section 4.1). Then, leveraging sparsity, we illustrate computational benefits of first-order optimization with respect to zero-order techniques for solving Problem (2) on Lasso, elastic net and multiclass logistic regression (Section 4.2).

• We release our implementation as a high-quality, documented and tested Python package: https://github.com/qb3/sparse-ho.

General notation. We write • the Euclidean norm on vectors. For a set S, we denote by S c its complement. We denote [p] = {1, . . . , p}. We denote by (e j ) p j=1 the vectors of the canonical basis of R p . We denote the coordinate-wise multiplication of two vectors u and v by u v, and by u M the row-wise multiplication between a vector and a matrix. The i-th line of the matrix M is M i: and its j-th column is M :j . The spectral radius of a matrix M ∈ R n×n is denoted ρ(M ) = max i |s i | where s 1 , . . . , s n are the eigenvalues of M . For a matrix M , we write that M 0 if M is positive definite. The regularization parameter, possibly multivariate, is denoted by λ = (λ 1 , . . . , λ r ) ∈ R r . Recall that for a locally integrable function f : x ∈ Ω → R, where Ω is an open subset of R n , its weak partial derivative [START_REF] Evans | Measure theory and fine properties of functions[END_REF] with respect to x i in Ω is the locally integrable function g i on Ω such that

Ω g i (x)φ(x)dx = - Ω f (x) ∂φ(x) ∂x i dx , (5)
holds for all functions φ that are continuously differentiable and of compact support. For vector valued function f , we denote by J (x) its weak Jacobian i.e., the matrix composed of weak partial derivatives. An important use of this notation in the rest of the paper is

Ĵ(λ) (∇ λ β(λ) 1 , . . . , ∇ λ β(λ) p ) ∈ R
p×r the weak Jacobian of β(λ) with respect to λ. Convex analysis. For a convex function h : R p → R, the proximal operator of h is defined, for Table 3: Partial derivatives of proximal operators used.

g j (β j , λ) prox g j (•,λ) (z j ) ∂ z prox g j (•,λ) (z j ) ∂ λ prox g j (•,λ) (z j ) e λ β 2 j /2 z j /(1 + e λ ) 1/(1 + e λ ) -z j e λ /(1 + e λ ) 2 e λ |β j | ST(z j , e λ ) | sign(ST(z j , e λ ))
| -e λ sign(ST(z j , e λ ))

e λ1 |β j | + 1 2 e λ2 β 2 j ST(z j ,e λ 1 ) 1+e λ 2 | sign(ST(z j ,e λ 1 ))| 1+e λ 2 -e λ 1 sign(ST(z j ,e λ 1 )) 1+e λ 2 , -ST(z j ,e λ 1 )e λ 2 (1+e λ 2 ) 2 ι [0,e λ ] (β j ) max(0, min(z j , e λ )) 1 ]0,e λ [ (z j ) e λ 1 z j >e λ any x ∈ R p , as: prox h (x) = arg min y∈R p 1 2 x -y 2 + h(y). The subdifferential of h at x is denoted ∂h(x) = u ∈ R p : ∀z ∈ R p , h(z) ≥ h(x) + u (z -x) . A function is said to be smooth if it has Lipschitz gradients. Let f be a L-smooth function. Lipschitz constants of the functions ∇ j f are denoted by L j ; hence for all x ∈ R p , h ∈ R |∇ j f (x + he j ) -∇ j f (x)| ≤ L j |h| .
For a function f , its gradient restricted to the indices in a set S is denoted ∇ S f . For a set Ξ ⊂ R p , its relative interior is noted ri(Ξ), and its indicator function is defined for any

x ∈ R p by ι Ξ (x) = 0 if x ∈ Ξ and +∞ otherwise. A function h : R → R ∪ {+∞} is said to be proper if dom(h) = {x ∈ R : h(x) < +∞} = ∅), and closed if for any α ∈ R, the sublevel set {x ∈ dom(h) : h(x) ≤ α} is a closed set.
For a function ψ : R p × R r → R p , we denote ∂ z ψ the weak Jacobian with respect to the first variable and ∂ λ ψ the weak Jacobian with respect to the second variable. The proximal operator of g(•, λ) can be seen as such a function ψ of β and λ (see Table 1 for examples)

R p × R r → R p (z, λ) → prox g(•,λ) (z) = ψ(z, λ) .
In this case we denote ∂ z prox g(•,λ) ∂ z ψ and ∂ λ prox g(•,λ) ∂ λ ψ. Since we consider only separable penalties g(•, λ), ∂ z prox g(•,λ) is a diagonal matrix, so to make notation lighter, we write ∂ z prox g(•,λ) for its diagonal. We thus have

∂ z prox g(•,λ) = (∂ z prox g j (•,λ) ) j∈[p] ∈ R p (by separability of g) ∂ λ prox g(•,λ) ∈ R p×r .
Explicit partial derivatives formulas for usual proximal operators can be found in Table 3.

Related Work

The main challenge to evaluate the hypergradient ∇ λ L(λ) is the computation of the Jacobian J (λ) . We first focus on the case where Φ(•, λ) is convex and smooth for any λ.

Implicit differentiation. We recall how the implicit differentiation3 formula of the gradient ∇ λ L(λ) is obtained for smooth inner optimization problems. We will provide a generalization to non-smooth optimization problems in Section 3.2.

Theorem 1 [START_REF] Bengio | Gradient-based optimization of hyperparameters[END_REF]. Let β(λ) ∈ arg min β∈R p Φ(β, λ) be a solution of Problem (1). Assume that for all λ > 0, Φ(•, λ) is a convex smooth function, ∇ 2 β Φ( β(λ) , λ) 0, and that for all β ∈ R p , Φ(β, •) is differentiable over ]0, +∞[. Then the hypergradient ∇ λ L(λ) reads

∇ λ L(λ) ∈R r = -∇ 2 β,λ Φ( β(λ) , λ) ∈R r×p ∇ 2 β Φ( β(λ) , λ) ∈R p×p -1 ∇C( β(λ) ) ∈R p . ( 6
)
Proof For a smooth convex function β → Φ(β, λ) the first-order condition writes:

∇ β Φ( β(λ) , λ) = 0 , (7)
for any β(λ) solution of the inner problem. Moreover, if λ → ∇ β Φ( β(λ) , λ) is differentiable, differentiating Equation ( 7) with respect to λ leads to

∇ 2 β,λ Φ( β(λ) , λ) + Ĵ (λ) ∇ 2 β Φ( β(λ) , λ) = 0 . (8)
The Jacobian Ĵ (λ) is computed by solving the following linear system

Ĵ (λ) = -∇ 2 β,λ Φ( β(λ) , λ) ∈R r×p ∇ 2 β Φ( β(λ) , λ) ∈R p×p -1 . (9)
Plugging Equation (9) into Equation (3) yields the desired result.

The computation of the gradient via implicit differentiation (Equation ( 6)) involves the resolution of a p×p linear system (Bengio, 2000, Sec. 4). This potentially large linear system can be solved using different algorithms such as conjugate gradient [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF]Stiefel 1952, as in Pedregosa 2016) or fixed point methods [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF][START_REF] Tseng | Block-coordinate gradient descent method for linearly constrained nonsmooth separable optimization[END_REF]Yun 2009, as in Grazzi et al. 2020). Implicit differentiation has been used for model selection of multiple estimators with smooth regularization term: kernel-based models [START_REF] Chapelle | Choosing multiple parameters for support vector machines[END_REF][START_REF] Seeger | Cross-validation optimization for large scale structured classification kernel methods[END_REF], weighted Ridge estimator [START_REF] Foo | Efficient multiple hyperparameter learning for log-linear models[END_REF], neural networks [START_REF] Lorraine | Optimizing millions of hyperparameters by implicit differentiation[END_REF] or meta-learning [START_REF] Rajeswaran | Meta-learning with implicit gradients[END_REF]. In addition to hyperparameter selection, it has been applied successfully in natural language processing [START_REF] Bai | Deep equilibrium models[END_REF] and computer vision [START_REF] Bai | Multiscale deep equilibrium models[END_REF].

Problem (1) is typically solved using iterative solvers. In practice, the number of iterations is limited to reduce computation time, and also since very precise solutions are generally not necessary for machine learning tasks. Thus, Equation ( 7) is not exactly satisfied at machine precision, and consequently the linear system to solve Equation (6) does not lead to the exact gradient ∇ λ L(λ), see [START_REF] Ablin | Super-efficiency of automatic differentiation for functions defined as a minimum[END_REF] for quantitative convergence results. However, [START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF] showed that one can resort to approximate gradients when the inner problem is smooth, justifying that implicit differentiation can be applied using an approximation of β. Interestingly, this approximation scheme was shown to yield significant practical speedups when solving Problem (2), while preserving theoretical properties of convergence toward the optimum. Practitioners now have access to powerful software to use implicit differentiation with smooth inner optimization problems [START_REF] Blondel | Efficient and modular implicit differentiation[END_REF].

Iterative differentiation. Iterative differentiation computes the gradient ∇ λ L(λ) by differentiating through the iterates of the algorithm used to solve Problem (1). Iterative differentiation can be applied using the forward-mode [START_REF] Wengert | A simple automatic derivative evaluation program[END_REF] or the reverse-mode [START_REF] Linnainmaa | The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors[END_REF]. Both rely on the chain rule, the gradient being decomposed as a large product of matrices, computed either in a forward or backward way. Note that forward and reverse modes are algorithm-dependent: in this section we illustrate iterative differentiation for proximal gradient descent (PGD, [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF][START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], using the forward-mode (Algorithm 1), and the reverse-mode (Algorithm 2).

The most popular method in automatic differentiation is the reverse-mode, a cornerstone of deep learning (Goodfellow et al., 2016, Chap. 8). Iterative differentiation for hyperparameter optimization can be traced back to [START_REF] Domke | Generic methods for optimization-based modeling[END_REF], who derived (for smooth loss functions) a reverse-mode with gradient descent, heavy ball and L-BFGS algorithms. It first computes the solution of the optimization Problem (1) using an iterative solver, but requires storing the iterates along the computation for a backward evaluation of the hypergradient (Algorithm 2). [START_REF] Maclaurin | Gradient-based hyperparameter optimization through reversible learning[END_REF] used the reverse-mode on stochastic gradient descent to select thousands of hyperparameters. Alternatively, the forward-mode computes jointly the solution along with the gradient ∇ λ L(λ). The forward-mode has been applied to hyperparameter optimization with smooth inner problems by [START_REF] Franceschi | Forward and reverse gradient-based hyperparameter optimization[END_REF]. [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF] paved the way for applying it to non-smooth optimization problems. The forward-mode is memory efficient (no iterates storage) but more computationally expensive when the number of hyperparameters (r) is large; see [START_REF] Baydin | Automatic differentiation in machine learning: a survey[END_REF] for a survey.

Resolution of the bilevel Problem (2). From a theoretical point of view, solving Problem (2) using gradient-based methods is also challenging, and results in the literature are quite scarce. [START_REF] Kunisch | A bilevel optimization approach for parameter learning in variational models[END_REF] studied the convergence of a semi-Newton algorithm where both the outer and inner problems are smooth. [START_REF] Franceschi | Bilevel programming for hyperparameter optimization and meta-learning[END_REF] gave similar results with weaker assumptions to unify hyperparameter optimization and meta-learning with a bilevel point of view. They required the inner problem to have a unique solution for all λ > 0 but do not have second-order assumptions on Φ. Recent results [START_REF] Ghadimi | Approximation methods for bilevel programming[END_REF][START_REF] Ji | Provably faster algorithms for bilevel optimization and applications to meta-learning[END_REF] have provided quantitative convergence toward a global solution of Problem (2), under the assumption that the inner problem is strongly convex and one has the exact knowledge of the hypergradient Lipschitz constant.

Algorithm 1 Forward-mode PGD input : λ ∈ R r , γ > 0, n iter ∈ N, β (0) ∈ R p , J (0) ∈ R p×r // jointly compute coef. & Jacobian for k = 1, . . . , n iter do // update the regression coefficients z (k) = β (k-1) -γ∇f (β (k-1) ) // GD step dz (k) = J (k-1) -γ∇ 2 f (β (k-1) )J (k-1) β (k) = prox γg(•,λ) (z (k) )
// prox. step // update the Jacobian

J (k) = ∂ z prox γg(•,λ) (z (k) ) dz (k) J (k) += ∂ λ prox γg(•,λ) (z (k) ) // O(pr) v = ∇C(β n iter ) return β n iter , J n iter v Algorithm 2 Reverse-mode PGD input : λ ∈ R r , γ > 0, n iter ∈ N, β (0) ∈ R p // computation of β for k = 1, . . . , n iter do z (k) =β (k-1) -γ∇f (β (k-1) ) // GD step β (k) = prox γg(•,λ) z (k) // prox. step // backward computation of the gradient g v = ∇C(β (n iter ) ), h = 0 R r for k = n iter , n iter -1, . . . , 1 do h += v ∂ λ prox γg(•,λ) (z (k) ) // O(pr) v ← ∂ z prox γg(•,λ) (z (k) ) v // O(p) v ← (Id -γ∇ 2 f (β (k) ))v // O(np)
return β n iter , h

Bilevel Optimization with Non-Smooth Inner Problems

We recalled above how to compute hypergradients when the inner optimization problem is smooth. In this section we tackle the bilevel optimization Problem (2) with non-smooth inner optimization Problem (1). Handling non-smooth inner problems requires specific tools detailed in Section 3.1. We then show how to compute gradients with non-smooth inner problems using implicit differentiation (Section 3.2) or iterative differentiation (Section 3.3). In Section 3.4 we tackle the problem of approximate gradient for a non-smooth inner optimization problem. Finally, we propose in Section 3.6 an algorithm to solve the bilevel optimization Problem (2).

Theoretical Framework

Differentiability of the regularization path. Before applying first-order methods to tackle Problem (2), one must ensure that the regularization path λ → β(λ) is almost everywhere differentiable (as in Figure 2). This is the case for the Lasso [START_REF] Mairal | Complexity analysis of the lasso regularization path[END_REF] and the SVM [START_REF] Pontil | Properties of support vector machines[END_REF] since solution paths are piecewise differentiable (see Figure 2). Results for nonquadratic datafitting terms are scarcer: [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF] address the practical resolution of sparse logistic regression, but stay evasive regarding the differentiability of the regularization path. In the general case for problems of the form Problem (1), we believe it is an open question and leave it for future work. Differentiability of proximal operators. The key point to obtain an implicit differentiation formula for non-smooth inner problems is to differentiate the fixed point equation of proximal gradient descent. From a theoretical point of view, ensuring this differentiability at the optimum is non-trivial: Poliquin and Rockafellar (1996, Thm. 3.8) showed that under a twice epi-differentiability

Sparse logistic regression

Figure 2: Regularization paths (coefficient values as a function of λ), on the diabetes and breast cancer data sets for the Lasso, the elastic net and sparse logistic regression. This illustrates the weak differentiability of the paths. We used diabetes for the Lasso and the elastic net, and the 10 first features of breast cancer for the sparse logistic regression.

condition the proximal operator is differentiable at optimum. For the convergence of forward and reverse modes in the non-smooth case, one has to ensure that, after enough iterations, the updates of the algorithms become differentiable. Deledalle et al. ( 2014) justified (weak) differentiability of proximal operators as they are non-expansive. However this may not be a sufficient condition, see Bolte and Pauwels (2020a,b). In our case, we show differentiability after support identification of the algorithms: active constraints are identified after a finite number of iterations by proximal gradient descent [START_REF] Liang | Local linear convergence of forward-backward under partial smoothness[END_REF][START_REF] Vaiter | Model consistency of partly smooth regularizers[END_REF] and proximal coordinate descent, see Nutini (2018, Sec. 6.2) or [START_REF] Klopfenstein | Model identification and local linear convergence of coordinate descent[END_REF]. Once these constraints have been identified convergence is linear towards the Jacobian (see Theorem 12 andFigures 3, 10 and11).

For the rest of this paper, we consider the bilevel optimization Problem (2) with the following assumptions on the inner Problem (1).

Assumption 2 Smoothness. The function f : R p → R is a convex, differentiable function, with a L-Lipschitz gradient.

Assumption 3 Proper, closed, convex. For all λ ∈ R r , for any j ∈ [p], the function g j (•, λ) : R → R is proper, closed and convex.

Assumption 4 Non-degeneracy The problem admits at least one solution arg min

β∈R p Φ(β, λ) = ∅ ,
and, for any β solution of Problem (1), we have

-∇f ( β) ∈ ri ∂ β g( β, λ) .
To be able to extend iterative and implicit differentiation to the non-smooth case, we need to introduce the notion of generalized support.

Definition 5 Generalized support, (Nutini et al., 2019, Def. 1).

For a solution β ∈ arg min

β∈R p Φ(β, λ), its generalized support Ŝ ⊆ [p] is the set of indices j ∈ [p] such that g j is differentiable at βj Ŝ {j ∈ [p] : ∂ β g j ( βj , λ) is a singleton} .
An iterative algorithm is said to achieve finite support identification if its iterates β (k) converge to β, and there exists K ≥ 0 such that for all j / ∈ Ŝ, for all k ≥ K, β

(k) j = βj .
Examples. For the 1 norm (promoting sparsity), g j ( βj , λ) = e λ | βj |, the generalized support is Ŝ {j ∈ [p] : βj = 0}. This set corresponds to the indices of the non-zero coefficients, which is the usual support definition. For the SVM estimator, g j ( βj , λ) = ι [0,e λ ] ( βj ). This function is non-differentiable at 0 and at e λ . The generalized support for the SVM estimator then corresponds to the set of indices such that βj ∈]0, e λ [.

Finally, to prove local linear convergence of the Jacobian we assume regularity and strong convexity on the generalized support.

Assumption 6 Locally C 2 and C 3 . The map β → f (β) is locally C 3 around β. For all λ ∈ R r , for all j ∈ Ŝ the map g j (•, λ) is locally C 2 around βj . Assumption 7 Restricted injectivity.
Let β be a solution of Problem (1) and Ŝ its generalized support. The solution β satisfies the following restricted injectivity condition

∇ 2 Ŝ, Ŝ f ( β) 0 .
Assumptions 2 and 3 are classical to ensure inner problems can be solved using proximal algorithms. Assumption 4 can be seen as a generalization of constraint qualifications (Hare and Lewis, 2007, Sec. 1) and is crucial to ensure support identification. Note that Assumption 4 is hard to verify in advance in practice, and hard to relax theoretically [START_REF] Fadili | Sensitivity analysis for mirror-stratifiable convex functions[END_REF]. Assumptions 6 and 7 are classical for the analysis [START_REF] Liang | Activity identification and local linear convergence of Forward-Backward-type Methods[END_REF] and sufficient to derive rates of convergence for the Jacobian of the inner problem once the generalized support has been identified. Assumption 6 is met for usual quadratic and logistic losses, as well as for usual penalties ( 1 , 1 + 2 -squared, box constraints). For instance for the Lasso Assumption 7 boils down to X S X S 0, which holds with probability one if the entries of X are drawn from a continuous distribution [START_REF] Tibshirani | The lasso problem and uniqueness[END_REF].

The next lemma guarantees uniqueness of Problem (1) under Assumptions 4 and 7.

Lemma 8 (Liang et al., 2017, Prop. 4.1). Assume that there exists a neighborhood Λ of λ such that Assumptions 4 and 7 are satisfied for every λ ∈ Λ. Then for every λ ∈ Λ, Problem (1) has a unique solution, and the map λ → β(λ) is well-defined on Λ.

We first show how implicit and iterative differentiation can be used with a non-smooth inner problem. [START_REF] Peyré | Learning analysis sparsity priors[END_REF] proposed to smooth the inner optimization problem, [START_REF] Ochs | Bilevel optimization with nonsmooth lower level problems[END_REF]; [START_REF] Frecon | Bilevel learning of the group lasso structure[END_REF] relied on the forward-mode combined with Bregman iterations to get differentiable steps. For non-smooth optimization problems, implicit differentiation has been considered for (constrained) convex optimization problems [START_REF] Gould | On differentiating parameterized argmin and argmax problems with application to bi-level optimization[END_REF][START_REF] Amos | Optnet: Differentiable optimization as a layer in neural networks[END_REF][START_REF] Agrawal | Differentiable convex optimization layers[END_REF], Lasso-type problems (Mairal et al., 2012;[START_REF] Bertrand | Implicit differentiation of Lasso-type models for hyperparameter optimization[END_REF], total variation penalties [START_REF] Cherkaoui | Learning to solve TV regularised problems with unrolled algorithms[END_REF] and generalized to strongly monotone operators [START_REF] Winston | Neural monotone operator equilibrium networks[END_REF].

Hypergradient Computation: Implicit Differentiation

The exact proof of Theorem 1 cannot be applied when β → Φ(β, λ) is non-smooth, as Equations ( 7) and ( 8) no longer hold. Nevertheless, instead of the optimality condition of smooth optimization, Equation ( 7), one can leverage the fixed point iteration of proximal gradient descent, which we will see in Equation ( 12). The main theoretical challenge is to show the differentiability of the function β → prox γg (β -γ∇f (β)). Besides, taking advantage of the generalized sparsity of the regression coefficients β(λ) , one can show that the Jacobian Ĵ is row-sparse, leading to substantial computational benefits when computing the hypergradient ∇ λ L(λ)) for Problem (1), Theorem 9 Non-smooth implicit formula Suppose Assumptions 2, 3 and 6 hold. Let

0 < γ ≤ 1/L, where L is the Lipschitz constant of ∇f . Let λ ∈ R r , Λ be a neighborhood of λ, and Γ Λ β(λ) -γ∇f ( β(λ) ) : λ ∈ Λ .
In addition, (H1) Suppose Assumptions 4 and 7 hold on Λ.

(H2) Suppose λ → β(λ) is continuously differentiable on Λ.

(H3) Suppose for all z ∈ Γ Λ , λ → prox γg(•,λ) (z) is continuously differentiable on Λ.

(H4) Suppose ∂ z prox γg(•,λ) and ∂ λ prox γg(•,λ) are Lipschitz continuous on Γ Λ × Λ.

Let β β(λ) be the solution of Problem (1), Ŝ its generalized support of cardinality ŝ. Then the Jacobian Ĵ of the inner Problem (1) is given by the following formula, ẑ = β -γ∇f ( β), and

A Id ŝ -∂ z prox γg(•,λ) (ẑ) Ŝ Id ŝ -γ∇ 2 Ŝ, Ŝ f ( β) Ĵ Ŝc : = ∂ λ prox γg(•,λ) (ẑ) Ŝc , (10)
Ĵ Ŝ: = A -1 ∂ λ prox γg(•,λ) (ẑ) Ŝ -γ∂ z prox γg(•,λ) (ẑ) Ŝ ∇ 2 Ŝ, Ŝc f ( β) Ĵ Ŝc . ( 11
)
Proof According to Theorem 8, Assumptions 4 and 7 ensure Problem (1) has a unique minimizer and λ → β(λ) is well-defined on Λ. We consider the proximal gradient descent fixed point equation:

β(λ) = prox γg (•,λ) β(λ) -γ∇f ( β(λ) ) . ( 12
)
Together with the conclusion of Theorem 8, Assumptions 2 and 6, and given (H2), ( H3) and (H4), we have that λ

→ ψ β (λ) -γ∇f ( β(λ) ), λ prox γg(•,λ) β(λ) -γ∇f ( β(λ)
) is differentiable at λ. One can thus differentiate Equation ( 12) with respect to λ, which leads to

Ĵ = ∂ z prox γg(•,λ) (ẑ) Id -γ∇ 2 f ( β) Ĵ + ∂ λ prox γg(•,λ) (ẑ) , (13)
with ẑ = β -γ∇f ( β). In addition to 0 < γ < 1/L ≤ 1/L j , the separability of g and Assumptions 2 to 4 and 6 ensure (see Theorem 20) that for any j ∈ Ŝc ,

∂ z prox γg j (•,λ) βj -γ∇ j f ( β) = 0 . ( 14
)
Plugging Equation ( 14) into Equation ( 13) ensures Equation ( 10) for all j ∈ Ŝc

Ĵj: = ∂ λ prox γg j (•,λ) βj -γ∇ j f ( β) . ( 15
)
Plugging Equations ( 14) and ( 15) into Equation ( 13) shows that the Jacobian restricted on the generalized support Ŝ satisfies the following linear system

Id ŝ -∂ z prox γg(•,λ) (ẑ) Ŝ Id ŝ -γ∇ 2 Ŝ, Ŝ f ( β) Ĵ Ŝ: = -γ∂ z prox γg(•,λ) (ẑ) Ŝ ∇ 2 Ŝ, Ŝc f ( β) Ĵ Ŝc : + ∂ λ prox g (ẑ) Ŝ: . Since 0 < γ ≤ 1/L, ∂ z prox γg(•,λ) (ẑ) Ŝ (Id ŝ -γ∇ 2 Ŝ, Ŝ f ( β)) 2 ≤ ∂ z prox γg(•,λ) (ẑ) Ŝ • Id ŝ -γ∇ 2 Ŝ, Ŝ f ( β) 2 < 1 . ( 16
)
Since Equation ( 16) holds, A

Id ŝ -∂ z prox γg(•,λ) (ẑ) Ŝ (Id ŝ -γ∇ 2 Ŝ, Ŝ f ( β))
is invertible, which leads to Equation (11).

Remark 10 In the smooth case a p × p linear system is needed to compute the Jacobian in Equation (9). For non-smooth problems this is reduced to an ŝ × ŝ linear system (ŝ ≤ p being the size of the generalized support, e.g., the number of non-zero coefficients for the Lasso). This leads to significant speedups in practice, especially for very sparse vector β(λ) .

Remark 11 To obtain Theorem 9 we differentiated the fixed point equation of proximal gradient descent, though one could differentiate other fixed point equations (such as the one from proximal coordinate descent). The value of the Jacobian Ĵ obtained with different fixed point equations would be the same, yet the associated systems could have different numerical stability properties. We leave this analysis to future work.

Hypergradient Computation: Iterative Differentiation

Instead of implicit differentiation, it is also possible to use iterative differentiation on proximal solvers. In section Section 2 we presented forward and reverse modes differentiation of proximal gradient descent (Algorithms 1 and 2). In this section we study the iterative differentiation of proximal coordinate descent (Algorithms 3 and 4). To instantiate algorithms easily on problems such as the Lasso, partial derivatives of usual proximal operators can be found in Table 3.

For coordinate descent, the computation of the iterative Jacobian in a forward way involves differentiating the following update

z j ← β j -γ j ∇ j f (β) β j ← prox γ j g j (β j -γ j ∇ j f (β)) J j: ← ∂ z prox γ j g j (•,λ) (z j ) ∈R J j: -γ j ∇ 2 j: f (β)J ∈R p + ∂ λ prox γ j g j (•,λ) (z j ) ∈R p
.

We address now the convergence of the iterative Jacobian scheme, a question which remained open in Deledalle et al. (2014, Section 4.1). We show next that the forward-mode converges to the Jacobian in the non-smooth separable setting of this paper. Moreover, we prove that the iterative Jacobian convergence is locally linear after support identification.

Theorem 12 Local linear convergence of the Jacobian. Let 0 < γ ≤ 1/L. Suppose Assumptions 2, 3 and 6 hold. Let λ ∈ R r , Λ be a neighborhood of λ, and Γ Λ β(λ) -γ∇f ( β(λ) ) : λ ∈ Λ . In addition, suppose hypotheses (H1) to (H4) from Theorem 9 are satisfied and the sequence (β (k) ) k∈N generated by Algorithm 1 (respectively by Algorithm 3) converges toward β.

Then, the sequence of Jacobians (J (k) ) k≥0 generated by the forward-mode differentiation of proximal gradient descent (Algorithm 1) (respectively by forward-mode differentiation of proximal coordinate descent, Algorithm 3) converges locally linearly towards Ĵ .

Proof of Theorem 12 can be found in Appendix C. Comments on Figure 3. We illustrate the results of Theorem 12 on SVM (for the Lasso and sparse logistic regression, see Figures 10 and11 in Appendix D) for multiple data sets (leukemia, rcv1, news20 and real-sim 4 ). The values of the hyperparameters λ are summarized in Table 6. Regression coefficients β(λ) were computed to machine precision (up to duality gap smaller than 10 -16 ) using a state-of-the-art coordinate descent solver implemented in Lightning [START_REF] Blondel | Lightning: large-scale linear classification, regression and ranking in python[END_REF]. The exact Jacobian was computed via implicit differentiation (Equation ( 11)). Once these quantities were obtained, we used the forward-mode differentiation of proximal coordinate descent (Algorithm 3) and monitored the distance between the iterates of the regression coefficients β (k) and the exact solution 4. Data available on the libsvm website: https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/ Algorithm 3 Forward-mode PCD input :

X ∈ R n×p , y ∈ R n , λ ∈ R r , n iter ∈ N, β ∈ R p , J ∈ R p×r ,
z j ← β j -γ j ∇ j f (β) // CD step dz j ← J j: -γ j ∇ 2 j: f (β)J β j ← prox γ j g j (•,λ) (z j ) // prox. step
// update the Jacobian // diff. with respect to λ

J j: ← ∂ z prox γ j g j (•,λ) (z j )dz j J j: += ∂ λ prox γ j g j (•,λ) (z j ) β (k) = β J (k) = J v = ∇C(β) return β n iter , J v Algorithm 4 Reverse-mode PCD input : X ∈ R n×p , y ∈ R n , λ ∈ R r , n iter ∈ N, β ∈ R p , γ 1 , . . . , γ p // compute coef.
for k = 1, . . . , n iter do for j = 1, . . . , p do // update the regression coefficients β. We also monitored the distance between the iterates of the Jacobian J (k) and the exact Jacobian Ĵ . The red vertical dashed line represents the iteration number where support identification happens. Once the support is identified, Figures 3, 10 and 11 illustrate the linear convergence of the Jacobian. However, the behavior of the iterative Jacobian before support identification is more erratic and not even monotone.

z j ← β j -γ j ∇ j f (β) // CD step β j ← prox γ j g j (•,λ) (z j ) // prox. step β (k,j) = β; z (k) j = z j // store iterates // compute gradient g in a backward way v = ∇C(β n iter ), h = 0 R r for k = n iter , n iter -1, . . . , 1 do for j = p, . . . , 1 do h -= γ j v j ∂ λ prox γ j g j (•,λ) z (k) j v j * = ∂ z prox γ j g j (•,λ) z (k) j v -= γ j v j ∇ 2 j: f (β (k,j) ) // O(np) return β n iter , h

Generalized support identification

Hypergradient Computation with Approximate Gradients

As mentioned in Section 2, relying on iterative algorithms to solve Problem (1), one only has access to an approximation of β(λ) : this may lead to numerical errors when computing the gradient in Theorem 9. Extending the result of Pedregosa (2016, Thm. 1), which states that hypergradients can be computed approximately, we give a stability result for the computation of approximate hypergradients in the case of non-smooth inner problems.

For this purpose we need to add several assumptions to the previous framework.

Theorem 13 Bound on the error of approximate hypergradient.

For λ ∈ R r , let β(λ) ∈ R p be the exact solution of the inner Problem (1), and Ŝ its generalized support. Suppose Assumptions 2, 3 and 6 hold. Let Λ be a neighborhood of λ, and Γ Λ β(λ) -γ∇f ( β(λ) ) : λ ∈ Λ . Suppose hypotheses (H1) to (H4) from Theorem 9 are satisfied. In addition suppose (H5) The application β → ∇ 2 f (β) is Lipschitz continuous.

(H6) The criterion β → ∇C(β) is Lipschitz continuous.

(H7) Both optimization problems in Algorithm 5 are solved up to precision with support identification: β (λ) -β(λ) ≤ , A is invertible, and

A -1 ∇ Ŝ C(β (λ) ) -v ≤ .
Then the error on the approximate hypergradient h returned by Algorithm 5 is of the order of magnitude of the error on β (λ) and v

∇L(λ) -h = O( ) .
Proof of Theorem 13 can be found in Appendix C.1. Following the analysis of Pedregosa (2016), two sources of approximation errors arise when computing the hypergradient: one from the inexact computation of β, and another from the approximate resolution of the linear system. Theorem 13 states that if the inner optimization problem and the linear system are solved up to precision , i.e., β(λ)β (λ) ≤ and A -1 ∇ S C(β (λ) )v ≤ , then the approximation on the hypergradient is also of the order of .

Remark 14 The Lipschitz continuity of the proximity operator with respect to λ (H4) is satisfied for usual proximal operators, in particular all the operators in Table 3. The Lipschitz continuity of the Hessian and the criterion, hypotheses (H5) and (H6), are satisfied for usual machine learning loss functions and criteria, such as the least squares and the logistic loss.

Remark 15 To simplify the analysis, we used the same tolerance for the resolution of the inner Problem (1) and the resolution of the linear system. Theorem 13 gives intuition on the fact that the inner problem does not need to be solved at high precision to lead to good hypergradients estimation. Note that in practice one does not easily control the distance between the approximate solution and the exact one β (k) -β : most softwares provide a solution up to a given duality gap (sometimes even other criteria), not β (k) -β .

Proposed Method for Hypergradient Computation

We now describe our proposed method to compute the hypergradient of Problem (2). In order to take advantage of the sparsity induced by the generalized support, we propose an implicit differentiation algorithm for non-smooth inner problem that can be found in Algorithm 5. First, we compute a solution of the inner Problem (1) using a solver identifying the generalized support [START_REF] Liang | Local linear convergence of forward-backward under partial smoothness[END_REF][START_REF] Klopfenstein | Model identification and local linear convergence of coordinate descent[END_REF]. Then, the hypergradient is computed by solving the linear system in Equation ( 11). This linear system, as mentioned in Section 2, can be solved using multiple algorithms, including conjugate gradient or fixed point methods. Table 4 summarizes the computational complexity in space and time of the described algorithms. Table 4: Cost in time and space for each method: p is the number of features, n the number of samples, r the number of hyperparameters, and ŝ is the size of the generalized support (Theorem 5, ŝ ≤ p and usually ŝ p). The number of iterations of the inner solver is noted n iter , the number of iterations of the solver of the linear system is noted n sys . From a practical point of view, once the hypergradient has been computed, first-order methods require the definition of a step size to solve the non-convex Problem (2). As the Lipschitz constant is not available for the outer problem, first-order methods need to rely on other strategies, such as:

Differentiation

• Gradient descent with manually adjusted fixed step sizes [START_REF] Frecon | Bilevel learning of the group lasso structure[END_REF][START_REF] Ji | Provably faster algorithms for bilevel optimization and applications to meta-learning[END_REF]. The main disadvantage of this technique is that it requires a careful tuning of the step size for each experiment. In addition to being potentially tedious, it does not lead to an automatic procedure.

• L-BFGS (as in [START_REF] Deledalle | Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection[END_REF]. L-BFGS is a quasi-Newton algorithm that exploits past iterates to approximate the Hessian and propose a better descent direction, which is combined with some line search [START_REF] Nocedal | Numerical optimization[END_REF]). Yet, due to the approximate gradient computation, we observed that L-BFGS did not always converge.

• ADAM [START_REF] Kingma | A method for stochastic optimization[END_REF]. It turned out to be inappropriate to the present setting. ADAM was very sensitive to the initial step size and required a careful tuning for each experiment.

• Iteration specific step sizes obtained by line search [START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF]. While the approach from Pedregosa (2016) requires no tuning, we observed that it could diverge when close to the optimum. The normalized gradient strategy (Watt et al., 2020, Sec. 3.9)5 proposed in Algorithm 6, used in all the experiments, turned out to be robust and efficient across problems and data sets.

Remark 16 Uniqueness. The solution of Problem (1) may be non-unique, leading to a multi-valued regularization path λ → β(λ) [START_REF] Liu | A generic first-order algorithmic framework for bi-level programming beyond lower-level singleton[END_REF] and requiring tools such as optimistic gradient (Dempe et al., 2015, Chap. 3.8). Though it is not possible to ensure uniqueness in practice, we did not face experimental issues due to potential non-uniqueness.

For the Lasso, this experimental observation can be theoretically justified [START_REF] Tibshirani | The lasso problem and uniqueness[END_REF]: when the design matrix is sampled from a continuous distribution, the solution of the Lasso is almost surely unique.

Remark 17 Initialization and warm start. One advantage of the non-smooth case with the 1 norm is that one can find a good initialization point: there exists a value λ max (see Table 1) such that the solution of Problem (1) vanishes for λ ≥ λ max . Hence, a convenient and robust initialization value can be chosen as e λ = e λmax /100. This is in contrast with the smooth case, where finding a good initialization heuristic is hard: starting in flat zones can lead to poor performance for gradient-based methods [START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF]. Algorithm 5 is called multiple times in Algorithm 6: several inner optimization problems and linear systems which are "similar" are solved successively. That is why we use warm-start to solve these problems.

Remark 18 Role of the step size γ in Algorithm 5. In all the convex penalties we used ( 1 -norm, 1 + 2 -squared norm, indicator function) the step size γ simplifies and does not appear in the implicit differentiation formula. Instantiations of Algorithm 5 for the Lasso, the elastic net, the weighted Lasso and the dual of the SVM can be found in Appendix A.

Algorithm 5 Implicit differentiation input : λ ∈ R, > 0 init : γ > 0 // compute the solution of inner problem

Find β such that:

Φ(β, λ) -Φ( β, λ) ≤ // compute the gradient Compute the generalized support S of β, z = β -γ∇f (β) J S c : = ∂ λ prox γg(•,λ) (z) S c s = |S| A=Id s -∂ z prox γg(•,λ) (z) S (Id s -γ∇ 2 S,S f (β)) Find v ∈ R s s.t. A -1 ∇ S C(β) -v ≤ B = ∂ λ prox γg(•,λ) (z) S -γ∂ z prox γg(•,λ) (z) S ∇ 2 S,S c f (β)J S c ∇L(λ) = J S c : ∇ S c C(β) + v B return L(λ) C(β), ∇L(λ)
Algorithm 6 Gradient descent with approximate gradient input : λ ∈ R r , ( i ) init : use_adaptive_step_size = True for i = 1, . . . , iter do λ old ← λ // compute the value and the gradient

L(λ), ∇L(λ) ← Algorithm 5(X, y, λ, i ) if use_adaptive_step_size then α = 1/ ∇L(λ) λ -= α∇L(λ) // gradient step if L(λ) > L(λ old ) then use_adaptive_step_size = False α /= 10 return λ

Experiments

In this section, we illustrate the benefits of our proposed Algorithm 5 to compute hypergradients and Algorithm 6 to solve Problem (2). Our package, sparse-ho, is implemented in Python. It relies on Numpy [START_REF] Harris | Array programming with NumPy[END_REF], Numba [START_REF] Lam | Numba: A LLVM-based Python JIT Compiler[END_REF] and SciPy [START_REF] Virtanen | Scipy 1.0: fundamental algorithms for scientific computing in python[END_REF]. Figures were plotted using matplotlib [START_REF] Hunter | Matplotlib: A 2d graphics environment[END_REF]. The package is available under BSD3 license at https://github.com/qb3/sparse-ho, with documentation and examples available at https://qb3.github.io/sparse-ho/. Online code includes scripts to reproduce all figures and experiments of the paper.

Hypergradient computation

Comparison with alternative approaches (Figure 4). First, we compare different methods to compute the hypergradient:

• Forward-mode differentiation of proximal coordinate descent (Algorithm 3).

• Reverse-mode differentiation of proximal coordinate descent (Algorithm 4).

• cvxpylayers [START_REF] Agrawal | Differentiable convex optimization layers[END_REF], a software based on cvxpy [START_REF] Diamond | CVXPY: A Python-embedded modeling language for convex optimization[END_REF], solving disciplined parametrized programming and providing derivatives with respect to the parameters of the program. It is thus possible to use cvxpylayers to compute gradients with respect to the regularization parameters.

Figure 4 compares the time taken by multiple methods to compute a single hypergradient ∇L(λ) for the Lasso (see Table 1), for multiple values of λ. It shows the time taken to compute the regression coefficients and the hypergradient, as a function of the number of columns, sampled from the design matrix from the gina data set. The columns were selected at random and 10 repetitions were performed for each point of the curves. In order to aim for good numerical precision, problems were solved up to a duality gap of 10 -6 for the forward-mode and the reverse-mode. cvxpylayers relies on cvxpy, solving Problem (1) using a splitting conic solver [START_REF] O'donoghue | SCS: Splitting conic solver[END_REF]. Since the termination criterion of the splitting conic solver is not exactly the duality gap (O'Donoghue et al., 2016, Sec. 3.5), we used the default tolerance of 10 -4 . The hypergradient ∇L(λ) was computed for hold-out mean squared error (see Table 2). The forward-mode differentiation of proximal coordinate descent is one order of magnitude faster than cvxpylayers and two orders of magnitude faster than the reverse-mode Implicit diff.

Implicit diff. + Celer Forward-mode PCD Time (s)

Figure 5: Lasso with hold-out criterion: absolute difference between the exact hypergradient (using β) and the iterate hypergradient (using β (k) ) of the Lasso as a function of time.

Results are for three data sets and two different regularization parameters. "Implicit diff. + Celer)" uses Celer [START_REF] Massias | Dual extrapolation for sparse generalized linear models[END_REF] instead of our proximal coordinate descent implementation.

differentiation of proximal coordinate descent. The larger the value of λ, the sparser the coefficients β are, leading to significant speedups in this regime. This performance is in accordance with the lower time cost of the forward mode in Table 4.

Combining implicit differentiation with state-of-the art solvers (Figures 5 and6). We now compare the different approaches described in Section 3:

• Forward-mode differentiation of proximal coordinate descent (Algorithm 3).

• Implicit differentiation (Algorithm 5) with proximal coordinate descent to solve the inner problem. For efficiency, this solver was coded in Numba [START_REF] Lam | Numba: A LLVM-based Python JIT Compiler[END_REF].

• Implicit differentiation (Algorithm 5) with state-of-the-art algorithm to solve the inner problem: we used Celer [START_REF] Massias | Dual extrapolation for sparse generalized linear models[END_REF] for the Lasso, and Lightning [START_REF] Blondel | Lightning: large-scale linear classification, regression and ranking in python[END_REF] for the SVM.

Figure 5 shows for three data sets and two values of regularization parameters the absolute difference between the exact hypergradient and the approximate hypergradient obtained via multiple algorithms as a function of time. Figure 6 reports similar results for the SVM, on the same data sets, except news20, which is not well suited for SVM, due to limited number of samples. Figure 6: SVM with hold-out criterion: absolute difference between the exact hypergradient (using β) and the iterate hypergradient (using β (k) ) of the SVM as a function of time. "Implicit diff. + Lightning" uses Lightning [START_REF] Blondel | Lightning: large-scale linear classification, regression and ranking in python[END_REF], instead of our proximal coordinate descent implementation.

First, it demonstrates that implicit differentiation methods are faster than the forwardmode of proximal coordinate descent (pink). This illustrates the benefits of restricting the gradient computation to the support of the Jacobian, as described in Section 3.5. Second, thanks to the flexibility of our approach, we obtain additional speed-ups by combining implicit differentiation with a state-of-the-art solver, Celer. The resulting method (orange) significantly improves over implicit differentiation using a vanilla proximal coordinate descent (green).

Resolution of the Bilevel Optimization Problem

In this section we compare multiple methods to find the optimal hyperparameters for the Lasso, elastic net and multiclass sparse logistic regression. The following methods are compared:

• Grid-search: for the Lasso and the elastic net, the number of hyperparameters is small, and grid-search is tractable. For the Lasso we chose a grid of 100 hyperparameters λ, uniformly spaced between λ maxln(10 4 ) and λ max . For the elastic net we chose for each of the two hyperparameters a grid of 10 values uniformly spaced between λ max and λ maxln(10 4 ). The product grid thus has 10 2 points.

• Random-search: we chose 30 values of λ sampled uniformly between λ max and λ maxln(10 4 ) for each hyperparameter. For the elastic net we chose 30 points sampled uniformly in [λ maxln(10 4 ), λ max ] × [λ maxln(10 4 ), λ max ].

• SMBO: this algorithm is SMBO using as criterion expected improvement (EI) and the Tree-structured Parzen Estimator (TPE) as model. First it evaluates L using 5 values of λ, chosen uniformly at random between λ max and λ maxln(10 4 ). Then a TPE model is fitted on the data points (λ (1) , L(λ (1) )), . . . , (λ (5) , L(λ (5) )). Iteratively, the EI is used to choose the next point to evaluate L at, and this value is used to update the model. We used the hyperopt implementation [START_REF] Bergstra | Hyperopt: A python library for optimizing the hyperparameters of machine learning algorithms[END_REF].

• 1st order : first-order method with exact gradient (Algorithm 6 with constant tolerances i = 10 -6 ), with λ maxln(10 2 ) as a starting point.

• 1st order approx : a first-order method using approximate gradient (Algorithm 6 with tolerances i , geometrically decreasing from 10 -2 to 10 -6 ), with λ maxln(10 2 ) as a starting point.

Outer criterion. In the Lasso and elastic net experiments, we pick a K-fold CV loss as outer criterion 6 . Hence, the data set (X, y) is partitioned into K hold-out data sets (X train k , y train k ), (X val k , y val k ). The bilevel optimization problems then write

arg min λ=(λ 1 ,λ 2 )∈R 2 L(λ) = 1 K K k=1 y val k -X val k β(λ,k) 2 2 s.t. β(λ,k) ∈ arg min β∈R p 1 2n y train k -X train k β 2 2 + e λ 1 β 1 + e λ 2 2 β 2 2 , ∀k ∈ [K] , (17)
while Lasso CV is obtained taking λ 2 → -∞ in the former. By considering an extended variable β ∈ R K×p , cross-validation can be cast as an instance of Problem (2). Figure 7 represents the cross-validation loss in Lasso CV as a function of the regularization parameter λ (black curve, three top rows) and as a function of time (bottom). Each point corresponds to the evaluation of the cross-validation criterion for one λ value. The top rows show cross-validation loss as a function of λ, for the grid-search, the SMBO optimizer and the first-order method. The lightest crosses correspond to the first iterations of the algorithm and the darkest, to the last ones. For instance, Lasso grid-search starts to evaluate the cross-validation function with λ = λ max and then decreases to λ = λ maxln(10 4 ). On all the data sets, first-order methods are faster to find the optimal regularization parameter, requiring only 5 iterations.

Figure 8 represents the level sets of the cross-validation loss for the elastic net (three top rows) and the cross-validation loss as a function of time (bottom). One can see that after 5 iterations the SMBO algorithm (blue crosses) suddenly slows down (bottom) as the hyperparameter suggested by the algorithm leads to a costly optimization problem to solve, while first-order methods converge quickly as for Lasso CV. In the present context, inner problems are slower to solve for low values of the regularization parameters.

Multiclass sparse logistic regression (# classes hyperparameters, Figure 9). We consider a multiclass classification problem with q classes. The design matrix is noted X ∈ R n×p , and the target variable y ∈ {1, . . . , q} n . We chose to use a one-versus-all model with q 6. In our experiments the default choice is K = 5. regularization parameters. We use a binary cross-entropy for the inner loss

ψ k (β, λ k ; X, y) - 1 n n i=1 (1 y i =k ln(σ(X i: β)) + (1 -1 y i =k ) ln(1 -σ(X i: β))) + e λ k β 1 ,
and a multiclass cross-entropy for the outer criterion C β(λ 1 ) , . . . , β(λq) ; X, y -

n i=1 q k=1 ln e X i: β(λ k )
q l=1 e X i: β(λ l ) 1 y i =k .

(

) 18
With a single train/test split, the bilevel problem to solve writes:

arg min λ (λ 1 ,...,λq)∈R q C β(λ 1 ) , . . . , β(λq) ; X test , y test s.t. β(λ k ) ∈ arg min

β∈R p ψ k (β, λ k ; X train , y train ) ∀k ∈ [q] . (19)
Figure 9 represents the multiclass cross-entropy (top), the accuracy on the validation set (middle) and the accuracy on the test set (unseen data, bottom). When the number of hyperparameter is moderate (q = 10, on mnist and usps), the multiclass cross-entropy reached by SMBO and random techniques is as good as first-order techniques. This is expected and follows the same conclusion as [START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF]; [START_REF] Frazier | A tutorial on Bayesian optimization[END_REF]: when the number of hyperparameters is moderate, SMBO and random techniques can be used efficiently. However, when the number of hyperparameters increases (rcv1, q = 53 and aloi, q = 1000), the hyperparameter space is too large: zero-order solvers simply fail. On the contrary, first-order techniques manage to find hyperparameters leading to significantly better accuracy.

Remark 19 On the data used in Figure 9, the model with one hyperparameter per class did not yield significantly better test accuracy compared to a multiclass logistic regression with only one regularization hyperparameter for all the classes. This may mean that the model with one hyperparameter per class is not well suited for this data. It can also be due to the fact that in this case, the bilevel optimization problem becomes highly non-convex, and only converges toward a poor local minima. We want to emphasize that we provide an efficient way to compute the hypergradient ∇ λ L(λ). Besides, to our knowledge, the perfect resolution of the full bilevel optimization problem with a non-smooth inner problem remains an open question.

Conclusion

In this work we considered the problem of hyperparameter optimization to select the regularization parameter of linear models with non-smooth objective. Casting this problem as a bilevel optimization problem, we proposed to use first-order methods. We showed that the usual automatic differentiation techniques, implicit differentiation, forward and reverse modes, can be used to compute the hypergradient, despite the non-smoothness of the inner problem. Experimentally, we showed the interest of first-order techniques to solve bilevel optimization on a wide range of estimators ( 1 penalized methods, SVM, etc.) and data sets. The presented techniques could also be extended to more general bilevel optimization problems, in particular implicit differentiation could be well suited for meta-learning problems, with a potentially large number of hyperparameters. Another important future direction would be to extend the work on stochastic hypergradients [START_REF] Grazzi | Convergence properties of stochastic hypergradients[END_REF] in the non-smooth case.

1st-order 1st-order approx Grid-search Random-search SMBO Then, for all j ∈ Ŝ, the map β → prox γ j g j (•,λ) is differentiable at β Ŝ . Moreover, for all j ∈ Ŝc , prox γ j g j (•,λ) is constant around βjγ j ∇ j f ( β). Thus, β → prox γ j g j (•,λ) (β jγ j ∇ j f (β))

is differentiable at β with gradient 0.

B.2 Linear Convergence

We now detail the following result: an asymptotic vector autoregressive sequence, with an error term vanishing linearly to 0, converges linearly to its limit. In a more formal way:

Lemma 21 Let A ∈ R p×p , b ∈ R p with ρ(A) < 1. Let (J (k) ) k∈N be a sequence of R p such that J (k+1) = AJ (k) + b + (k) , (20)
with ( (k) ) k∈N a sequence which converges linearly to 0, then (J (k) ) k∈N converges linearly to its limit Ĵ (Id -A) -1 b.

Proof Assume ( (k) ) k∈N converges linearly. Then, there exists c 1 > 0, 0 < ν < 1 such that

(k) ≤ c 1 ν k .
Applying a standard result on spectral norms (see [START_REF] Polyak | Introduction to optimization[END_REF], Chapter 2, Lemma 1) yields a bound on A k 2 . More precisely, for every δ > 0 there is a constant c 2 (δ) = c 2 such that

A k 2 ≤ c 2 (ρ(A) + δ) k .
Without loss of generality, we consider from now on a choice of δ such that ρ(A) + δ < 1.

Since Ĵ = (Id -A) -1 b the limit Ĵ of the sequence satisfies

Ĵ = A Ĵ + b . (21)
Taking the difference between Equations ( 20) and ( 21) yields:

J (k+1) -Ĵ = A(J (k) -Ĵ ) + (k) . (22)
Unrolling Equation ( 22) yields J -k ) . Taking the norm on both sides and using the triangle inequality leads to

(k+1) -Ĵ = A k+1 (J (0) -J ) + k k =0 A k (k
J (k+1) -Ĵ 2 ≤ A k+1 (J (0) -J ) 2 + k k =0 A k 2 (k-k ) ≤ A k+1 2 • J (0) -Ĵ 2 + c 1 k k =0 A k 2 • ν k-k ≤ c 2 (ρ(A) + δ) k+1 • J (0) -Ĵ 2 + c 1 k k =0 c 2 (ρ(A) + δ) k ν k-k
Assumptions 2 to 4 and 6 and the convergence of (β (k) ) toward β ensure proximal gradient descent algorithm has finite identification property (Liang et al., 2014, Thm. 3.1): we note K the iteration when identification is achieved. As before, the separability of g, Assumptions 2 to 4 and 6 ensure (see Theorem 20) ∂ z prox γg(•,λ) (z k ) Ŝc = 0, for all k ≥ K. Thus, for all k ≥ K,

J (k) Ŝc : = Ĵ Ŝc : = ∂ λ prox γg(•,λ) (z (k) ) Ŝc : .
The updates of the Jacobian then become

J (k+1) Ŝ: = ∂ z prox γg(•,λ) (z (k) ) Ŝ Id -γ∇ 2 Ŝ, Ŝ f (β (k) ) J (k) Ŝ: + ∂ λ prox γg(•,λ) (z (k) ) Ŝ: .
From Assumption 6, we have that f is locally

C 3 at β, g(•, λ) is locally C 2 at β hence prox g(•,λ) is locally C 2 . The function β → ∂ z prox γg(•,λ) (β -γ∇f (β)) Ŝ (Id -γ∇ 2 Ŝ, Ŝ f (β)) is differentiable at β. Using (H4) we have that β → ∂ λ prox γg(•,λ) (β -γ∇f (β)) Ŝ:
is also differentiable at β. Using the Taylor expansion of the previous functions yields

J (k+1) Ŝ: = ∂ z prox γg(•,λ) (ẑ) Ŝ Id -γ∇ 2 Ŝ, Ŝ f ( β) A J (k) Ŝ: + ∂ λ prox γg(•,λ) (ẑ) Ŝ: b + o( β (k) -β ) (k)
.

Thus, for 0

< γ ≤ 1/L, ρ(A) ≤ A 2 ≤ ∂ z prox γg(•,λ) (ẑ) Ŝ ≤1 (non-expansiveness) • Id -γ∇ 2 Ŝ, Ŝ f ( β) 2 <1 (Assumption 7 and 0 < γ ≤ 1/L) < 1 . ( 26
)
The inequality on the derivative of the proximal operator comes from the non-expansiveness of proximal operators. The second inequality comes from Assumption 7 and 0 < γ ≤ 1/L. Assumptions 2 to 4, 6 and 7 and the convergence of (β (k) ) toward β ensure (β (k) ) k∈N converges locally linearly (Liang et al., 2014, Thm. 3.1). The asymptotic autoregressive sequence in Equation (25), ρ(A) < 1, and the local linear convergence of ( (k) ) k∈N , yield our result using Theorem 21.

We now prove Theorem 12 for proximal coordinate descent. Proximal coordinate descent. Compared to proximal gradient descent, the analysis of coordinate descent requires studying functions defined as a the composition of p applications, each of them only modifying one coordinate.

Coordinate descent updates read as follows

β (k,j) j = prox γ j g j (•,λ) β (k,j-1) j -γ j ∇ j f (β (k,j-1) ) z (k,j-1) j . ( 27
)
We consider the following sequence J (k,j) j:

= ∂ z prox γ j g j (•,λ) (z (k,j-1) j

) J

(k,j-1) j:

γ j ∇ 2 j: f (β (k,j-1) )J (k,j-1)

+ ∂ λ prox γ j g j (•,λ) (z (k,j-1) j
) .

Note that if prox γg(•,λ) is not differentiable with respect to the first variable at z (k) (respectively with respect to the second variable λ), any weak Jacobian can be used. When (H3) holds, differentiating Equation ( 27) with respect to λ yields exactly Equation (28). Assumptions 2 to 4 and 6 and the convergence of (β (k) ) k∈N toward β ensure proximal coordinate descent has finite identification property (Klopfenstein et al., 2020, Thm. 1): we note K the iteration when identification is achieved. Once the generalized support Ŝ (of cardinality ŝ) has been identified, we have that for all k ≥ K, β

Ŝc = β Ŝ and for any j ∈ Ŝc , ∂ z prox γ j g j (•,λ) (z

(k,j-1) j ) = 0. Thus J (k,j) j: = ∂ λ prox γ j g j (•,λ) (z (k,j-1) j
). Then, we have that for any j ∈ Ŝ and for all k ≥ K:

J (k,j) j: = ∂ z prox γ j g j (•,λ) (z (k,j-1) j ) J (k,j-1) j: -γ j ∇ 2 j, Ŝ f (β (k,j-1) )J (k,j-1) Ŝ: + ∂ λ prox γ j g j (•,λ) (z (k,j-1) j ) -γ j ∂ z prox γ j g j (•,λ) (z (k,j-1) j )∇ 2 j, Ŝc f (β (k,j-1) )J (k,j -1) Ŝc :
. Let e 1 , . . . , e ŝ be the vectors of the canonical basis of R ŝ. We can consider the applications

R p → R ŝ β → ∂ z prox γ j g j (•,λ) (β j -γ j ∇ j f (β)) e j -γ j ∇ 2 j, Ŝ f (β) , and
R p → R ŝ×r β → ∂ λ prox γ j g j (•,λ) (β j -γ j ∇ j f (β)) -γ j ∂ z prox γ j g j (•,λ) (β j -γ j ∇ j f (β)) ∇ 2 j, Ŝc f (β) Ĵ Ŝc : ,
which are both differentiable at β using Assumption 6 and (H4). The Taylor expansion of the previous functions yields:

J (k,j) j: = ∂ z prox γ j g j (•,λ) (ẑ j ) e j -γ j ∇ 2 j, Ŝ f ( β) J (k,j-1) Ŝ: + ∂ λ prox γ j g j (•,λ) (ẑ j ) -γ j ∂ z prox γ j g j (•,λ) (ẑ j ) ∇ 2 j, Ŝc f ( β)J (k,j-1) Ŝc :
+ o(||β (k,j-1) -β||) .

Let j 1 , . . . , j ŝ be the indices of the generalized support of β. When considering a full epoch of coordinate descent, the Jacobian is obtained as the product of matrices of the form A s = e 1 . . . e s-1 v js e s+1 . . . e ŝ ∈ R ŝ×ŝ ,

• Approximation errors from the inexact computation of β. Dropping the dependency with respect to λ, we denote β the approximate solution and suppose the problem is solved to precision with support identification (H7)

β Ŝc = β Ŝc β Ŝ -β Ŝ ≤ .
• Approximation errors from the approximate resolution of the linear system, using (H7) yields:

A -1 ∇ Ŝ C(β) -v ≤ .
The exact solution of the exact linear system v satisfies

v = Â-1 ∇ Ŝ C( β) , with A Id | Ŝ| -∂ z prox γg(•,λ) (β -γ∇f (β)) Ŝ C Id | Ŝ| -γ∇ 2 Ŝ, Ŝ f (β) D , Â Id | Ŝ| -∂ z prox γg(•,λ) β -γ∇f ( β) Ŝ Ĉ Id | Ŝ| -γ∇ 2 Ŝ, Ŝ f ( β) D .
• Using the last two points, the goal is to bound the difference between the exact hypergradient and the approximate hypergradient, ∇L(λ)h . Following Algorithm 5, the exact hypergradient reads -γ∂ z prox γg(•,λ) β -γ∇f ( β) Ŝ ∇ 2 Ŝ, Ŝc f ( β) Ĵ Ŝc : .

We can exploit these decompositions to bound the difference between the exact hypergradient and the approximate hypergradient

Let ṽ be the exact solution of the approximate system A ṽ ∇ Ŝ C(β). The following conditions are met:

• v is the exact solution of the exact linear system and ṽ is the exact solution of the approximate linear system

 v ∇ Ŝ C( β)
A ṽ ∇ Ŝ C(β) .

• One can control the difference between the exact matrix in the linear system  and the approximate matrix

A A -Â 2 ≤ δ β -β ,
for a certain δ > 0 (Equation ( 31)).

• One can control the difference between the two right-hand side of the linear systems Figures 10 and 11 are the counterparts of Figure 3 for the Lasso and sparse logistic regression. It shows the local linear convergence of the Jacobian for the Lasso, obtained by the forward-mode differentiation of coordinate descent. The solvers used to determine the exact solution up to machine precision are Celer [START_REF] Massias | Celer: a fast solver for the lasso with dual extrapolation[END_REF][START_REF] Massias | Dual extrapolation for sparse generalized linear models[END_REF] for the Lasso and Blitz [START_REF] Johnson | Blitz: A principled meta-algorithm for scaling sparse optimization[END_REF] for the sparse logistic regression. Table 6 summarizes the values of the hyperparameters λ used in Figures 3, 10 and 11.

∇ Ŝ C(β) -∇ Ŝ C( β) ≤ L C β -β ,

D.2 Hypergradient Computation Time

The experimental setting for Figure 12 is the same as for Figure 4, but with a weighted Lasso (i. of features. The regularization amounts were chosen uniformly at random in the interval [0, e λ max ] and each point represents 10 repetitions. Figure 12 shows that when the number of hyperparameters is large the implicit differentiation outperforms the reverse-mode, ans the reverse-mode outperforms the forward-mode, by one or more orders of magnitude. This corroborates the complexities summarized in Table 4.

Figure 3 :

 3 Figure 3: Local linear convergence of the Jacobian for the SVM. Distance to optimum for the coefficients β (top) and the Jacobian J (bottom) of the forward-mode differentiation of proximal coordinate descent (Algorithm 3) on multiple data sets. One epoch corresponds to one pass over the data, i.e., one iteration with proximal gradient descent.

Figure 4 :

 4 Figure 4: Lasso with hold-out criterion: time comparison on the gina data set to compute a single hypergradient as a function of the number of features, for two values of λ, e λ = e λmax /10 (left) and e λ = e λmax /100 (right).

Figure 7 :

 7 Figure 7: Lasso with cross-validation criterion: cross-validation loss as a function of λ (black line, top) and as a function of time (bottom). Lighter markers correspond to earlier iterations of the algorithm.

Figure 8 :

 8 Figure 8: Elastic net cross-validation, time comparison (2 hyperparameters).Level sets of the cross-validation loss (black lines, top) and cross-validation loss as a function of time (bottom) on rcv1, real-sim and news20 data sets.

  ∇L(λ) = Bv + Ĵ Ŝc : ∇ Ŝc C( β) ,and similarly for the approximate versionsh = Bv + J Ŝc : ∇ Ŝc C(β) , with B ∂ λ prox γg(•,λ) (β -γ∇f (β)) Ŝ: -γ∂ z prox γg(•,λ) (β -γ∇f (β)) Ŝ ∇ 2 Ŝ, Ŝc f (β) Ĵ Ŝc : B ∂ λ prox γg(•,λ) β -γ∇f ( β) Ŝ:

  ∇L(λ) -h = Bv -Bv + Ĵ Ŝc : ∇ Ŝc C( β) -Ĵ Ŝc : ∇ Ŝc C(β) ≤ Bv -Bv + Ĵ Ŝc : ∇ Ŝc C( β) -Ĵ Ŝc : ∇ Ŝc C(β) ≤ Bv -Bv + Bv -Bv + Ĵ Ŝc : (∇ Ŝc C( β) -∇ Ŝc C(β)) ≤ v • B -B + B • vv + L C Ĵ Ŝc : • β -β . (30)Bounding vv and B -B in Equation (30) yields the desired result which is bounding the difference between the exact hypergradient and the approximate hypergradient∇L(λ)h . Bound on vv . We first prove that A -Â = O( ).Let L H be the Lipschitz constant of the application β → ∇ 2 f (β), then we haveA -Â 2 = CD -Ĉ D 2 ≤ CD -C D 2 + C D -Ĉ D 2 β ) using (H4) ≤ L H β -β + O( β -β ) = O( β -β ) .

  since β → ∇C(β) is L C -Lipschitz continuous (H6).

Figure 10 :

 10 Figure 10: Local linear convergence of the Jacobian for the Lasso. Distance to optimum for the coefficients β (top) and the Jacobian J (bottom) of the forward-mode differentiation of proximal coordinate descent (Algorithm 3) on multiple data sets.

Figure 11 :

 11 Figure 11: Local linear convergence of the Jacobian for sparse logistic regression. Distance to optimum for the coefficients β (top) and the Jacobian J (bottom) of the forward-mode differentiation of proximal coordinate descent (Algorithm 3) on multiple data sets.

Figure 12 :

 12 Figure 12: Weighted Lasso with hold-out criterion. Time to compute a single hypergradient as a function of the number of features on the gina data set. The regularization parameters have been chosen uniformly at random in the range [0, e λ max ].

Table 1 :

 1 Examples of non-smooth inner problems as in (1).

	Inner problem, Φ	f (β)	g j (β j , λ)	e λmax
	Lasso elastic net	1 2n y -Xβ 2 1 2n y -Xβ 2	e λ |β j |	1 n X y ∞

Table 2 :

 2 Examples of outer criteria used for hyperparameter selection.

	Criterion Hold-out mean squared error Stein unbiased risk estimate (SURE) 1 Hold-out logistic loss	Problem type Regression Regression Classification	Criterion C(β) n y val -X val β 2 1 y -Xβ 2 -nσ 2 + 2σ 2 dof(β) 1 n n i=1 ln(1 + e -y val i X val i: β )
	Hold-out smoothed Hinge loss 2	Classification	1 n	n i=1 (y val i , X val i: β)

  1 and 0 else.

		SMBO	1st order method
	0.17 0.21 0.26 0.32 0.39 0.48 0.60 0.73 0.90 1.11	λ 1 -λ max	λ 1 -λ max -10	0	Iterations

  γ 1 , . . . , γ p

	// jointly compute coef. & Jacobian
	for k = 1, . . . , n iter do
	for j = 1, . . . , p do
	// update the regression coefficients

Table 5 :

 5 Characteristics of the data sets used for the experiments.

			name	# samples n # features p # classes q	density
			breast cancer diabetes leukemia gina agnostic rcv1 real-sim news20 mnist usps	569 442 72 3,468 20,242 72,309 19,996 60.000 7,291	30 10 7,129 970 19,960 20,958 632,983 683 256	-------10 10	1 1 1 1 3.7 × 10 -3 2.4 × 10 -3 6.1 × 10 -4 2.2 × 10 -1 1
			rcv1 (multiclass) aloi	15,564 108,000	16,245 128	53 1,000	4.0 × 10 -3 2.4 × 10 -1
			Forward-mode PCD	Reverse-mode PCD	Cvxpylayers
	Time (s)	10 0 10 2
		10 -2	10 2		10 3	10 2	10 3
			# features p		# features p

e λ = e λ max / 10 e λ = e λ max / 100

Note that implicit refers to the implicit function theorem, but leads to an explicit formula for the gradient.

https://jermwatt.github.io/machine_learning_refined/notes/3_First_order_methods/3_9_ Normalized.html

Acknowledgments

This work was partially funded by the ERC Starting Grant SLAB ERC-StG-676943, the ANR BrAIN ANR-20-CHIA-0016, the ANR CaMeLOt ANR-20-CHIA-0001-01, and the ANR grant GraVa ANR-18-CE40-0005. Part of this work has been carried out at the Machine Learning Genoa (MaLGa) center, Università di Genova (IT). M. M. acknowledges the financial support of the European Research Council (grant SLING 819789).

Random

A. Implicit Differentiation Examples

For exposition purpose we provide instantiations of the implicit differentiation (Algorithm 5) for multiple optimization problems. For the Lasso (Algorithm 7)

the weighted Lasso (Algorithm 9)

and the dual of the SVM ŵ ∈ arg min

Algorithm 7 Lasso implicit differentiation input : λ ∈ R, > 0 // compute the solution of inner problem

Find β such that:

B. Additional Lemmas

B.1 Differentiability of the Proximal Operator

Here we recall results on the differentiability of the proximal operator at the optimum.

Lemma 20 [START_REF] Klopfenstein | Model identification and local linear convergence of coordinate descent[END_REF], Lemmas 2 and 3) Let 0 < γ j ≤ 1/L j . Let λ ∈ R r and Λ a neighborhood of λ. Consider a solution β ∈ arg min β∈R p Φ(β, λ) and Ŝ its generalized support. Suppose 1. Assumptions 2, 3 and 6 hold.

Assumption 4 hold on Λ.

We can now split the last summand in two parts and obtain the following bound, reminding that ρ(A) + δ < 1

Thus, (J (k) ) k∈N converges linearly towards its limit Ĵ .

C. Proof of Theorem 12

Theorem 12 Local linear convergence of the Jacobian. Let 0 < γ ≤ 1/L. Suppose Assumptions 2, 3 and 6 hold. Let λ ∈ R r , Λ be a neighborhood of λ, and Γ Λ β(λ) -γ∇f ( β(λ) ) : λ ∈ Λ . In addition, suppose hypotheses (H1) to (H4) from Theorem 9 are satisfied and the sequence (β (k) ) k∈N generated by Algorithm 1 (respectively by Algorithm 3) converges toward β.

Then, the sequence of Jacobians (J (k) ) k≥0 generated by the forward-mode differentiation of proximal gradient descent (Algorithm 1) (respectively by forward-mode differentiation of proximal coordinate descent, Algorithm 3) converges locally linearly towards Ĵ .

Proof We first prove Theorem 12 for proximal gradient descent.

Proximal gradient descent case. Solving Problem (1) with proximal gradient descent leads to the following updates:

) .

(

Consider the following sequence (J (k) ) k∈N defined by

Note that if prox γg(•,λ) is not differentiable with respect to the first variable at z (k) (respectively with respect to the second variable λ), any weak Jacobian can be used. When (H3) holds, differentiating Equation ( 23) with respect to λ yields exactly Equation ( 24).

where

A full epoch can then be written

The spectral radius of A is strictly bounded by 1 (Klopfenstein et al., 2020, Lemma 8): ρ(A) < 1. Assumptions 2 to 4 and 6 and the convergence of (β (k) ) k∈N toward β ensure local linear convergence of (β (k) ) k∈N (Klopfenstein et al., 2020, Thm. 2). Hence, we can write the update for the Jacobian after an update of the coordinates from 1 to p

with ( (k) ) k∈N converging linearly to 0. Recalling ρ(A) < 1, Theorem 21 and the last display yield our result using.

C.1 Proof of Theorem 13 (Approximate Hypergradients)

Theorem 13 Bound on the error of approximate hypergradient.

For λ ∈ R r , let β(λ) ∈ R p be the exact solution of the inner Problem (1), and Ŝ its generalized support. Suppose Assumptions 2, 3 and 6 hold. Let Λ be a neighborhood of λ, and Γ Λ β(λ) -γ∇f ( β(λ) ) : λ ∈ Λ . Suppose hypotheses (H1) to (H4) from Theorem 9 are satisfied. In addition suppose (H5) The application β → ∇ 2 f (β) is Lipschitz continuous.

(H6) The criterion β → ∇C(β) is Lipschitz continuous.

(H7) Both optimization problems in Algorithm 5 are solved up to precision with support identification: β (λ) -β(λ) ≤ , A is invertible, and

Then the error on the approximate hypergradient h returned by Algorithm 5 is of the order of magnitude of the error on β (λ) and v

Proof Overview of the proof. Our goal is to bound the error between the approximate hypergradient h returned by Algorithm 5 and the true hypergradient ∇L(λ). Following the analysis of [START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF], two sources of approximation errors arise when computing the hypergradient:

• One can control the product of the perturbations

Conditions are met to apply the result by Higham (2002, Thm 7.2), which leads to

The bound on ṽv finally yields a bound on the first quantity in Equation ( 3), vv

Bound on B -B 2 . We now bound the second quantity in Equation ( 3 (34)

Plugging Equations ( 33) and (34) into Equation (3) yields the desired result: ∇L(λ)-h = O( ).