# Implicit differentiation for fast hyperparameter selection in non-smooth
convex learning

Quentin Bertrand∗<EMAIL_ADDRESS>
Université Paris-Saclay, Inria, CEA, Palaiseau, France Quentin Klopfenstein∗
<EMAIL_ADDRESS>
Institut Mathématique de Bourgogne, Université de Bourgogne, Dijon, France
Mathurin Massias<EMAIL_ADDRESS>
MaLGa, DIBRIS, Università degli Studi di Genova, Genova, Italy Mathieu Blondel
<EMAIL_ADDRESS>
Google Research, Brain team, Paris, France Samuel Vaiter
<EMAIL_ADDRESS>
CNRS and Institut Mathématique de Bourgogne, Université de Bourgogne, Dijon,
France Alexandre Gramfort<EMAIL_ADDRESS>
Université Paris-Saclay, Inria, CEA, Palaiseau, France Joseph Salmon
<EMAIL_ADDRESS>
IMAG, Université de Montpellier, CNRS, Montpellier, France

###### Abstract

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.

##  1 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 the regularization hyperparameter controls the
trade-off between a data fidelity term and a regularization term. Such
estimators, including Ridge regression (Hoerl and Kennard, 1970), Lasso
(Tibshirani, 1996; Chen et al., 1998), elastic net (Zou and Hastie, 2005),
sparse logistic regression (Koh et al., 2007), support-vector machine/SVM
(Boser et al., 1992; Platt, 1999) are often cast as an optimization problem
(Table 1)

Table 1: Examples of non-smooth inner problems as in (1). Inner problem, $\Phi$ | $f(\beta)$ | $g_{j}(\beta_{j},\lambda)$ | $e^{\lambda_{\max}}$
---|---|---|---
Lasso | $\frac{1}{2n}\|y-X\beta\|^{2}$ | $e^{\lambda}|\beta_{j}|$ | $\tfrac{1}{n}\lVert X^{\top}y\rVert_{\infty}$
elastic net | $\frac{1}{2n}\|y-X\beta\|^{2}$ | $e^{\lambda_{1}}|\beta_{j}|+\tfrac{1}{2}e^{\lambda_{2}}\beta_{j}^{2}$ | $\tfrac{1}{n}\lVert X^{\top}y\rVert_{\infty}$
sparse log. reg. | $\frac{1}{n}\sum_{i=1}^{n}\ln(1+e^{-y_{i}X_{i:}\beta})$ | $e^{\lambda}|\beta_{j}|$ | $\tfrac{1}{2n}\lVert X^{\top}y\rVert_{\infty}$
dual SVM | $\frac{1}{2}\lVert(y\odot X)^{\top}\beta\rVert^{2}-\sum_{j=1}^{p}\beta_{j}$ | $\iota_{[0,e^{\lambda}]}(\beta_{j})$ | $-$

$\displaystyle\hat{\beta}^{(\lambda)}\in\operatorname*{\mathrm{arg\,min}}_{\beta\in\mathbb{R}^{p}}\Phi(\beta,\lambda)\triangleq
f(\beta)+\underbrace{\sum_{j=1}^{p}g_{j}(\beta_{j},\lambda)}_{\triangleq
g(\beta,\lambda)}\enspace,$ (1)

with smooth $f:\mathbb{R}^{p}\rightarrow\mathbb{R}$ (i.e., with Lipschitz
gradient), proper closed convex (possibly non-smooth) functions
$g_{j}(\cdot,\lambda)$, and a regularization hyperparameter
$\lambda\in\mathbb{R}^{r}$. In the examples of Table 1, the computation of $f$
involves a design matrix $X\in\mathbb{R}^{n\times p}$; and the cost of
computing $\nabla f(\beta)$ is $\mathcal{O}(np)$. In the SVM example, since we
consider the dual problem, we chose to reverse the roles of $n$ and $p$ to
enforce $\beta\in\mathbb{R}^{p}$. We often drop the $\lambda$ dependency and
write $\hat{\beta}$ instead of $\hat{\beta}^{(\lambda)}$ when it is clear from
context.

For a fixed $\lambda$, the issue of solving efficiently Equation 1 has been
largely explored. If the functions $g_{j}$ are smooth, one can use solvers
such as L-BFGS (Liu and Nocedal, 1989), SVRG (Johnson and Zhang, 2013; Zhang
et al., 2013), or SAGA (Defazio et al., 2014). When the functions $g_{j}$ are
non-smooth, Equation 1 can be tackled efficiently with stochastic algorithms
(Pedregosa et al., 2017) or using working set methods (Fan and Lv, 2008;
Tibshirani et al., 2012) combined with coordinate descent (Tseng and Yun,
2009), see overview by Massias et al. (2020). The question of _model
selection_ , i.e., how to select the hyperparameter $\lambda\in\mathbb{R}^{r}$
(potentially multidimensional), is more open, especially when the dimension
$r$ of the regularization hyperparameter $\lambda$ 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 $\lambda$ (Lounici,
2008; Bickel et al., 2009; Belloni et al., 2011). 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 $\lambda$.

A popular approach for hyperparameter selection is _hyperparameter
optimization_ (Kohavi and John, 1995; Hutter et al., 2015; Feurer and Hutter,
2019): one selects the hyperparameter $\lambda$ such that the regression
coefficients $\hat{\beta}^{(\lambda)}$ minimize a given criterion
$\mathcal{C}:\mathbb{R}^{p}\rightarrow\mathbb{R}$. Here $\mathcal{C}$ should
ensure good generalization, or avoid overcomplex models. Common examples (see
Table 2) include the hold-out loss (Devroye and Wagner, 1979), the cross-
validation loss (CV, Stone and Ramer 1965, see Arlot and Celisse 2010 for a
survey), the AIC (Akaike, 1974), BIC (Schwarz, 1978) or SURE (Stein, 1981)
criteria. Formally, the hyperparameter optimization problem is a bilevel
optimization problem (Colson et al., 2007):

$\displaystyle\begin{aligned}
&\operatorname*{\mathrm{arg\,min}}_{\lambda\in\mathbb{R}^{r}}\left\\{\mathcal{L}(\lambda)\triangleq\mathcal{C}\left(\hat{\beta}^{(\lambda)}\right)\right\\}\\\
&{s.t.\leavevmode\nobreak\
}\hat{\beta}^{(\lambda)}\in\operatorname*{\mathrm{arg\,min}}_{\beta\in\mathbb{R}^{p}}\Phi(\beta,\lambda)\enspace.\end{aligned}$
(2)

Table 2: Examples of outer criteria used for hyperparameter selection.

Criterion | Problem type | Criterion $\mathcal{C}(\beta)$
---|---|---
Hold-out mean squared error | Regression | $\frac{1}{n}\|y^{\text{val}}-X^{\text{val}}\beta\|^{2}$
Stein unbiased risk estimate (SURE)111 For a linear model $y=X\beta+\varepsilon$, with $\varepsilon\sim~\mathcal{N}(0,\sigma^{2})$, the degree of freedom (dof, Efron 1986) is defined as $\text{dof}(\beta)=\sum_{i=1}^{n}\text{cov}(y_{i},(X\beta)_{i})/\sigma^{2}$. | Regression | $\|y-X\beta\|^{2}-n\sigma^{2}+2\sigma^{2}\text{dof}(\beta)$
Hold-out logistic loss | Classification | $\frac{1}{n}\sum_{i=1}^{n}\ln(1+e^{-y^{\text{val}}_{i}X_{i:}^{\text{val}}\beta})$
Hold-out smoothed Hinge loss222The smoothed Hinge loss is given by $\ell(x)=\frac{1}{2}-x$ if $x\leq 0$, $\frac{1}{2}(1-x)^{2}$ if $0\leq x\leq 1$ and 0 else. | Classification | $\frac{1}{n}\sum_{i=1}^{n}\ell(y^{\text{val}}_{i},X_{i:}^{\text{val}}\beta)$

Popular approaches to solve (the generally non-convex) Equation 2 include
zero-order optimization (gradient-free) techniques such as grid-search,
random-search (Rastrigin, 1963; Bergstra and Bengio, 2012; Bergstra et al.,
2013) or Sequential Model-Based Global Optimization (SMBO), often referred to
as Bayesian optimization (Mockus, 1989; Jones et al., 1998; Forrester et al.,
2008; Brochu et al., 2010; Snoek et al., 2012). Grid-search is a naive
discretization of Equation 2. It consists in evaluating the outer function
$\mathcal{L}$ on a grid of hyperparameters, solving one inner optimization
Equation 1 for each $\lambda$ in the grid (see Figure 1). For each inner
problem solution $\hat{\beta}^{(\lambda)}$, the criterion
$\mathcal{C}(\hat{\beta}^{(\lambda)})$ 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 $\mathcal{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 Bergstra et al. (2011). 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
$\lambda\mapsto\hat{\beta}^{(\lambda)}$ is well-defined and almost everywhere
differentiable, first-order optimization methods are well suited to solve the
bilevel optimization Equation 2. Using the chain rule, the gradient of
$\mathcal{L}$ with respect to $\lambda$, also referred to as the
_hypergradient_ , evaluates to

$\displaystyle\nabla_{\lambda}\mathcal{L}(\lambda)$
$\displaystyle=\hat{\mathcal{J}}^{\top}_{(\lambda)}\nabla\mathcal{C}(\hat{\beta}^{(\lambda)})\enspace,$
(3)

with $\hat{\mathcal{J}}_{(\lambda)}\in\mathbb{R}^{p\times r}$ the _Jacobian_
of the function $\lambda\mapsto\hat{\beta}^{(\lambda)}$,

$\displaystyle\hat{\mathcal{J}}_{(\lambda)}\triangleq\begin{pmatrix}&\tfrac{\partial\hat{\beta}^{(\lambda)}_{1}}{\partial\lambda_{1}}&\ldots&\tfrac{\partial\hat{\beta}^{(\lambda)}_{1}}{\partial\lambda_{r}}\\\
&\vdots&\ldots&\vdots\\\
&\tfrac{\partial\hat{\beta}^{(\lambda)}_{p}}{\partial\lambda_{1}}&\ldots&\tfrac{\partial\hat{\beta}^{(\lambda)}_{p}}{\partial\lambda_{r}}\end{pmatrix}\enspace.$
(4)

An important challenge of applying first-order methods to solve Equation 2 is
evaluating the hypergradient in Equation 3. There are three main algorithms to
compute the hypergradient $\nabla_{\lambda}\mathcal{L}(\lambda)$: implicit
differentiation (Larsen et al., 1996; Bengio, 2000) and automatic
differentiation using the reverse-mode (Linnainmaa, 1970; LeCun et al., 1998)
or the forward-mode (Wengert, 1964; Deledalle et al., 2014; Franceschi et al.,
2017). As illustrated in Figure 1, once the hypergradient in Equation 3 has
been computed, one can solve Equation 2 with first-order schemes, e.g.,
gradient descent.

Figure 1: 5-fold cross-validation error $\mathcal{C}(\beta^{(\lambda)})$:
(top) Lasso CV error with respect to $\lambda$ for multiple hyperparameter
optimization methods on the _real-sim_ dataset, and (bottom) elastic net CV
error with respect to $\lambda_{1}$ and $\lambda_{2}$ on the _rcv1_ dataset.
Crosses represent the $10$ (top) or $25$ (bottom) first error evaluations for
each method.

#### Contributions.

We are interested in tackling the bilevel optimization Equation 2, with a non-
smooth inner optimization Equation 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 Equation
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 state-of-the-art solvers (Algorithm 5). We combine in
Algorithm 6 this fast hypergradient computation with a gradient descent scheme
to solve Equation 2.

  * •

We provide extensive experiments on diverse datasets 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 Equation 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 $\lVert\cdot\rVert$ the Euclidean norm on vectors. For a set $S$, we
denote by $S^{c}$ its complement. We denote $[p]=\\{1,\dots,p\\}$. We denote
by $(e_{j})_{j=1}^{p}$ the vectors of the canonical basis of $\mathbb{R}^{p}$.
We denote the coordinate-wise multiplication of two vectors $u$ and $v$ by
$u\odot v$, and by $u\odot 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\in\mathbb{R}^{n\times n}$ is
denoted $\rho(M)=\max_{i}|s_{i}|$ where $s_{1},\ldots,s_{n}$ are the
eigenvalues of $M$. For a matrix $M$, we write that $M\succ 0$ if $M$ is
positive definite. The regularization parameter, possibly multivariate, is
denoted by $\lambda=(\lambda_{1},\dots,\lambda_{r})^{\top}\in\mathbb{R}^{r}$.
We denote
$\hat{\mathcal{J}}_{(\lambda)}\triangleq(\nabla_{\lambda}\hat{\beta}_{1}^{(\lambda)},\dots,\nabla_{\lambda}\hat{\beta}_{p}^{(\lambda)})^{\top}\in\mathbb{R}^{p\times
r}$ the weak Jacobian (Evans and Gariepy, 1992) of $\hat{\beta}^{(\lambda)}$
with respect to $\lambda$.

#### Convex analysis.

For a convex function $h:\mathbb{R}^{p}\to\mathbb{R}$, the proximal operator
of $h$ is defined, for any $x\in\mathbb{R}^{p}$, as:
$\operatorname{prox}_{h}(x)=\operatorname*{\mathrm{arg\,min}}_{y\in\mathbb{R}^{p}}\frac{1}{2}\|x-y\|^{2}+h(y)$.
The subdifferential of $h$ at $x$ is denoted $\partial
h(x)=\left\\{u\in\mathbb{R}^{p}\,:\,\forall z\in\mathbb{R}^{p},h(z)\geq
h(x)+u^{\top}(z-x)\right\\}$. A function is said to be _smooth_ if it has
Lipschitz gradients. Let $f$ be a $L$-smooth function. Lipschitz constants of
the functions $\nabla_{j}f$ are denoted by $L_{j}$; hence for all
$x\in\mathbb{R}^{p}$, $h\in\mathbb{R}$:

$\displaystyle|\nabla_{j}f(x+he_{j})-\nabla_{j}f(x)|\leq L_{j}|h|\enspace.$

For a function $f$, its gradient restricted to the indices in a set $S$ is
denoted $\nabla_{S}f$. For a set $\Xi\subset\mathbb{R}^{p}$, its relative
interior is noted $\operatorname{ri}(\Xi)$, and its indicator function is
defined for any $x\in\mathbb{R}^{p}$ by $\iota_{\Xi}(x)=0$ if $x\in\Xi$ and
$+\infty$ otherwise. A function
$h:\mathbb{R}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ is said to be proper if
$\operatorname{dom}(h)=\\{x\in\mathbb{R}:h(x)<+\infty\\}\neq\emptyset$), and
closed if for any $\alpha\in\mathbb{R}$, the sublevel set
$\\{x\in\operatorname{dom}(h):h(x)\leq\alpha\\}$ is a closed set.

For a function $\psi:\mathbb{R}^{p}\times\mathbb{R}^{r}\mapsto\mathbb{R}^{p}$,
we denote $\partial_{z}\psi$ the weak Jacobian with respect to the first
variable and $\partial_{\lambda}\psi$ the weak Jacobian with respect to the
second variable. The proximal operator of $g(\cdot,\lambda)$ can be seen as
such a function $\psi$ of $\beta$ and $\lambda$ (see Table 1 for examples):

$\displaystyle\mathbb{R}^{p}\times\mathbb{R}^{r}$
$\displaystyle\rightarrow\mathbb{R}^{p}$ $\displaystyle(z,\lambda)$
$\displaystyle\mapsto\operatorname{prox}_{g(\cdot,\lambda)}(z)=\psi(z,\lambda)\enspace.$

In this case we denote
$\partial_{z}\operatorname{prox}_{g(\cdot,\lambda)}\triangleq\partial_{z}\psi$
and
$\partial_{\lambda}\operatorname{prox}_{g(\cdot,\lambda)}\triangleq\partial_{\lambda}\psi$.
Since we consider only separable penalties $g(\cdot,\lambda)$,
$\partial_{z}\operatorname{prox}_{g(\cdot,\lambda)}$ is a diagonal matrix, so
to make notation lighter, we write
$\partial_{z}\operatorname{prox}_{g(\cdot,\lambda)}$ for its diagonal. We thus
have

$\displaystyle\partial_{z}\operatorname{prox}_{g(\cdot,\lambda)}=(\partial_{z}\operatorname{prox}_{g_{j}(\cdot,\lambda)})_{j\in[p]}$
$\displaystyle\in\mathbb{R}^{p}\quad(\text{by separability of $g$})$
$\displaystyle\partial_{\lambda}\operatorname{prox}_{g(\cdot,\lambda)}$
$\displaystyle\in\mathbb{R}^{p\times r}\enspace.$

Explicit partial derivatives formulas for usual proximal operators can be
found in Table 3.

Table 3: Partial derivatives of proximal operators used.

$g_{j}(\beta_{j},\lambda)$ | $\operatorname{prox}_{g_{j}(\cdot,\lambda)}(z_{j})$ | $\partial_{z}\operatorname{prox}_{g_{j}(\cdot,\lambda)}(z_{j})$ | $\partial_{\lambda}\operatorname{prox}_{g_{j}(\cdot,\lambda)}(z_{j})$
---|---|---|---
$e^{\lambda}\beta_{j}^{2}/2$ | $z_{j}/(1+e^{\lambda})$ | $1/(1+e^{\lambda})$ | $-z_{j}e^{\lambda}/(1+e^{\lambda})^{2}$
$e^{\lambda}|\beta_{j}|$ | $\operatorname{ST}(z_{j},e^{\lambda})$ | $|\operatorname{sign}(\operatorname{ST}(z_{j},e^{\lambda}))|$ | $-e^{\lambda}\operatorname{sign}(\operatorname{ST}(z_{j},e^{\lambda}))$
$e^{\lambda_{1}}|\beta_{j}|+\tfrac{1}{2}e^{\lambda_{2}}\beta_{j}^{2}$ | $\frac{\operatorname{ST}(z_{j},e^{\lambda_{1}})}{1+e^{\lambda_{2}}}$ | $\frac{|\operatorname{sign}(\operatorname{ST}(z_{j},e^{\lambda_{1}}))|}{1+e^{\lambda_{2}}}$ | $\left(\frac{-e^{\lambda_{1}}\operatorname{sign}(\operatorname{ST}(z_{j},e^{\lambda_{1}}))}{1+e^{\lambda_{2}}},\frac{-\operatorname{ST}(z_{j},e^{\lambda_{1}})e^{\lambda_{2}}}{(1+e^{\lambda_{2}})^{2}}\right)$
$\iota_{[0,e^{\lambda}]}(\beta_{j})$ | $\max(0,\min(z_{j},e^{\lambda}))$ | $\mathds{1}_{]0,e^{\lambda}[}(z_{j})$ | $e^{\lambda}\mathds{1}_{z_{j}>e^{\lambda}}$

##  2 Related work

The main challenge to evaluate the hypergradient
$\nabla_{\lambda}\mathcal{L}(\lambda)$ is the computation of the Jacobian
$\mathcal{J}_{(\lambda)}$. We first focus on the case where
$\Phi(\cdot,\lambda)$ is convex and smooth for any $\lambda$.

#### Implicit differentiation.

We recall how the _implicit differentiation_ 333Note that _implicit_ refers to
the implicit function theorem, but leads to an _explicit_ formula for the
gradient. formula of the gradient $\nabla_{\lambda}\mathcal{L}(\lambda)$ is
obtained for smooth inner optimization problems. We will provide a
generalization to non-smooth optimization problems in Section 3.2.

######  Theorem 1 (Bengio 2000)

Let
$\hat{\beta}^{(\lambda)}\in\operatorname*{\mathrm{arg\,min}}_{\beta\in\mathbb{R}^{p}}\Phi(\beta,\lambda)$
be a solution of Equation 1. Assume that for all $\lambda>0$,
$\Phi(\cdot,\lambda)$ is a convex smooth function,
$\nabla_{\beta}^{2}\Phi(\hat{\beta}^{(\lambda)},\lambda)\succ 0$, and that for
all $\beta\in\mathbb{R}^{p}$, $\Phi(\beta,\cdot)$ is differentiable over
$]0,+\infty[$. Then the hypergradient $\nabla_{\lambda}\mathcal{L}(\lambda)$
reads:

$\underbrace{\nabla_{\lambda}\mathcal{L}(\lambda)}_{\in\mathbb{R}^{r}}=\underbrace{-\nabla_{\beta,\lambda}^{2}\Phi(\hat{\beta}^{(\lambda)},\lambda)}_{\in\mathbb{R}^{r\times
p}}{\underbrace{\left(\nabla_{\beta}^{2}\Phi(\hat{\beta}^{(\lambda)},\lambda)\right)}_{\in\mathbb{R}^{p\times
p}}}^{-1}\underbrace{\nabla\mathcal{C}(\hat{\beta}^{(\lambda)})}_{\in\mathbb{R}^{p}}\enspace.$
(5)

Proof  For a smooth convex function $\beta\mapsto\Phi(\beta,\lambda)$ the
first-order condition writes:

$\displaystyle\nabla_{\beta}\Phi(\hat{\beta}^{(\lambda)},\lambda)=0\enspace,$
(6)

for any $\hat{\beta}^{(\lambda)}$ solution of the inner problem. Moreover, if
$\lambda\mapsto\nabla_{\beta}\Phi(\hat{\beta}^{(\lambda)},\lambda)$ is
differentiable, differentiating Equation 6 with respect to $\lambda$ leads to:

$\displaystyle\nabla_{\beta,\lambda}^{2}\Phi(\hat{\beta}^{(\lambda)},\lambda)+\hat{\mathcal{J}}^{\top}_{(\lambda)}\nabla_{\beta}^{2}\Phi(\hat{\beta}^{(\lambda)},\lambda)=0\enspace.$
(7)

The Jacobian $\hat{\mathcal{J}}^{\top}_{(\lambda)}$ is computed by solving the
following linear system:

$\displaystyle\hat{\mathcal{J}}_{(\lambda)}^{\top}=-\underbrace{\nabla_{\beta,\lambda}^{2}\Phi(\hat{\beta}^{(\lambda)},\lambda)}_{\in\mathbb{R}^{r\times
p}}{\underbrace{\left(\nabla_{\beta}^{2}\Phi(\hat{\beta}^{(\lambda)},\lambda)\right)}_{\in\mathbb{R}^{p\times
p}}}^{-1}\enspace.$ (8)

Plugging Equation 8 into Equation 3 yields the desired result.

The computation of the gradient via implicit differentiation (Equation 5)
involves the resolution of a $p\times p$ linear system (Bengio, 2000, Sec. 4).
This potentially large linear system can be solved using different algorithms
such as conjugate gradient (Hestenes and Stiefel 1952, as in Pedregosa 2016)
or fixed point methods (Lions and Mercier 1979; Tseng and 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 (Chapelle et al., 2002; Seeger, 2008), weighted Ridge estimator (Foo et
al., 2008), neural networks (Lorraine et al., 2019) or meta-learning
(Franceschi et al., 2018; Rajeswaran et al., 2019). In addition to
hyperparameter selection, it has been applied successfully in natural language
processing (Bai et al., 2019) and computer vision (Bai et al., 2020).

Equation 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 6 is not exactly satisfied at machine precision, and
consequently the linear system to solve Equation 5 does not lead to the exact
gradient $\nabla_{\lambda}\mathcal{L}(\lambda)$, see Ablin et al. (2020) for
quantitative convergence results. However, Pedregosa (2016) 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 $\hat{\beta}$. Interestingly, this approximation scheme was shown to yield
significant practical speedups when solving Equation 2, while preserving
theoretical properties of convergence toward the optimum.

#### Iterative differentiation.

Iterative differentiation computes the gradient
$\nabla_{\lambda}\mathcal{L}(\lambda)$ by differentiating through the iterates
of the algorithm used to solve Equation 1. Iterative differentiation can be
applied using the forward-mode (Wengert 1964; Deledalle et al. 2014;
Franceschi et al. 2017) or the reverse-mode (Linnainmaa 1970; LeCun et al.
1998; Domke 2012). 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,
Lions and Mercier 1979; Combettes and Wajs 2005), 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 Domke
(2012), 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 Equation 1 using an iterative solver, but requires storing
the iterates along the computation for a backward evaluation of the
hypergradient (Algorithm 2). Alternatively, the forward-mode computes jointly
the solution along with the gradient $\nabla_{\lambda}\mathcal{L}(\lambda)$.
It is memory efficient (no iterates storage) but more computationally
expensive when the number of hyperparameters ($r$) is large; see Baydin et al.
(2018) for a survey.

#### Resolution of the bilevel Equation 2.

From a theoretical point of view, solving Equation 2 using gradient-based
methods is also challenging, and results in the literature are quite scarce.
Kunisch and Pock (2013) studied the convergence of a semi-Newton algorithm
where both the outer and inner problems are smooth. Franceschi et al. (2018)
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 $\lambda>0$ but do not have
second-order assumptions on $\Phi$. Recent results (Ghadimi and Wang, 2018; Ji
et al., 2020; Mehmood and Ochs, 2021) have provided quantitative convergence
toward a global solution of Equation 2, but under global joint convexity
assumption and exact knowledge of the gradient Lipschitz constant.

input : $\lambda\in\mathbb{R}^{r},\gamma>0,n_{\mathrm{iter}}\in\mathbb{N}$,
$\beta^{(0)}\in\mathbb{R}^{p}$, $\mathcal{J}^{(0)}\in\mathbb{R}^{p\times r}$

// jointly compute coef. & Jacobian

for  _$k=1,\dots,n_{\mathrm{iter}}$_ do

// update the regression coefficients

$z^{(k)}=\beta^{(k-1)}-\gamma\nabla f(\beta^{(k-1)})$ ;

// GD step

$\mathrm{d}z^{(k)}=\mathcal{J}^{(k-1)}-\gamma\nabla^{2}f(\beta^{(k-1)})\mathcal{J}^{(k-1)}$

$\beta^{(k)}=\operatorname{prox}_{\gamma g(\cdot,\lambda)}(z^{(k)})$ ;

// prox. step

// update the Jacobian

$\mathcal{J}^{(k)}=\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z^{(k)})\odot\mathrm{d}z^{(k)}$

$\mathcal{J}^{(k)}\mathrel{+}=\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z^{(k)})$ ;

// $\mathcal{O}(pr)$

$v=\nabla\mathcal{C}(\beta^{n_{\mathrm{iter}}})$

return _$\beta^{n_{\mathrm{iter}}},\mathcal{J}^{n_{\mathrm{iter}}\top}v$_

Algorithm 1 Forward-mode PGD

input : $\lambda\in\mathbb{R}^{r},\gamma>0,n_{\mathrm{iter}}\in\mathbb{N}$,
$\beta^{(0)}\in\mathbb{R}^{p}$

// computation of $\hat{\beta}$

for  _$k=1,\dots,n_{\mathrm{iter}}$_ do

$z^{(k)}=\beta^{(k-1)}-\gamma\nabla f(\beta^{(k-1)})$ ;

// GD step

$\beta^{(k)}=\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(z^{(k)}\right)$ ;

// prox. step

// backward computation of the gradient $g$

$v=\nabla\mathcal{C}(\beta^{(n_{\mathrm{iter}})})$, $h=0_{\mathbb{R}^{r}}$

for  _$k=n_{\mathrm{iter}},n_{\mathrm{iter}}-1,\dots,1$_ do

$h\mathrel{+}=v^{\top}\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z^{(k)})$ ;

// $\mathcal{O}(pr)$

$v\leftarrow\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z^{(k)})\odot v$ ;

// $\mathcal{O}(p)$

$v\leftarrow(\operatorname{Id}-\gamma\nabla^{2}f(\beta^{(k)}))v$ ;

// $\mathcal{O}(np)$

return _ $\beta^{n_{\mathrm{iter}}},h$ _

Algorithm 2 Reverse-mode PGD

##  3 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 Equation
2 with non-smooth inner optimization Equation 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 Equation 2.

###  3.1 Theoretical framework

Differentiability of the regularization path. Before applying first-order
methods to tackle Equation 2, one must ensure that the regularization path
$\lambda\mapsto\hat{\beta}^{(\lambda)}$ is almost everywhere differentiable
(as in Figure 2). This is the case for the Lasso (Mairal and Yu, 2012) and the
SVM (Hastie et al., 2004; Rosset and Zhu, 2007) since solution paths are
piecewise differentiable (see Figure 2). Results for nonquadratic datafitting
terms are scarcer: Friedman et al. (2010) 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
Equation 1, we believe it is an open question and leave it for future work.

Figure 2:  Regularization paths (coefficient values as a function of
$\lambda$), on the _diabetes_ and _breast cancer_ datasets 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.

#### 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_ 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 (Liang et al., 2014; Vaiter et al., 2018) and proximal
coordinate descent, see Nutini (2018, Sec. 6.2) or Klopfenstein et al. (2020).
Once these constraints have been identified convergence is linear towards the
Jacobian (see Theorems 12, 10, 11 and 3).

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

######  Assumption 2 (Smoothness)

The function $f:\mathbb{R}^{p}\rightarrow\mathbb{R}$ is a convex,
differentiable function, with a $L$-Lipschitz gradient.

######  Assumption 3 (Proper, closed, convex)

For all $\lambda\in\mathbb{R}^{r}$, for any $j\in[p]$, the function
$g_{j}(\cdot,\lambda):\mathbb{R}\rightarrow\mathbb{R}$ is proper, closed and
convex.

######  Assumption 4 (Non-degeneracy)

The problem admits at least one solution:

$\operatorname*{\mathrm{arg\,min}}_{\beta\in\mathbb{R}^{p}}\Phi(\beta,\lambda)\neq\emptyset\enspace,$

and, for any $\hat{\beta}$ solution of Equation 1, we have

$-\nabla
f(\hat{\beta})\in\operatorname{ri}\left(\partial_{\beta}g(\hat{\beta},\lambda)\right)\enspace.$

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
$\hat{\beta}\in\operatorname*{\mathrm{arg\,min}}_{\beta\in\mathbb{R}^{p}}\Phi(\beta,\lambda)$,
its _generalized support_ $\hat{S}\subseteq[p]$ is the set of indices
$j\in[p]$ such that $g_{j}$ is differentiable at $\hat{\beta}_{j}$:

$\displaystyle\hat{S}\triangleq\\{j\in[p]:\partial_{\beta}g_{j}(\hat{\beta}_{j},\lambda)\text{
is a singleton}\\}\enspace.$

An iterative algorithm is said to achieve finite support identification if its
iterates $\beta^{(k)}$ converge to $\hat{\beta}$, and there exists $K\geq 0$
such that for all $j\notin\hat{S}$, for all $k\geq
K,\beta_{j}^{(k)}=\hat{\beta}_{j}$.

Examples. For the $\ell_{1}$ norm (promoting sparsity),
$g_{j}(\hat{\beta}_{j},\lambda)=e^{\lambda}|\hat{\beta}_{j}|$, the generalized
support is $\hat{S}\triangleq\\{j\in[p]:\hat{\beta}_{j}\neq 0\\}$. This set
corresponds to the indices of the non-zero coefficients, which is the usual
support definition. For the SVM estimator,
$g_{j}(\hat{\beta}_{j},\lambda)=\iota_{[0,e^{\lambda}]}(\hat{\beta_{j}})$.
This function is non-differentiable at $0$ and at $e^{\lambda}$. The
generalized support for the SVM estimator then corresponds to the set of
indices such that $\hat{\beta}_{j}\in]0,e^{\lambda}[$.

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

######  Assumption 6 (Locally $\mathcal{C}^{2}$ and $\mathcal{C}^{3}$)

The map $\beta\mapsto f(\beta)$ is locally $\mathcal{C}^{3}$ around
$\hat{\beta}$. For all $\lambda\in\mathbb{R}^{r}$, for all $j\in\hat{S}$ the
map $g_{j}(\cdot,\lambda)$ is locally $\mathcal{C}^{2}$ around
$\hat{\beta}_{j}$.

######  Assumption 7 (Restricted injectivity)

Let $\hat{\beta}$ be a solution of Equation 1 and $\hat{S}$ its generalized
support. The solution $\hat{\beta}$ satisfies the following restricted
injectivity condition:

$\displaystyle\nabla^{2}_{\hat{S},\hat{S}}f(\hat{\beta})\succ 0\enspace.$

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_. Assumptions 6 and 7 are classical for the
analysis (Liang et al., 2017) and sufficient to derive rates of convergence
for the Jacobian of the inner problem once the generalized support has been
identified.

The next lemma guarantees uniqueness of Equation 1 under Assumptions 4 and 7.

######  Lemma 8 (Liang et al. 2017, Prop. 4.1)

Assume that there exists a neighborhood $\Lambda$ of $\lambda$ such that
Assumptions 4 and 7 are satisfied for every $\lambda\in\Lambda$. Then for
every $\lambda\in\Lambda$, Equation 1 has a unique solution, and the map
$\lambda\mapsto\hat{\beta}^{(\lambda)}$ is well-defined on $\Lambda$.

We first show how implicit and iterative differentiation can be used with a
non-smooth inner problem. Peyré and Fadili (2011) proposed to smooth the inner
optimization problem, Ochs et al. (2015); Frecon et al. (2018) 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 (Gould et al., 2016; Amos and
Kolter, 2017; Agrawal et al., 2019), Lasso-type problems (Mairal et al., 2012;
Bertrand et al., 2020), total variation penalties (Cherkaoui et al., 2020) and
generalized to strongly monotone operators (Winston and Kolter, 2020).

###  3.2 Hypergradient computation: implicit differentiation

The exact proof of Theorem 1 cannot be applied when
$\beta\mapsto\Phi(\beta,\lambda)$ is non-smooth, as Equations 7 and 6 no
longer hold. Nevertheless, instead of the optimality condition of smooth
optimization, Equation 6, one can leverage the fixed point iteration of
proximal gradient descent, which we will see in Equation 11. The main
theoretical challenge is to show the differentiability of the function
$\beta\mapsto\operatorname{prox}_{\gamma g}(\beta-\gamma\nabla f(\beta))$.
Besides, taking advantage of the generalized sparsity of the regression
coefficients $\hat{\beta}^{(\lambda)}$, one can show that the Jacobian
$\hat{\mathcal{J}}$ is row-sparse, leading to substantial computational
benefits when computing the hypergradient
$\nabla_{\lambda}\mathcal{L}(\lambda)$) for Equation 1,

######  Theorem 9 (Non-smooth implicit formula)

Suppose Assumptions 2, 3 and 6 hold. Let $0<\gamma\leq 1/L$, where $L$ is the
Lipschitz constant of $\nabla f$. Let $\lambda\in\mathbb{R}^{r}$, $\Lambda$ be
a neighborhood of $\lambda$, and
$\Gamma^{\Lambda}\triangleq\left\\{\hat{\beta}^{(\lambda)}-\gamma\nabla
f(\hat{\beta}^{(\lambda)})\,:\,\lambda\in\Lambda\right\\}$. In addition,

  1. (H1)

Suppose Assumptions 4 and 7 hold on $\Lambda$.

  2. (H2)

Suppose $\lambda\mapsto\hat{\beta}^{(\lambda)}$ is continuously differentiable
on $\Lambda$.

  3. (H3)

Suppose for all $z\in\Gamma^{\Lambda}$,
$\lambda\mapsto\operatorname{prox}_{\gamma g(\cdot,\lambda)}(z)$ is
continuously differentiable on $\Lambda$.

  4. (H4)

Suppose $\partial_{z}\operatorname{prox}_{\gamma g(\cdot,\lambda)}$ and
$\partial_{\lambda}\operatorname{prox}_{\gamma g(\cdot,\lambda)}$ are
Lipschitz continuous on $\Gamma^{\Lambda}\times\Lambda$.

Let $\hat{\beta}\triangleq\hat{\beta}^{(\lambda)}$ be the solution of Equation
1, $\hat{S}$ its generalized support of cardinality $\hat{s}$. Then the
Jacobian $\hat{\mathcal{J}}$ of the inner Equation 1 is given by the following
formula,

$\hat{z}=\hat{\beta}-\gamma\nabla f(\hat{\beta})$, and
$A\triangleq\operatorname{Id}_{\hat{s}}-\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})_{\hat{S}}\odot\left(\operatorname{Id}_{\hat{s}}-\gamma\nabla^{2}_{\hat{S},\hat{S}}f(\hat{\beta})\right)$:

$\displaystyle\hat{\mathcal{J}}_{\hat{S}^{c}:}$
$\displaystyle=\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(\hat{z}\right)_{\hat{S}^{c}}\enspace,$ (9)
$\displaystyle\hat{\mathcal{J}}_{\hat{S}:}$
$\displaystyle=A^{-1}\left(\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})_{\hat{S}}-\gamma\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})_{\hat{S}}\odot\nabla^{2}_{\hat{S},\hat{S}^{c}}f(\hat{\beta})\hat{\mathcal{J}}_{\hat{S}^{c}}\right)\enspace.$
(10)

Proof  According to Lemma 8, Assumptions 7 and 4 ensure Equation 1 has a
unique minimizer and $\lambda\mapsto\hat{\beta}^{(\lambda)}$ is well-defined
on $\Lambda$. We consider the proximal gradient descent fixed point equation:

$\displaystyle\hat{\beta}^{(\lambda)}=\operatorname{prox}_{\gamma
g_{(\cdot,\lambda)}}\left(\hat{\beta}^{(\lambda)}-\gamma\nabla
f(\hat{\beta}^{(\lambda)})\right)\enspace.$ (11)

Together with the conclusion of Lemma 8, Assumptions 2 and 6, and given (H2),
(H3) and (H4), we have that
$\lambda\mapsto\psi\left(\beta^{(\lambda)}-\gamma\nabla
f(\hat{\beta}^{(\lambda)}),\lambda\right)\triangleq\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(\hat{\beta}^{(\lambda)}-\gamma\nabla
f(\hat{\beta}^{(\lambda)})\right)$ is differentiable at $\lambda$. One can
thus differentiate Equation 11 with respect to $\lambda$, which leads to:

$\displaystyle\hat{\mathcal{J}}=\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})\odot\left(\operatorname{Id}-\gamma\nabla^{2}f(\hat{\beta})\right)\hat{\mathcal{J}}+\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(\hat{z}\right)\enspace,$ (12)

with $\hat{z}=\hat{\beta}-\gamma\nabla f(\hat{\beta})$. In addition to
$0<\gamma<1/L\leq 1/L_{j}$, the separability of $g$ and Assumptions 2, 3, 6
and 4 ensure (see LABEL:{lemma:diff_prox}) that for any $j\in\hat{S}^{c}$,

$\displaystyle\partial_{z}\operatorname{prox}_{\gamma
g_{j}(\cdot,\lambda)}\left(\hat{\beta}_{j}-\gamma\nabla_{j}f(\hat{\beta})\right)=0\enspace.$
(13)

Plugging Equation 13 into Equation 12 ensures Equation 9 for all
$j\in\hat{S}^{c}$:

$\displaystyle\hat{\mathcal{J}}_{j:}=\partial_{\lambda}\operatorname{prox}_{\gamma
g_{j}(\cdot,\lambda)}\left(\hat{\beta}_{j}-\gamma\nabla_{j}f(\hat{\beta})\right)\enspace.$
(14)

Plugging Equations 13 and 14 into Equation 12 shows that the Jacobian
restricted on the generalized support $\hat{S}$ satisfies the following linear
system:

$\displaystyle\left(\text{Id}_{\hat{s}}-\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(\hat{z}\right)_{\hat{S}}\odot\big{(}\text{Id}_{\hat{s}}-\gamma\nabla^{2}_{\hat{S},\hat{S}}f(\hat{\beta})\big{)}\right)$
$\displaystyle\hat{\mathcal{J}}_{\hat{S}:}=$
$\displaystyle-\gamma\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})_{\hat{S}}\odot$
$\displaystyle\nabla^{2}_{\hat{S},\hat{S}^{c}}f(\hat{\beta})\hat{\mathcal{J}}_{\hat{S}^{c}:}+\partial_{\lambda}\operatorname{prox}_{g}(\hat{z})_{\hat{S}:}\enspace.$

Since $0<\gamma\leq 1/L$,

$\displaystyle\lVert\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})_{\hat{S}}\odot(\operatorname{Id}_{\hat{s}}-\gamma\nabla^{2}_{\hat{S},\hat{S}}f(\hat{\beta}))\rVert_{2}$
$\displaystyle\leq\lVert\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})_{\hat{S}}\rVert\cdot\lVert\operatorname{Id}_{\hat{s}}-\gamma\nabla^{2}_{\hat{S},\hat{S}}f(\hat{\beta})\rVert_{2}$
$\displaystyle<1\enspace.$ (15)

Since Equation 15 holds,
$A\triangleq\operatorname{Id}_{\hat{s}}-\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})_{\hat{S}}\odot(\operatorname{Id}_{\hat{s}}-\gamma\nabla^{2}_{\hat{S},\hat{S}}f(\hat{\beta}))$
is invertible, which leads to Equation 10.

###### Remark 10

In the smooth case a $p\times p$ linear system is needed to compute the
Jacobian in Equation 8. For non-smooth problems this is reduced to an
$\hat{s}\times\hat{s}$ linear system ($\hat{s}\leq 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 $\hat{\beta}^{(\lambda)}$.

###### 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
$\hat{\mathcal{J}}$ 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.

###  3.3 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:

$\displaystyle z_{j}$
$\displaystyle\leftarrow\beta_{j}-\gamma_{j}\nabla_{j}f(\beta)$
$\displaystyle\beta_{j}$
$\displaystyle\leftarrow\operatorname{prox}_{\gamma_{j}g_{j}}\left(\beta_{j}-\gamma_{j}\nabla_{j}f(\beta)\right)$
$\displaystyle\mathcal{J}_{j:}$
$\displaystyle\leftarrow\underbrace{\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j})}_{\in\mathbb{R}}\underbrace{\left(\mathcal{J}_{j:}-\gamma_{j}\nabla_{j:}^{2}f(\beta)\mathcal{J}\right)}_{\in\mathbb{R}^{p}}+\underbrace{\partial_{\lambda}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j})}_{\in\mathbb{R}^{p}}\enspace.$

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.

input : $X\in\mathbb{R}^{n\times
p},y\in\mathbb{R}^{n},\lambda\in\mathbb{R}^{r},n_{\mathrm{iter}}\in\mathbb{N}$,
$\beta\in\mathbb{R}^{p}$, $\mathcal{J}\in\mathbb{R}^{p\times
r},\gamma_{1},\dots,\gamma_{p}$

// jointly compute coef. & Jacobian

for  _$k=1,\dots,n_{\mathrm{iter}}$_ do

for  _$j=1,\ldots,p$_ do

// update the regression coefficients

$z_{j}\leftarrow\beta_{j}-\gamma_{j}\nabla_{j}f(\beta)$ ;

// CD step

$\mathrm{d}z_{j}\leftarrow\mathcal{J}_{j:}-\gamma_{j}\nabla^{2}_{j:}f(\beta)\mathcal{J}$

$\beta_{j}\leftarrow\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j})$
;

// prox. step

// update the Jacobian

// diff. with respect to $\lambda$

$\mathcal{J}_{j:}\leftarrow\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j})\mathrm{d}z_{j}$

$\mathcal{J}_{j:}\mathrel{+}=\partial_{\lambda}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j})$

$\beta^{(k)}=\beta$

$\mathcal{J}^{(k)}=\mathcal{J}$

$v=\nabla C(\beta)$

return _$\beta^{n_{\mathrm{iter}}},\mathcal{J}^{\top}v$_

Algorithm 3 Forward-mode PCD

input : $X\in\mathbb{R}^{n\times
p},y\in\mathbb{R}^{n},\lambda\in\mathbb{R}^{r},n_{\mathrm{iter}}\in\mathbb{N}$,
$\beta\in\mathbb{R}^{p},\gamma_{1},\dots,\gamma_{p}$

// compute coef.

for  _$k=1,\dots,n_{\mathrm{iter}}$_ do

for  _$j=1,\ldots,p$_ do

// update the regression coefficients

$z_{j}\leftarrow\beta_{j}-\gamma_{j}\nabla_{j}f(\beta)$ ;

// CD step

$\beta_{j}\leftarrow\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j})$
;

// prox. step

$\beta^{(k,j)}=\beta;z_{j}^{(k)}=z_{j}$ ;

// store iterates

// compute gradient $g$ in a backward way

$v=\nabla C(\beta^{n_{\mathrm{iter}}})$, $h=0_{\mathbb{R}^{r}}$

for  _$k=n_{\mathrm{iter}},n_{\mathrm{iter}}-1,\dots,1$_ do

for  _$j=p,\dots,1$_ do

$h\mathrel{-}=\gamma_{j}v_{j}\partial_{\lambda}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}\big{(}z_{j}^{(k)}\big{)}$

$v_{j}\mathrel{*}=\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}\big{(}z_{j}^{(k)}\big{)}$

$v\mathrel{-}=\gamma_{j}v_{j}\nabla_{j:}^{2}f(\beta^{(k,j)})$ ;

// $\mathcal{O}(np)$

return _ $\beta^{n_{\mathrm{iter}}},h$ _

Algorithm 4 Reverse-mode PCD

######  Theorem 12 (Local linear convergence of the Jacobian)

Let $0<\gamma\leq 1/L$. Suppose Assumptions 2, 3 and 6 hold. Let
$\lambda\in\mathbb{R}^{r}$, $\Lambda$ be a neighborhood of $\lambda$, and
$\Gamma^{\Lambda}\triangleq\left\\{\hat{\beta}^{(\lambda)}-\gamma\nabla
f(\hat{\beta}^{(\lambda)})\,:\,\lambda\in\Lambda\right\\}$. In addition,
suppose hypotheses (H1) to (H4) from Theorem 9 are satisfied and the sequence
$(\beta^{(k)})_{k\in\mathbb{N}}$ generated by Algorithm 1 (respectively by
Algorithm 3) converges toward $\hat{\beta}$.

Then, the sequence of Jacobians $(\mathcal{J}^{(k)})_{k\geq 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 $\hat{\mathcal{J}}$.

Proof of Theorem 12 can be found in Section B.

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

#### Comments on Figure 3.

We illustrate the results of Theorem 12 on SVM (for the Lasso and sparse
logistic regression, see Figures 10 and 11 in Section C) for multiple datasets
(_leukemia_ , _rcv1_ , _news20_ and _real-sim_ 444Data available on the
_libsvm_ website: https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/).
The values of the hyperparameters $\lambda$ are summarized in Table 6.
Regression coefficients $\hat{\beta}^{(\lambda)}$ 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 (Blondel and Pedregosa,
2016). The exact Jacobian was computed via implicit differentiation (Equation
10). 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 $\beta^{(k)}$ and
the exact solution $\hat{\beta}$. We also monitored the distance between the
iterates of the Jacobian $\mathcal{J}^{(k)}$ and the exact Jacobian
$\hat{\mathcal{J}}$. The red vertical dashed line represents the iteration
number where support identification happens. Once the support is identified,
Figures 10, 11 and 3 illustrate the linear convergence of the Jacobian.
However, the behavior of the iterative Jacobian before support identification
is more erratic and not even monotone.

###  3.4 Hypergradient computation with approximate gradients

As mentioned in Section 2, relying on iterative algorithms to solve Equation
1, one only has access to an approximation of $\hat{\beta}^{(\lambda)}$: 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 $\lambda\in\mathbb{R}^{r}$, let $\hat{\beta}^{(\lambda)}\in\mathbb{R}^{p}$
be the exact solution of the inner Equation 1, and $\hat{S}$ its generalized
support. Suppose Assumptions 2, 3 and 6 hold. Let $\Lambda$ be a neighborhood
of $\lambda$, and
$\Gamma^{\Lambda}\triangleq\left\\{\hat{\beta}^{(\lambda)}-\gamma\nabla
f(\hat{\beta}^{(\lambda)})\,:\,\lambda\in\Lambda\right\\}$. Suppose hypotheses
(H1) to (H4) from Theorem 9 are satisfied. In addition suppose

  1. (H5)

The application $\beta\mapsto\nabla^{2}f(\beta)$ is Lipschitz continuous.

  2. (H6)

The criterion $\beta\mapsto\nabla\mathcal{C}(\beta)$ is Lipschitz continuous.

  3. (H7)

Both optimization problems in Algorithm 5 are solved up to precision
$\epsilon$ with support identification:
$\lVert\beta^{(\lambda)}-\hat{\beta}^{(\lambda)}\rVert\leq\epsilon$,
$A^{\top}$ is invertible, and $\lVert
A^{-1\top}\nabla_{\hat{S}}\mathcal{C}(\beta^{(\lambda)})-v\rVert\leq\epsilon$.

Then the error on the approximate hypergradient $h$ returned by Algorithm 5 is
of the order of magnitude of the error $\epsilon$ on $\beta^{(\lambda)}$ and
$v$:

$\lVert\nabla\mathcal{L}(\lambda)-h\rVert=\mathcal{O}(\epsilon)\enspace.$

Proof of Theorem 13 can be found in Section B.1. Following the analysis of
Pedregosa (2016), two sources of approximation errors arise when computing the
hypergradient: one from the inexact computation of $\hat{\beta}$, 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 $\epsilon$, i.e.,
$\lVert\hat{\beta}^{(\lambda)}-\beta^{(\lambda)}\rVert\leq\epsilon$ and
$\lVert
A^{-1\top}\nabla_{S}\mathcal{C}(\beta^{(\lambda)})-v\rVert\leq\epsilon$, then
the approximation on the hypergradient is also of the order of $\epsilon$.

###### Remark 14

The Lipschitz continuity of the proximity operator with respect to $\lambda$
(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 Equation 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 $\lVert\beta^{(k)}-\hat{\beta}\rVert$: most
softwares provide a solution up to a given duality gap (sometimes even other
criteria), not $\lVert\beta^{(k)}-\hat{\beta}\rVert$.

###  3.5 Proposed method for hypergradient computation

We now describe our proposed method to compute the hypergradient of Equation
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 Equation 1 using a solver identifying the generalized support (Liang et
al., 2014; Klopfenstein et al., 2020). Then, the hypergradient is computed by
solving the linear system in Equation 10. 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 $\hat{s}$ is the size of the generalized support (Definition 5, $\hat{s}\leq p$ and usually $\hat{s}\ll p$). The number of iterations of the inner solver is noted $n_{\mathrm{iter}}$, the number of iterations of the solver of the linear system is noted $n_{\text{sys}}$. Differentiation | Algorithm | Space | Time
---|---|---|---
Forward-mode PGD | Algorithm 1 | $\mathcal{O}(p\,r)$ | $\mathcal{O}(n\,p\,r\,n_{\mathrm{iter}})$
Reverse-mode PGD | Algorithm 2 | $\mathcal{O}(p\,n_{\mathrm{iter}})$ | $\mathcal{O}(n\,p\,n_{\mathrm{iter}}+n\,p\,n_{\mathrm{iter}})$
Forward-mode PCD | Algorithm 3 | $\mathcal{O}(p\,r)$ | $\mathcal{O}(n\,p\,r\,n_{\mathrm{iter}})$
Reverse-mode PCD | Algorithm 4 | $\mathcal{O}(p\,n_{\mathrm{iter}})$ | $\mathcal{O}(n\,p\,n_{\mathrm{iter}}+n\,p^{2}\,n_{\mathrm{iter}})$
Implicit differentiation | Algorithm 5 | $\mathcal{O}(p+\hat{s})$ | $\mathcal{O}(n\,p\,n_{\mathrm{iter}}+n\,\hat{s}\,n_{\mathrm{sys}})$

###  3.6 Resolution of the bilevel optimization Equation 2

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 Equation 2. As the Lipschitz constant is not available for the outer
problem, first-order methods need to rely on other strategies, such as:

  * •

Gradient descent with manually adjusted fixed step sizes (Frecon et al., 2018;
Ji et al., 2020). 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 Deledalle et al. 2014). 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 (Nocedal and Wright, 2006).
Yet, due to the approximate gradient computation, we observed that L-BFGS did
not always converge.

  * •

ADAM (Kingma and Ba, 2014). 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 (Pedregosa, 2016). While
the approach from Pedregosa (2016) requires no tuning, we observed that it
could diverge when close to the optimum. The adaptive step size strategy
proposed in Algorithm 6, used in all the experiments, turned out to be robust
and efficient across problems and datasets.

######  Remark 16 (Uniqueness)

The solution of Equation 1 may be non-unique, leading to a multi-valued
regularization path $\lambda\mapsto\hat{\beta}^{(\lambda)}$ (Liu et al., 2020)
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 (Tibshirani, 2013):
when the design matrix is sampled from a continuous distribution, the solution
of the Lasso is almost surely unique.

######  Remark 17 (Initialization)

One advantage of the non-smooth case with the $\ell_{1}$ norm is that one can
find a good initialization point: there exists a value $\lambda_{\max}$ (see
Table 1) such that the solution of Equation 1 vanishes for
$\lambda\geq\lambda_{\max}$. Hence, a convenient and robust initialization
value can be chosen as $e^{\lambda}=e^{\lambda_{\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 (Pedregosa, 2016).

input : $\lambda\in\mathbb{R},\epsilon>0$

init : $\gamma>0$

;

// compute the solution of inner problem

Find $\beta$ such that:
$\Phi(\beta,\lambda)-\Phi(\hat{\beta},\lambda)\leq\epsilon$

;

// compute the gradient

Compute the generalized support $S$ of $\beta$,

$z={\beta}-\gamma\nabla f(\beta)$

$\mathcal{J}_{S^{c}:}=\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z)_{S^{c}}$

$s=|S|$

$A=\operatorname{Id}_{s}-\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z)_{S}\odot(\operatorname{Id}_{s}-\gamma\nabla^{2}_{S,S}f(\beta))$

Find $v\in\mathbb{R}^{s}$ s.t. $\lVert
A^{-1\top}\nabla_{S}\mathcal{C}(\beta)-v\rVert\leq\epsilon$

$None$

-γ∂zproxγg(⋅,λ)(z)S⊙∇2S,Scf(β)JSc

$\nabla\mathcal{L}(\lambda)=\mathcal{J}_{S^{c}:}^{\top}\nabla_{S^{c}}\mathcal{C}(\beta)+v^{\top}B$

return _
$\mathcal{L}(\lambda)\triangleq\mathcal{C}(\beta),\nabla\mathcal{L}(\lambda)$
_

Algorithm 5 Implicit differentiation

input : $\lambda\in\mathbb{R}^{r},(\epsilon_{i})$

init : $\textrm{use\\_adaptive\\_step\\_size}=\textrm{True}$

for  _$i=1,\dots,\mathrm{iter}$_ do

$\lambda^{\mathrm{old}}\leftarrow\lambda$

// compute the value and the gradient
$\mathcal{L}(\lambda),\nabla\mathcal{L}(\lambda)\leftarrow{\rm\lx@cref{creftypecap~refnum}{alg:implicit}}(X,y,\lambda,\epsilon_{i})$

if  _$\mathrm{use\\_adaptive\\_step\\_size}$_ then

$\alpha=1/\lVert\nabla\mathcal{L}(\lambda)\rVert$

$\lambda\mathrel{-}=\alpha\nabla\mathcal{L}(\lambda)$ ;

// gradient step

if  _$\mathcal{L}(\lambda) >\mathcal{L}(\lambda^{\mathrm{old}})$_ then

$\mathrm{use\\_adaptive\\_step\\_size}=\mathrm{False}$

$\alpha\mathrel{/}=10$

return _ $\lambda$ _

Algorithm 6 Gradient descent with approximate gradient

##  4 Experiments

In this section, we illustrate the benefits of our proposed Algorithm 5 to
compute hypergradients and Algorithm 6 to solve Equation 2. Our package,
sparse-ho, is implemented in Python. It relies on Numpy (Harris et al., 2020),
Numba (Lam et al., 2015) and SciPy (Virtanen et al., 2020). Figures were
plotted using matplotlib (Hunter, 2007). 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.

Table 5:  Characteristics of the datasets used for the experiments. name | $\\#\text{ samples }n$ | $\\#\text{ features }p$ | $\\#\text{ classes }q$ | density
---|---|---|---|---
_breast cancer_ | $569$ | $30$ | $-$ | $1$
_diabetes_ | $442$ | $10$ | $-$ | $1$
_leukemia_ | $72$ | $7129$ | $-$ | $1$
_gina agnostic_ | $3468$ | $970$ | $-$ | $1$
_rcv1_ | $20242$ | $19960$ | $-$ | $3.7\times 10^{-3}$
_real-sim_ | $72309$ | $20958$ | $-$ | $2.4\times 10^{-3}$
_news20_ | $19996$ | $632983$ | $-$ | $6.1\times 10^{-4}$
_mnist_ | $60.000$ | $683$ | $10$ | $2.2\times 10^{-1}$
_usps_ | $7291$ | $256$ | $10$ | $1$
_rcv1 (multiclass)_ | $15564$ | $16245$ | $53$ | $4.0\times 10^{-3}$
_aloi_ | $108000$ | $128$ | $1000$ | $2.4\times 10^{-1}$

###  4.1 Hypergradient computation

#### Comparison with alternative approaches (Figure 4).

Figure 4:  Lasso with hold-out criterion: time comparison on the gina dataset
to compute a single hypergradient as a function of the number of features, for
two values of $\lambda$, $e^{\lambda}=e^{\lambda_{\max}}/10$ (left) and
$e^{\lambda}=e^{\lambda_{\max}}/100$ (right).

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 (Agrawal et al., 2019), a software based on cvxpy (Diamond and
Boyd, 2016), 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 $\nabla\mathcal{L}(\lambda)$ for the Lasso (see Table 1), for
multiple values of $\lambda$. 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_ dataset. 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 Equation 1 using a splitting conic solver
(O’Donoghue et al., 2019). 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
$\nabla\mathcal{L}(\lambda)$ 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 differentiation of proximal coordinate descent. The larger
the value of $\lambda$, the sparser the coefficients $\beta$ 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 and 6).

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 (Lam
et al., 2015).

  * •

Implicit differentiation (Algorithm 5) with state-of-the-art algorithm to
solve the inner problem: we used Celer (Massias et al., 2020) for the Lasso,
and Lightning (Blondel and Pedregosa, 2016) for the SVM.

Figure 5:  Lasso with hold-out criterion: absolute difference between the
exact hypergradient (using $\hat{\beta}$) and the iterate hypergradient (using
$\beta^{(k)}$) of the Lasso as a function of time. Results are for three
datasets and two different regularization parameters. “Implicit diff. +
Celer)” uses Celer (Massias et al., 2020) instead of our proximal coordinate
descent implementation.

Figure 6:  SVM with hold-out criterion: absolute difference between the exact
hypergradient (using $\hat{\beta}$) and the iterate hypergradient (using
$\beta^{(k)}$) of the SVM as a function of time. “Implicit diff. + Lightning”
uses Lightning (Blondel and Pedregosa, 2016), instead of our proximal
coordinate descent implementation.

Figure 5 shows for three datasets 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 datasets, except _news20_ ,
which is not well suited for SVM, due to limited number of samples.

First, it demonstrates that implicit differentiation methods are faster than
the forward-mode 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).

###  4.2 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 $\lambda$, uniformly spaced between
$\lambda_{\max}-\ln(10^{4})$ and $\lambda_{\max}$. For the elastic net we
chose for each of the two hyperparameters a grid of 10 values uniformly spaced
between $\lambda_{\max}$ and $\lambda_{\max}-\ln(10^{4})$. The product grid
thus has $10^{2}$ points.

  * •

Random-search: we chose $30$ values of $\lambda$ sampled uniformly between
$\lambda_{\max}$ and $\lambda_{\max}-\ln(10^{4})$ for each hyperparameter. For
the elastic net we chose $30$ points sampled uniformly in
$[\lambda_{\max}-\ln(10^{4}),\lambda_{\max}]\times[\lambda_{\max}-\ln(10^{4}),\lambda_{\max}]$.

  * •

SMBO: this algorithm is SMBO using as criterion expected improvement (EI) and
the Tree-structured Parzen Estimator (TPE) as model. First it evaluates
$\mathcal{L}$ using $5$ values of $\lambda$, chosen uniformly at random
between $\lambda_{\max}$ and $\lambda_{\max}-\ln(10^{4})$. Then a TPE model is
fitted on the data points
$(\lambda^{(1)},\mathcal{L}(\lambda^{(1)})),\dots,(\lambda^{(5)},\mathcal{L}(\lambda^{(5)}))$.
Iteratively, the EI is used to choose the next point to evaluate $\mathcal{L}$
at, and this value is used to update the model. We used the hyperopt
implementation (Bergstra et al., 2013).

  * •

1st order: first-order method with exact gradient (Algorithm 6 with constant
tolerances $\epsilon_{i}=10^{-6}$), with $\lambda_{\max}-\ln(10^{2})$ as a
starting point.

  * •

1st order approx: a first-order method using approximate gradient (Algorithm 6
with tolerances $\epsilon_{i}$, geometrically decreasing from $10^{-2}$ to
$10^{-6}$), with $\lambda_{\max}-\ln(10^{2})$ as a starting point.

Outer criterion. In the Lasso and elastic net experiments, we pick a $K$-fold
CV loss as outer criterion555In our experiments the default choice is $K=5$..
Hence, the dataset $(X,y)$ is partitioned into $K$ hold-out datasets
$(X^{\text{train}_{k}},y^{\text{train}_{k}}),(X^{\text{val}_{k}},y^{\text{val}_{k}})$.
The bilevel optimization problems then write:

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

$\displaystyle\operatorname*{\mathrm{arg\,min}}_{\lambda=(\lambda_{1},\lambda_{2})\in\mathbb{R}^{2}}\mathcal{L}(\lambda)=\frac{1}{K}\sum_{k=1}^{K}\lVert
y^{\text{val}_{k}}-X^{\text{val}_{k}}\hat{\beta}^{(\lambda,k)}\rVert^{2}_{2}$
(16) $\displaystyle{s.t.\leavevmode\nobreak\
}\hat{\beta}^{(\lambda,k)}\in\operatorname*{\mathrm{arg\,min}}_{\beta\in\mathbb{R}^{p}}\tfrac{1}{2n}\left\lVert
y^{\text{train}_{k}}-X^{\text{train}_{k}}\beta\right\rVert^{2}_{2}+e^{\lambda_{1}}\lVert\beta\rVert_{1}+\frac{e^{\lambda_{2}}}{2}\lVert\beta\rVert_{2}^{2},\quad\forall
k\in[K]\enspace,$

while Lasso CV is obtained taking $\lambda_{2}\to-\infty$ in the former. By
considering an extended variable $\beta\in\mathbb{R}^{K\times p}$, cross-
validation can be cast as an instance of Equation 2.

Figure 7 represents the cross-validation loss in Lasso CV as a function of the
regularization parameter $\lambda$ (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 $\lambda$ value. The top rows show cross-
validation loss as a function of $\lambda$, 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 $\lambda=\lambda_{\max}$ and then decreases to
$\lambda=\lambda_{\max}-\ln(10^{4})$. On all the datasets, first-order methods
are faster to find the optimal regularization parameter, requiring only $5$
iterations.

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_ datasets.

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\in\mathbb{R}^{n\times p}$, and the target variable
$y\in\\{1,\dots,q\\}^{n}$. We chose to use a one-versus-all model with $q$
regularization parameters. We use a binary cross-entropy for the inner loss:

$\displaystyle\psi^{k}(\beta,\lambda_{k};X,y)\triangleq-\frac{1}{n}\sum_{i=1}^{n}\left(\mathbbm{1}_{y_{i}=k}\ln(\sigma(X_{i:}\beta))+(1-\mathbbm{1}_{y_{i}=k})\ln(1-\sigma(X_{i:}\beta))\right)+e^{\lambda_{k}}\lVert\beta\rVert_{1}\enspace,$

and a multiclass cross-entropy for the outer criterion:

$\displaystyle\mathcal{C}\left(\hat{\beta}^{(\lambda_{1})},\dots,\hat{\beta}^{(\lambda_{q})};X,y\right)\triangleq-\sum_{i=1}^{n}\sum_{k=1}^{q}\ln\left(\frac{e^{X_{i:}\hat{\beta}^{(\lambda_{k})}}}{\sum_{l=1}^{q}e^{X_{i:}\hat{\beta}^{(\lambda_{l})}}}\right)\mathbbm{1}_{y_{i}=k}\enspace.$
(17)

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

$\displaystyle\operatorname*{\mathrm{arg\,min}}_{\lambda\triangleq(\lambda_{1},\dots,\lambda_{q})\in\mathbb{R}^{q}}$
$\displaystyle\mathcal{C}\left(\hat{\beta}^{(\lambda_{1})},\dots,\hat{\beta}^{(\lambda_{q})};X^{\mathrm{test}},y^{\mathrm{test}}\right)$
(18) $\displaystyle{s.t.\leavevmode\nobreak}$
$\displaystyle\hat{\beta}^{(\lambda_{k})}\in\operatorname*{\mathrm{arg\,min}}_{\beta\in\mathbb{R}^{p}}\psi^{k}(\beta,\lambda_{k};X^{\mathrm{train}},y^{\mathrm{train}})\quad\forall
k\in[q]\enspace.$

Figure 9:  Multiclass sparse logistic regression hold-out, time comparison
($\\#$ classes hyperparameters). Multiclass cross-entropy (top), accuracy on
the validation set (middle), and accuracy on the test set (bottom) as a
function of time on _mnist_ , _usps_ ($q=10$ classes), _rcv1_ ($q=53$
classes), _aloi_ ($q=1000$ classes).

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 Bergstra and Bengio (2012); Frazier (2018): 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.

##  5 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 ($\ell_{1}$ penalized
methods, SVM, etc.) and datasets. 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.

## Acknowledgements

This work was partially funded by the ERC Starting Grant SLAB ERC-StG-676943 ,
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).

## References

  * Ablin et al. (2020) P. Ablin, G. Peyré, and T. Moreau.  Super-efficiency of automatic differentiation for functions defined as a minimum.  In _International Conference on Machine Learning_ , pages 32–41. PMLR, 2020.
  * Agrawal et al. (2019) A. Agrawal, B. Amos, S. Barratt, S. Boyd, S. Diamond, and J. Z. Kolter.  Differentiable convex optimization layers.  In _NeurIPS_ , pages 9558–9570, 2019.
  * Akaike (1974) H. Akaike.  A new look at the statistical model identification.  _IEEE Trans. Automat. Control_ , AC-19:716–723, 1974.
  * Amos and Kolter (2017) B. Amos and J. Z. Kolter.  Optnet: Differentiable optimization as a layer in neural networks.  In _ICML_ , volume 70, pages 136–145, 2017.
  * Arlot and Celisse (2010) S. Arlot and A. Celisse.  A survey of cross-validation procedures for model selection.  _Statistics surveys_ , 4:40–79, 2010.
  * Bai et al. (2019) S. Bai, J. Z. Kolter, and V. Koltun.  Deep equilibrium models.  _NeurIPS_ , 2019.
  * Bai et al. (2020) S. Bai, V. Koltun, and J. Z. Kolter.  Multiscale deep equilibrium models.  _NeurIPS_ , 2020.
  * Baydin et al. (2018) A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind.  Automatic differentiation in machine learning: a survey.  _Journal of Machine Learning Research_ , 18(153):1–43, 2018.
  * Belloni et al. (2011) A. Belloni, V. Chernozhukov, and L. Wang.  Square-root Lasso: pivotal recovery of sparse signals via conic programming.  _Biometrika_ , 98(4):791–806, 2011.
  * Bengio (2000) Y. Bengio.  Gradient-based optimization of hyperparameters.  _Neural computation_ , 12(8):1889–1900, 2000\.
  * Bergstra and Bengio (2012) J. Bergstra and Y. Bengio.  Random search for hyper-parameter optimization.  _Journal of Machine Learning Research_ , 13(2), 2012.
  * Bergstra et al. (2011) J. Bergstra, R. Bardenet, Y. Bengio, and B. Kégl.  Algorithms for hyper-parameter optimization.  In _NeurIPS_ , 2011.
  * Bergstra et al. (2013) J. Bergstra, D. Yamins, and D. D. Cox.  Hyperopt: A python library for optimizing the hyperparameters of machine learning algorithms.  In _Proceedings of the 12th Python in science conference_ , pages 13–20, 2013.
  * Bertrand et al. (2020) Q. Bertrand, Q. Klopfenstein, M. Blondel, S. Vaiter, A. Gramfort, and J. Salmon.  Implicit differentiation of Lasso-type models for hyperparameter optimization.  _ICML_ , 2020.
  * Bickel et al. (2009) P. J. Bickel, Y. Ritov, and A. B. Tsybakov.  Simultaneous analysis of Lasso and Dantzig selector.  _Ann. Statist._ , 37(4):1705–1732, 2009.
  * Blondel and Pedregosa (2016) M. Blondel and F. Pedregosa.  Lightning: large-scale linear classification, regression and ranking in python, 2016.
  * Bolte and Pauwels (2020a) J. Bolte and E. Pauwels.  Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning.  _Mathematical Programming_ , pages 1–33, 2020a.
  * Bolte and Pauwels (2020b) J. Bolte and E. Pauwels.  A mathematical model for automatic differentiation in machine learning.  _arXiv preprint arXiv:2006.02080_ , 2020b.
  * Boser et al. (1992) B. E. Boser, I. M. Guyon, and V. N. Vapnik.  A training algorithm for optimal margin classifiers.  In _Proceedings of the fifth annual workshop on Computational learning theory_ , pages 144–152. ACM, 1992.
  * Brochu et al. (2010) E. Brochu, V. M. Cora, and N. De Freitas.  A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning.  2010\.
  * Chapelle et al. (2002) O. Chapelle, V. Vapnik, O. Bousquet, and S. Mukherjee.  Choosing multiple parameters for support vector machines.  _Machine learning_ , 46(1-3):131–159, 2002.
  * Chen et al. (1998) S. S. Chen, D. L. Donoho, and M. A. Saunders.  Atomic decomposition by basis pursuit.  _SIAM J. Sci. Comput._ , 20(1):33–61, 1998.
  * Cherkaoui et al. (2020) H. Cherkaoui, J. Sulam, and T. Moreau.  Learning to solve TV regularised problems with unrolled algorithms.  _NeurIPS_ , 33, 2020.
  * Colson et al. (2007) B. Colson, P. Marcotte, and G. Savard.  An overview of bilevel optimization.  _Annals of operations research_ , 153(1):235–256, 2007.
  * Combettes and Wajs (2005) P. L. Combettes and V. R. Wajs.  Signal recovery by proximal forward-backward splitting.  _Multiscale Modeling & Simulation_, 4(4):1168–1200, 2005.
  * Defazio et al. (2014) A. Defazio, F. Bach, and S. Lacoste-Julien.  Saga: A fast incremental gradient method with support for non-strongly convex composite objectives.  In _NeurIPS_ , pages 1646–1654, 2014.
  * Deledalle et al. (2014) C.-A. Deledalle, S. Vaiter, J. Fadili, and G. Peyré.  Stein Unbiased GrAdient estimator of the Risk (SUGAR) for multiple parameter selection.  _SIAM J. Imaging Sci._ , 7(4):2448–2487, 2014\.
  * Dempe et al. (2015) S. Dempe, V. Kalashnikov, G. A. Pérez-Valdés, and N. Kalashnykova.  Bilevel programming problems.  _Energy Systems. Springer, Berlin_ , 2015.
  * Devroye and Wagner (1979) L. Devroye and T. Wagner.  Distribution-free performance bounds for potential function rules.  _IEEE Transactions on Information Theory_ , 25(5):601–604, 1979.
  * Diamond and Boyd (2016) S. Diamond and S. Boyd.  CVXPY: A Python-embedded modeling language for convex optimization.  _Journal of Machine Learning Research_ , 17(83):1–5, 2016.
  * Domke (2012) J. Domke.  Generic methods for optimization-based modeling.  In _AISTATS_ , volume 22, pages 318–326, 2012.
  * Efron (1986) B. Efron.  How biased is the apparent error rate of a prediction rule?  _J. Amer. Statist. Assoc._ , 81(394):461–470, 1986.
  * Evans and Gariepy (1992) L. C. Evans and R. F. Gariepy.  _Measure theory and fine properties of functions_.  CRC Press, 1992.
  * Fan and Lv (2008) J. Fan and J. Lv.  Sure independence screening for ultrahigh dimensional feature space.  _J. R. Stat. Soc. Ser. B Stat. Methodol._ , 70(5):849–911, 2008.
  * Feurer and Hutter (2019) M. Feurer and F. Hutter.  Hyperparameter optimization.  In _Automated Machine Learning_ , pages 3–33. Springer, Cham, 2019\.
  * Foo et al. (2008) C. S. Foo, C. B. Do, and A. Y. Ng.  Efficient multiple hyperparameter learning for log-linear models.  In _NeurIPS_ , pages 377–384, 2008.
  * Forrester et al. (2008) A. Forrester, A. Sobester, and A. Keane.  _Engineering design via surrogate modelling: a practical guide_.  John Wiley & Sons, 2008.
  * Franceschi et al. (2017) L. Franceschi, M. Donini, P. Frasconi, and M. Pontil.  Forward and reverse gradient-based hyperparameter optimization.  In _ICML_ , pages 1165–1173, 2017.
  * Franceschi et al. (2018) L. Franceschi, P. Frasconi, S. Salzo, and M. Pontil.  Bilevel programming for hyperparameter optimization and meta-learning.  In _ICML_ , pages 1563–1572, 2018.
  * Frazier (2018) P.I. Frazier.  A tutorial on Bayesian optimization.  _arXiv preprint arXiv:1807.02811_ , 2018.
  * Frecon et al. (2018) J. Frecon, S. Salzo, and M. Pontil.  Bilevel learning of the group lasso structure.  In _NeurIPS_ , pages 8301–8311, 2018.
  * Friedman et al. (2010) J. Friedman, T. J. Hastie, and R. Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _J. Stat. Softw._ , 33(1):1–22, 2010.
  * Ghadimi and Wang (2018) S. Ghadimi and M. Wang.  Approximation methods for bilevel programming.  _arXiv preprint arXiv:1802.02246_ , 2018.
  * Goodfellow et al. (2016) I. Goodfellow, A. Courville, and Y. Bengio.  _Deep learning_ , volume 1.  MIT press Cambridge, 2016.
  * Gould et al. (2016) S. Gould, B. Fernando, A. Cherian, P. Anderson, R. S. Cruz, and E. Guo.  On differentiating parameterized argmin and argmax problems with application to bi-level optimization.  _arXiv preprint arXiv:1607.05447._ , 2016.
  * Grazzi et al. (2020) R. Grazzi, L. Franceschi, M. Pontil, and S. Salzo.  On the iteration complexity of hypergradient computation.  _ICML_ , 2020.
  * Hare and Lewis (2007) W. L. Hare and A. S. Lewis.  Identifying active manifolds.  _Algorithmic Operations Research_ , 2(2):75–75, 2007.
  * Harris et al. (2020) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. Fern’andez del R’ıo, M. Wiebe, P. Peterson, P. G’erard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, 2020.
  * Hastie et al. (2004) T. Hastie, S. Rosset, R. Tibshirani, and J. Zhu.  The entire regularization path for the support vector machine.  _Journal of Machine Learning Research_ , 5(Oct):1391–1415, 2004.
  * Hestenes and Stiefel (1952) M. R. Hestenes and E. Stiefel.  _Methods of conjugate gradients for solving linear systems_ , volume 49.  NBS Washington, DC, 1952.
  * Higham (2002) N. J. Higham.  _Accuracy and stability of numerical algorithms_.  SIAM, 2002.
  * Hoerl and Kennard (1970) A. E. Hoerl and R. W. Kennard.  Ridge regression: Biased estimation for nonorthogonal problems.  _Technometrics_ , 12(1):55–67, 1970.
  * Hunter (2007) J. D. Hunter.  Matplotlib: A 2d graphics environment.  _IEEE Annals of the History of Computing_ , 9(03):90–95, 2007.
  * Hutter et al. (2015) F. Hutter, J. Lücke, and L. Schmidt-Thieme.  Beyond manual tuning of hyperparameters.  _KI-Künstliche Intelligenz_ , 29(4):329–337, 2015.
  * Ji et al. (2020) K. Ji, J. Yang, and Y. Liang.  Provably faster algorithms for bilevel optimization and applications to meta-learning.  _arXiv preprint arXiv:2010.07962_ , 2020.
  * Johnson and Zhang (2013) R. Johnson and T. Zhang.  Accelerating stochastic gradient descent using predictive variance reduction.  In _NeurIPS_ , pages 315–323, 2013.
  * Johnson and Guestrin (2015) T. B. Johnson and C. Guestrin.  Blitz: A principled meta-algorithm for scaling sparse optimization.  In _ICML_ , pages 1171–1179, 2015.
  * Jones et al. (1998) D. R. Jones, M. Schonlau, and W. J. Welch.  Efficient global optimization of expensive black-box functions.  _Journal of Global optimization_ , 13(4):455–492, 1998.
  * Kingma and Ba (2014) D. P. Kingma and J. L. Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Klopfenstein et al. (2020) Q. Klopfenstein, Q. Bertrand, A. Gramfort, J. Salmon, and S. Vaiter.  Model identification and local linear convergence of coordinate descent.  _arXiv preprint arXiv:2010.11825_ , 2020.
  * Koh et al. (2007) K. Koh, S.-J. Kim, and S. Boyd.  An interior-point method for large-scale l1-regularized logistic regression.  _Journal of Machine Learning Research_ , 8(8):1519–1555, 2007.
  * Kohavi and John (1995) R. Kohavi and G. H. John.  Automatic parameter selection by minimizing estimated error.  In _Machine Learning Proceedings 1995_ , pages 304–312. Elsevier, 1995.
  * Kunisch and Pock (2013) K. Kunisch and T. Pock.  A bilevel optimization approach for parameter learning in variational models.  _SIAM J. Imaging Sci._ , 6(2):938–983, 2013.
  * Lam et al. (2015) S. K. Lam, A. Pitrou, and S. Seibert.  Numba: A LLVM-based Python JIT Compiler.  In _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_ , pages 1–6. ACM, 2015.
  * Larsen et al. (1996) J. Larsen, L. K. Hansen, C. Svarer, and M. Ohlsson.  Design and regularization of neural networks: the optimal use of a validation set.  In _Neural Networks for Signal Processing VI. Proceedings of the 1996 IEEE Signal Processing Society Workshop_ , 1996.
  * LeCun et al. (1998) Y. A. LeCun, L. Bottou, G. B. Orr, and K-R. Müller.  Efficient backprop.  In _Neural networks: Tricks of the trade_ , pages 9–48. Springer, 1998.
  * Liang et al. (2014) J. Liang, J. Fadili, and G. Peyré.  Local linear convergence of forward–backward under partial smoothness.  In _NeurIPS_ , pages 1970–1978, 2014.
  * Liang et al. (2017) J. Liang, J. Fadili, and G. Peyré.  Activity identification and local linear convergence of Forward–Backward-type Methods.  _SIAM J. Optim._ , 27(1):408–437, 2017.
  * Linnainmaa (1970) S. Linnainmaa.  The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors.  _Master’s Thesis (in Finnish), Univ. Helsinki_ , pages 6–7, 1970\.
  * Lions and Mercier (1979) P-L. Lions and B. Mercier.  Splitting algorithms for the sum of two nonlinear operators.  _SIAM Journal on Numerical Analysis_ , 16(6):964–979, 1979.
  * Liu and Nocedal (1989) D. C. Liu and J. Nocedal.  On the limited memory BFGS method for large scale optimization.  _Mathematical programming_ , 45(1-3):503–528, 1989.
  * Liu et al. (2020) R. Liu, P. Mu, X. Yuan, S. Zeng, and J. Zhang.  A generic first-order algorithmic framework for bi-level programming beyond lower-level singleton.  _ICML_ , 2020.
  * Lorraine et al. (2019) J. Lorraine, P. Vicol, and D. Duvenaud.  Optimizing millions of hyperparameters by implicit differentiation.  _arXiv preprint arXiv:1911.02590_ , 2019.
  * Lounici (2008) K. Lounici.  Sup-norm convergence rate and sign concentration property of Lasso and Dantzig estimators.  _Electron. J. Stat._ , 2:90–102, 2008.
  * Mairal and Yu (2012) J. Mairal and B. Yu.  Complexity analysis of the lasso regularization path.  In _ICML_ , pages 353–360, 2012.
  * Mairal et al. (2012) J. Mairal, F. Bach, and J. Ponce.  Task-driven dictionary learning.  _IEEE Trans. Pattern Anal. Mach. Intell._ , 34(4):791–804, 2012.
  * Massias et al. (2018) M. Massias, A. Gramfort, and J. Salmon.  Celer: a fast solver for the lasso with dual extrapolation.  In _ICML_ , volume 80, pages 3315–3324, 2018.
  * Massias et al. (2020) M. Massias, S. Vaiter, A. Gramfort, and J. Salmon.  Dual extrapolation for sparse generalized linear models.  _Journal of Machine Learning Research_ , 21(234):1–33, 2020.
  * Mehmood and Ochs (2021) S. Mehmood and P. Ochs.  Differentiating the value function by using convex duality.  In _AISTATS_ , pages 3871–3879. PMLR, 2021.
  * Mockus (1989) J. Mockus.  The bayesian approach to local optimization.  In _Bayesian Approach to Global Optimization_ , pages 125–156. Springer, 1989.
  * Nesterov (2004) Y. Nesterov.  _Introductory lectures on convex optimization_ , volume 87 of _Applied Optimization_.  Kluwer Academic Publishers, Boston, MA, 2004.
  * Nocedal and Wright (2006) J. Nocedal and S. J. Wright.  _Numerical optimization_.  Springer Series in Operations Research and Financial Engineering. Springer, New York, second edition, 2006.
  * Nutini (2018) J. Nutini.  _Greed is good: greedy optimization methods for large-scale structured problems_.  PhD thesis, University of British Columbia, 2018.
  * Nutini et al. (2019) J. Nutini, M. Schmidt, and W. Hare.  “active-set complexity” of proximal gradient: How long does it take to find the sparsity pattern?  _Optimization Letters_ , 13(4):645–655, 2019\.
  * Ochs et al. (2015) P. Ochs, R. Ranftl, T. Brox, and T. Pock.  Bilevel optimization with nonsmooth lower level problems.  In _SSVM_ , pages 654–665, 2015.
  * O’Donoghue et al. (2016) B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd.  Conic optimization via operator splitting and homogeneous self-dual embedding.  _Journal of Optimization Theory and Applications_ , 169(3):1042–1068, 2016.
  * O’Donoghue et al. (2019) B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd.  SCS: Splitting conic solver, version 2.1.2, 2019.
  * Pedregosa (2016) F. Pedregosa.  Hyperparameter optimization with approximate gradient.  In _ICML_ , volume 48, pages 737–746, 2016.
  * Pedregosa et al. (2017) F. Pedregosa, R. Leblond, and S. Lacoste-Julien.  Breaking the nonsmooth barrier: A scalable parallel method for composite optimization.  _NeurIPS_ , pages 56–65, 2017.
  * Peyré and Fadili (2011) G. Peyré and J. M. Fadili.  Learning analysis sparsity priors.  In _Sampta_ , 2011.
  * Platt (1999) J. C. Platt.  Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods.  In _Advances in large margin classifiers_ , pages 61–74. MIT Press, 1999.
  * Poliquin and Rockafellar (1996) R. A. Poliquin and R. T. Rockafellar.  Generalized hessian properties of regularized nonsmooth functions.  _SIAM Journal on Optimization_ , 6(4):1121–1137, 1996.
  * Polyak (1987) B. T. Polyak.  Introduction to optimization.  _Inc., Publications Division, New York_ , 1, 1987.
  * Rajeswaran et al. (2019) A. Rajeswaran, C. Finn, S. M. Kakade, and S. Levine.  Meta-learning with implicit gradients.  In _NeurIPS_ , pages 113–124, 2019.
  * Rastrigin (1963) L. A. Rastrigin.  The convergence of the random search method in the extremal control of a many parameter system.  _Automaton & Remote Control_, 24:1337–1342, 1963.
  * Rosset and Zhu (2007) S. Rosset and J. Zhu.  Piecewise linear regularized solution paths.  _Ann. Statist._ , 35(3):1012–1030, 2007.
  * Schwarz (1978) G. Schwarz.  Estimating the dimension of a model.  _Ann. Statist._ , 6(2):461–464, 1978.
  * Seeger (2008) M. W. Seeger.  Cross-validation optimization for large scale structured classification kernel methods.  _Journal of Machine Learning Research_ , 9:1147–1178, 2008\.
  * Snoek et al. (2012) J. Snoek, H. Larochelle, and R. P. Adams.  Practical bayesian optimization of machine learning algorithms.  In _NeurIPS_ , pages 2960–2968, 2012.
  * Stein (1981) C. M. Stein.  Estimation of the mean of a multivariate normal distribution.  _Ann. Statist._ , 9(6):1135–1151, 1981.
  * Stone and Ramer (1965) L. R. A. Stone and J.C. Ramer.  Estimating WAIS IQ from Shipley Scale scores: Another cross-validation.  _Journal of clinical psychology_ , 21(3):297–297, 1965.
  * Tibshirani (1996) R. Tibshirani.  Regression shrinkage and selection via the lasso.  _J. R. Stat. Soc. Ser. B Stat. Methodol._ , 58(1):267–288, 1996.
  * Tibshirani et al. (2012) R. Tibshirani, J. Bien, J. Friedman, T. J. Hastie, N. Simon, J. Taylor, and R. J. Tibshirani.  Strong rules for discarding predictors in lasso-type problems.  _J. R. Stat. Soc. Ser. B Stat. Methodol._ , 74(2):245–266, 2012.
  * Tibshirani (2013) R. J. Tibshirani.  The lasso problem and uniqueness.  _Electron. J. Stat._ , 7:1456–1490, 2013.
  * Tseng and Yun (2009) P. Tseng and S. Yun.  Block-coordinate gradient descent method for linearly constrained nonsmooth separable optimization.  _J. Optim. Theory Appl._ , 140(3):513, 2009.
  * Vaiter et al. (2018) S. Vaiter, G. Peyré, and J. Fadili.  Model consistency of partly smooth regularizers.  _IEEE Trans. Inf. Theory_ , 64(3):1725–1737, 2018.
  * Virtanen et al. (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, et al.  Scipy 1.0: fundamental algorithms for scientific computing in python.  _Nature methods_ , 17(3):261–272, 2020.
  * Wengert (1964) R. E. Wengert.  A simple automatic derivative evaluation program.  _Communications of the ACM_ , 7(8):463–464, 1964\.
  * Winston and Kolter (2020) E. Winston and Z. Kolter.  Neural monotone operator equilibrium networks.  _NeurIPS_ , 2020.
  * Zhang et al. (2013) L. Zhang, M. Mahdavi, and R. Jin.  Linear convergence with condition number independent access of full gradients.  _NeurIPS_ , 26:980–988, 2013.
  * Zou and Hastie (2005) H. Zou and T. J. Hastie.  Regularization and variable selection via the elastic net.  _J. R. Stat. Soc. Ser. B Stat. Methodol._ , 67(2):301–320, 2005.

##  A Additional lemmas

###  A.1 Differentiability of the proximal operator

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

######  Lemma 18 (Klopfenstein et al. 2020, Lemmas 2 and 3)

Let $0<\gamma_{j}\leq 1/L_{j}$. Let $\lambda\in\mathbb{R}^{r}$ and $\Lambda$ a
neighborhood of $\lambda$. Consider a solution
$\hat{\beta}\in\operatorname*{\mathrm{arg\,min}}_{\beta\in\mathbb{R}^{p}}\Phi(\beta,\lambda)$
and $\hat{S}$ its generalized support. Suppose

  1. 1.

Assumptions 2, 3 and 6 hold.

  2. 2.

Assumption 4 hold on $\Lambda$.

Then, for all $j\in\hat{S}$, the map
$\beta\mapsto\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}$ is
differentiable at $\hat{\beta}_{\hat{S}}$. Moreover, for all
$j\in\hat{S}^{c}$, $\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}$ is
constant around $\hat{\beta}_{j}-\gamma_{j}\nabla_{j}f(\hat{\beta})$. Thus,
$\beta\mapsto\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(\beta_{j}-\gamma_{j}\nabla_{j}f(\beta))$
is differentiable at $\hat{\beta}$ with gradient $0$.

###  A.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 19

Let $A\in\mathbb{R}^{p\times p},b\in\mathbb{R}$ with $\rho(A)<1$. Let
$(\mathcal{J}^{(k)})_{k\in\mathbb{N}}$ be a sequence of $\mathbb{R}^{p}$ such
that:

$\mathcal{J}^{(k+1)}=A\mathcal{J}^{(k)}+b+\epsilon^{(k)}\enspace,$ (19)

with $(\epsilon^{(k)})_{k\in\mathbb{N}}$ a sequence which converges linearly
to $0$, then $(\mathcal{J}^{(k)})_{k\in\mathbb{N}}$ converges linearly to its
limit $\hat{\mathcal{J}}\triangleq(\operatorname{Id}-A)^{-1}b$.

Proof  Assume $(\epsilon^{(k)})_{k\in\mathbb{N}}$ converges linearly. Then,
there exists $c_{1}>0,0<\nu<1$ such that:

$\lVert\epsilon^{(k)}\rVert\leq c_{1}\nu^{k}\enspace.$

Applying a standard result on spectral norms (see (Polyak, 1987, Chapter 2,
Lemma 1)) yields a bound on $\lVert A^{k}\rVert_{2}$. More precisely, for
every $\delta>0$ there is a constant $c_{2}(\delta)=c_{2}$ such that

$\displaystyle\lVert A^{k}\rVert_{2}\leq c_{2}(\rho(A)+\delta)^{k}\enspace.$

Without loss of generality, we consider from now on a choice of $\delta$ such
that $\rho(A)+\delta<1$. Since $\hat{\mathcal{J}}=(\operatorname{Id}-A)^{-1}b$
the limit $\hat{\mathcal{J}}$ of the sequence satisfies:

$\hat{\mathcal{J}}=A\hat{\mathcal{J}}+b\enspace.$ (20)

Taking the difference between Equations 19 and 20 yields:

$\mathcal{J}^{(k+1)}-\hat{\mathcal{J}}=A(\mathcal{J}^{(k)}-\hat{\mathcal{J}})+\epsilon^{(k)}\enspace.$
(21)

Unrolling Equation 21 yields
$\mathcal{J}^{(k+1)}-\hat{\mathcal{J}}=A^{k+1}(\mathcal{J}^{(0)}-\mathcal{J})+\sum_{k^{\prime}=0}^{k}A^{k^{\prime}}\epsilon^{(k-k^{\prime})}$.
Taking the norm on both sides and using the triangle inequality leads to

$\displaystyle\lVert\mathcal{J}^{(k+1)}-\hat{\mathcal{J}}\rVert_{2}\leq$
$\displaystyle\leavevmode\nobreak\ \lVert
A^{k+1}(\mathcal{J}^{(0)}-\mathcal{J})\rVert_{2}+\sum_{k^{\prime}=0}^{k}\lVert
A^{k^{\prime}}\rVert_{2}\lVert\epsilon^{(k-k^{\prime})}\rVert$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ \lVert
A^{k+1}\rVert_{2}\cdot\lVert\mathcal{J}^{(0)}-\hat{\mathcal{J}}\rVert_{2}+c_{1}\sum_{k^{\prime}=0}^{k}\lVert
A^{k^{\prime}}\rVert_{2}\cdot\nu^{k-k^{\prime}}$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\
c_{2}(\rho(A)+\delta)^{k+1}\cdot\lVert\mathcal{J}^{(0)}-\hat{\mathcal{J}}\rVert_{2}+c_{1}\sum_{k^{\prime}=0}^{k}c_{2}(\rho(A)+\delta)^{k^{\prime}}\nu^{k-k^{\prime}}$

We can now split the last summand in two parts and obtain the following bound,
reminding that $\rho(A)+\delta<1$:

$\displaystyle\lVert\mathcal{J}^{(k+1)}-\hat{\mathcal{J}}\rVert_{2}\leq$
$\displaystyle\leavevmode\nobreak\
c_{2}(\rho(A)+\delta)^{k+1}\cdot\lVert\mathcal{J}^{(0)}-\hat{\mathcal{J}}\rVert_{2}$
$\displaystyle+c_{1}c_{2}\left(\sum_{k^{\prime}=0}^{k/2}(\rho(A)+\delta)^{k^{\prime}}\nu^{k-k^{\prime}}+\sum_{k^{\prime}=k/2}^{k}(\rho(A)+\delta)^{k^{\prime}}\nu^{k-k^{\prime}}\right)$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
c_{2}(\rho(A)+\delta)^{k+1}\cdot\lVert\mathcal{J}^{(0)}-\hat{\mathcal{J}}\rVert_{2}+\frac{c_{1}c_{2}(\rho(A)+\delta)}{1-\rho(A)-\delta}\sqrt{\nu}^{k}$
$\displaystyle+\frac{c_{1}c_{2}\nu}{1-\nu}\sqrt{(\rho(A)+\delta)}^{k}\enspace.$

Thus, $(\mathcal{J}^{(k)})_{k\in\mathbb{N}}$ converges linearly towards its
limit $\hat{\mathcal{J}}$.

##  B Proof of Theorem 12

See 12

Proof  We first prove Theorem 12 for proximal gradient descent.

Proximal gradient descent case. Solving Equation 1 with proximal gradient
descent leads to the following updates:

$\displaystyle\beta^{(k+1)}=\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\underbrace{\beta^{(k)}-\gamma\nabla
f(\beta^{(k)})}_{z^{(k)}})\enspace.$ (22)

Consider the following sequence $(\mathcal{J}^{(k)})_{k\in\mathbb{N}}$ defined
by:

$\displaystyle\mathcal{J}^{(k+1)}=\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z^{(k)})\odot\left(\operatorname{Id}-\gamma\nabla^{2}f(\beta^{(k)})\right)\mathcal{J}^{(k)}+\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z^{(k)})\enspace.$ (23)

Note that if $\operatorname{prox}_{\gamma g(\cdot,\lambda)}$ is not
differentiable with respect to the first variable at $z^{(k)}$ (respectively
with respect to the second variable $\lambda$), any weak Jacobian can be used.
When (H3) holds, differentiating Equation 22 with respect to $\lambda$ yields
exactly Equation 23.

Assumptions 2, 3, 4 and 6 and the convergence of $(\beta^{(k)})$ toward
$\hat{\beta}$ 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, 3, 6 and 4 ensure (see LABEL:{lemma:diff_prox})
$\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z^{k})_{\hat{S}^{c}}=0$, for all $k\geq K$. Thus, for all
$k\geq K$,

$\displaystyle\mathcal{J}^{(k)}_{\hat{S}^{c}:}=\hat{\mathcal{J}}_{\hat{S}^{c}:}=\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z^{(k)})_{\hat{S}^{c}:}\enspace.$

The updates of the Jacobian then become:

$\displaystyle\mathcal{J}^{(k+1)}_{\hat{S}:}=\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z^{(k)})_{\hat{S}}\odot\left(\operatorname{Id}-\gamma\nabla^{2}_{\hat{S},\hat{S}}f(\beta^{(k)})\right)\mathcal{J}_{\hat{S}:}^{(k)}+\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(z^{(k)})_{\hat{S}:}\enspace.$

From Assumption 6, we have that $f$ is locally $\mathcal{C}^{3}$ at
$\hat{\beta}$, $g(\cdot,\lambda)$ is locally $\mathcal{C}^{2}$ at
$\hat{\beta}$ hence $\operatorname{prox}_{g(\cdot,\lambda)}$ is locally
$\mathcal{C}^{2}$. The function
$\beta\mapsto\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\beta-\gamma\nabla
f(\beta))_{\hat{S}}\odot(\operatorname{Id}-\gamma\nabla^{2}_{\hat{S},\hat{S}}f(\beta))$
is differentiable at $\hat{\beta}$. Using (H4) we have that
$\beta\mapsto\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\beta-\gamma\nabla f(\beta))_{\hat{S}:}$ is also
differentiable at $\hat{\beta}$. Using the Taylor expansion of the previous
functions yields:

$\displaystyle\mathcal{J}^{(k+1)}_{\hat{S}:}=\underbrace{\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})_{\hat{S}}\odot\left(\operatorname{Id}-\gamma\nabla^{2}_{\hat{S},\hat{S}}f(\hat{\beta})\right)}_{A}\mathcal{J}_{\hat{S}:}^{(k)}+\underbrace{\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})_{\hat{S}:}}_{b}+\underbrace{o(\lVert\beta^{(k)}-\hat{\beta}\rVert)}_{\epsilon^{(k)}}\enspace.$
(24)

Thus, for $0<\gamma\leq 1/L$,

$\displaystyle\rho(A)\leq\lVert A\rVert_{2}$
$\displaystyle\leq\underbrace{\lVert\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{z})_{\hat{S}}\rVert}_{\leq 1\text{ (non-
expansiveness)}}\,\cdot\underbrace{\lVert\operatorname{Id}-\gamma\nabla_{\hat{S},\hat{S}}^{2}f(\hat{\beta})\rVert_{2}}_{\quad\quad\quad\quad<1\text{
(\lx@cref{creftypecap~refnum}{ass:restricted_injectivity} and $0<\gamma\leq
1/L$)}}<1\enspace.$ (25)

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<\gamma\leq 1/L$.

Assumptions 2, 3, 4, 6 and 7 and the convergence of $(\beta^{(k)})$ toward
$\hat{\beta}$ ensure $(\beta^{(k)})_{k\in\mathbb{N}}$ converges locally
linearly (Liang et al., 2014, Thm. 3.1). The asymptotic autoregressive
sequence in Equation 24, $\rho(A)<1$, and the local linear convergence of
$(\epsilon^{(k)})_{k\in\mathbb{N}}$, yield our result using Lemma 19.

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:

$\displaystyle\beta_{j}^{(k,j)}=\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}\underbrace{\left(\beta_{j}^{(k,j-1)}-\gamma_{j}\nabla_{j}f(\beta^{(k,j-1)})\right)}_{\triangleq
z_{j}^{(k,j-1)}}\enspace.$ (26)

We consider the following sequence:

$\displaystyle\mathcal{J}^{(k,j)}_{j:}=\leavevmode\nobreak\ $
$\displaystyle\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j}^{(k,j-1)})\left(\mathcal{J}^{(k,j-1)}_{j:}-\gamma_{j}\nabla^{2}_{j:}f(\beta^{(k,j-1)})\mathcal{J}^{(k,j-1)}\right)$
$\displaystyle+\partial_{\lambda}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j}^{(k,j-1)})\enspace.$
(27)

Note that if $\operatorname{prox}_{\gamma g(\cdot,\lambda)}$ is not
differentiable with respect to the first variable at $z^{(k)}$ (respectively
with respect to the second variable $\lambda$), any weak Jacobian can be used.
When (H3) holds, differentiating Equation 26 with respect to $\lambda$ yields
exactly Section B.

Assumptions 3, 2, 4 and 6 and the convergence of
$(\beta^{(k)})_{k\in\mathbb{N}}$ toward $\hat{\beta}$ 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 $\hat{S}$ (of cardinality $\hat{s}$) has been
identified, we have that for all $k\geq K$,
$\beta_{\hat{S}^{c}}^{(k)}=\hat{\beta}_{\hat{S}}$ and for any
$j\in\hat{S}^{c}$,
$\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j}^{(k,j-1)})=0$.
Thus
$\mathcal{J}_{j:}^{(k,j)}=\partial_{\lambda}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j}^{(k,j-1)})$.
Then, we have that for any $j\in\hat{S}$ and for all $k\geq K$:

$\displaystyle\mathcal{J}^{(k,j)}_{j:}$
$\displaystyle=\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j}^{(k,j-1)})\left(\mathcal{J}^{(k,j-1)}_{j:}-\gamma_{j}\nabla^{2}_{j,\hat{S}}f(\beta^{(k,j-1)})\mathcal{J}^{(k,j-1)}_{\hat{S}:}\right)$
$\displaystyle\hskip
10.00002pt+\partial_{\lambda}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j}^{(k,j-1)})-\gamma_{j}\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}(z_{j}^{(k,j-1)})\nabla^{2}_{j,\hat{S}^{c}}f(\beta^{(k,j-1)})\mathcal{J}^{(k,j-1)}_{\hat{S}^{c}:}\enspace.$

Let $e_{1},\dots,e_{\hat{s}}$ be the vectors of the canonical basis of
$\mathbb{R}^{\hat{s}}$. We can consider the applications

$\displaystyle\mathbb{R}^{p}$ $\displaystyle\rightarrow\mathbb{R}^{\hat{s}}$
$\displaystyle\beta$
$\displaystyle\mapsto\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}\left(\beta_{j}-\gamma_{j}\nabla_{j}f(\beta)\right)\left(e_{j}-\gamma_{j}\nabla^{2}_{j,\hat{S}}f(\beta)\right)\enspace,$

and

$\displaystyle\mathbb{R}^{p}$
$\displaystyle\rightarrow\mathbb{R}^{\hat{s}\times r}$ $\displaystyle\beta$
$\displaystyle\mapsto\partial_{\lambda}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}\left(\beta_{j}-\gamma_{j}\nabla_{j}f(\beta)\right)-\gamma_{j}\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}\left(\beta_{j}-\gamma_{j}\nabla_{j}f(\beta)\right)\nabla^{2}_{j,\hat{S}^{c}}f(\beta)\hat{\mathcal{J}}_{\hat{S}^{c}:}\enspace,$

which are both differentiable at $\hat{\beta}$ using Assumption 6 and (H4).
The Taylor expansion of the previous functions yields:

$\displaystyle\mathcal{J}^{(k,j)}_{j:}$
$\displaystyle=\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}\left(\hat{z}_{j}\right)\left(e_{j}-\gamma_{j}\nabla^{2}_{j,\hat{S}}f(\hat{\beta})\right)\mathcal{J}^{(k,j-1)}_{\hat{S}:}$
$\displaystyle\hskip
10.00002pt+\partial_{\lambda}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}\left(\hat{z}_{j}\right)-\gamma_{j}\partial_{z}\operatorname{prox}_{\gamma_{j}g_{j}(\cdot,\lambda)}\left(\hat{z}_{j}\right)\nabla^{2}_{j,\hat{S}^{c}}f(\hat{\beta})\mathcal{J}^{(k,j-1)}_{\hat{S}^{c}:}$
$\displaystyle\hskip 10.00002pt+o(||\beta^{(k,j-1)}-\hat{\beta}||)\enspace.$

Let $j_{1},\dots,j_{\hat{s}}$ be the indices of the generalized support of
$\hat{\beta}$. When considering a full epoch of coordinate descent, the
Jacobian is obtained as the product of matrices of the form

$\displaystyle
A_{s}^{\top}=\left(\begin{array}[]{c|c|c|c|c|c|c}e_{1}&\ldots&e_{s-1}&v_{j_{s}}&e_{s+1}&\ldots&e_{\hat{s}}\end{array}\right)\in\mathbb{R}^{\hat{s}\times\hat{s}}\enspace,$
(29)

where
$v_{j_{s}}=\partial_{z}\operatorname{prox}_{\gamma_{j_{s}}g_{j_{s}}}\left(\hat{z}_{j_{s}}\right)\left(e_{s}-\gamma_{j_{s}}\nabla^{2}_{j_{s},\hat{S}}f(\hat{\beta})\right)\in\mathbb{R}^{\hat{s}}.$
A full epoch can then be written

$\displaystyle\mathcal{J}^{(k+1)}_{\hat{S}:}=\underbrace{A_{\hat{s}}A_{\hat{s}-1}\ldots
A_{1}}_{A}\mathcal{J}^{(k)}_{\hat{S}:}+b+\epsilon^{(k)}\enspace,$

for a certain $b\in\mathbb{R}^{\hat{s}}$.

The spectral radius of $A$ is strictly bounded by $1$ (Klopfenstein et al.,
2020, Lemma 8): $\rho(A)<1$. Assumptions 3, 2, 4 and 6 and the convergence of
$(\beta^{(k)})_{k\in\mathbb{N}}$ toward $\hat{\beta}$ ensure local linear
convergence of $(\beta^{(k)})_{k\in\mathbb{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$:

$\mathcal{J}_{\hat{S}:}^{(k+1)}=A\mathcal{J}_{\hat{S}:}^{(k)}+b+\epsilon^{(k)}\enspace,$
(30)

with $(\epsilon^{(k)})_{k\in\mathbb{N}}$ converging linearly to 0.

Recalling $\rho(A)<1$, Lemma 19 and the last display yield our result using.

###  B.1 Proof of Theorem 13 (approximate hypergradients)

See 13

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
$\nabla\mathcal{L}(\lambda)$. Following the analysis of Pedregosa (2016), two
sources of approximation errors arise when computing the hypergradient:

  * •

Approximation errors from the inexact computation of $\hat{\beta}$. Dropping
the dependency with respect to $\lambda$, we denote $\beta$ the approximate
solution and suppose the problem is solved to precision $\epsilon$ with
support identification (H7):

$\displaystyle\begin{cases}&\beta_{\hat{S}^{c}}=\hat{\beta}_{\hat{S}^{c}}\\\
&\lVert\beta_{\hat{S}}-\hat{\beta}_{\hat{S}}\rVert\leq\epsilon\enspace.\end{cases}$

  * •

Approximation errors from the approximate resolution of the linear system,
using (H7) yields:

$\displaystyle\lVert
A^{-1\top}\nabla_{\hat{S}}\mathcal{C}(\beta)-v\rVert\leq\epsilon\enspace.$

The exact solution of the exact linear system $\hat{v}$ satisfies:

$\displaystyle\hat{v}=\hat{A}^{-1\top}\nabla_{\hat{S}}\mathcal{C}(\hat{\beta})\enspace,$

with

$\displaystyle A$
$\displaystyle\triangleq\operatorname{Id}_{|\hat{S}|}-\underbrace{\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(\beta-\gamma\nabla
f(\beta)\right)_{\hat{S}}}_{\triangleq
C}\underbrace{\left(\operatorname{Id}_{|\hat{S}|}-\gamma\nabla^{2}_{\hat{S},\hat{S}}f(\beta)\right)}_{\triangleq
D}\enspace,$ $\displaystyle\hat{A}$
$\displaystyle\triangleq\operatorname{Id}_{|\hat{S}|}-\underbrace{\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(\hat{\beta}-\gamma\nabla
f(\hat{\beta})\right)_{\hat{S}}}_{\triangleq\hat{C}}\underbrace{\left(\operatorname{Id}_{|\hat{S}|}-\gamma\nabla^{2}_{\hat{S},\hat{S}}f(\hat{\beta})\right)}_{\triangleq\hat{D}}\enspace.$

  * •

Using the last two points, the goal is to bound the difference between the
exact hypergradient and the approximate hypergradient,
$\lVert\nabla\mathcal{L}(\lambda)-h\rVert$. Following Algorithm 5, the exact
hypergradient reads

$\displaystyle\nabla\mathcal{L}(\lambda)=\hat{B}\hat{v}+\hat{\mathcal{J}}^{\top}_{\hat{S}^{c}:}\nabla_{\hat{S}^{c}}\mathcal{C}(\hat{\beta})\enspace,$

and similarly for the approximate versions:

$\displaystyle
h=Bv+\mathcal{J}^{\top}_{\hat{S}^{c}:}\nabla_{\hat{S}^{c}}\mathcal{C}(\beta)\enspace,$

with

$\displaystyle B$
$\displaystyle\triangleq\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(\beta-\gamma\nabla
f(\beta)\right)_{\hat{S}:}-\gamma\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(\beta-\gamma\nabla
f(\beta)\right)_{\hat{S}}\odot\left(\nabla^{2}_{\hat{S},\hat{S}^{c}}f(\beta)\right)\hat{\mathcal{J}}_{\hat{S}^{c}:}$
$\displaystyle\hat{B}$
$\displaystyle\triangleq\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(\hat{\beta}-\gamma\nabla
f(\hat{\beta})\right)_{\hat{S}:}-\gamma\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}\left(\hat{\beta}-\gamma\nabla
f(\hat{\beta})\right)_{\hat{S}}\odot\left(\nabla^{2}_{\hat{S},\hat{S}^{c}}f(\hat{\beta})\right)\hat{\mathcal{J}}_{\hat{S}^{c}:}\enspace.$

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

$\displaystyle\lVert\nabla\mathcal{L}(\lambda)-h\rVert$
$\displaystyle=\lVert\hat{B}\hat{v}-Bv+\hat{\mathcal{J}}^{\top}_{\hat{S}^{c}:}\nabla_{\hat{S}^{c}}\mathcal{C}(\hat{\beta})-\hat{\mathcal{J}}^{\top}_{\hat{S}^{c}:}\nabla_{\hat{S}^{c}}\mathcal{C}(\beta)\rVert$
$\displaystyle\leq\lVert\hat{B}\hat{v}-Bv\rVert+\lVert\hat{\mathcal{J}}^{\top}_{\hat{S}^{c}:}\nabla_{\hat{S}^{c}}\mathcal{C}(\hat{\beta})-\hat{\mathcal{J}}^{\top}_{\hat{S}^{c}:}\nabla_{\hat{S}^{c}}\mathcal{C}(\beta)\rVert$
$\displaystyle\leq\lVert\hat{B}\hat{v}-B\hat{v}+B\hat{v}-Bv\rVert+\lVert\hat{\mathcal{J}}^{\top}_{\hat{S}^{c}:}(\nabla_{\hat{S}^{c}}\mathcal{C}(\hat{\beta})-\nabla_{\hat{S}^{c}}\mathcal{C}(\beta))\rVert$
$\displaystyle\leq\lVert\hat{v}\rVert\cdot\lVert\hat{B}-B\rVert+\lVert
B\rVert\cdot\lVert\hat{v}-v\rVert+L_{\mathcal{C}}\lVert\hat{\mathcal{J}}^{\top}_{\hat{S}^{c}:}\rVert\cdot\lVert\beta-\hat{\beta}\rVert\enspace.$
(31)

Bounding $\lVert\hat{v}-v\rVert$ and $\lVert\hat{B}-B\rVert$ in Equation 31
yields the desired result which is bounding the difference between the exact
hypergradient and the approximate hypergradient
$\lVert\nabla\mathcal{L}(\lambda)-h\rVert$.

#### Bound on $\lVert\hat{v}-v\rVert$.

We first prove that $\lVert A-\hat{A}\rVert=\mathcal{O}(\epsilon)$. Let
$L_{H}$ be the Lipschitz constant of the application
$\beta\mapsto\nabla^{2}f(\beta)$, then we have:

$\displaystyle\lVert A-\hat{A}\rVert_{2}$ $\displaystyle=\lVert
CD-\hat{C}\hat{D}\rVert_{2}$ $\displaystyle\leq\lVert
CD-C\hat{D}\rVert_{2}+\lVert C\hat{D}-\hat{C}\hat{D}\rVert_{2}$
$\displaystyle\leq\underbrace{\left\lVert C\right\rVert_{2}}_{\leq 1\text{
(non-expansiveness)}}\underbrace{\lVert D-\hat{D}\rVert_{2}}_{\leq
L_{H}\lVert\beta-\hat{\beta}\rVert\text{ using
\ref{hyp:lip_Hessian}}}+\underbrace{\lVert\hat{D}\rVert_{2}}_{\leq
1}\underbrace{\lVert
C-\hat{C}\rVert_{2}}_{\mathcal{O}(\lVert\beta-\hat{\beta}\rVert)\text{ using
\ref{hyp:prox_lip}}}$ $\displaystyle\leq
L_{H}\lVert\beta-\hat{\beta}\rVert+\mathcal{O}(\lVert\beta-\hat{\beta}\rVert)$
$\displaystyle=\mathcal{O}(\lVert\beta-\hat{\beta}\rVert)\enspace.$ (32)

Let $\tilde{v}$ be the exact solution of the approximate system
$A^{\top}\tilde{v}\triangleq\nabla_{\hat{S}}\mathcal{C}(\beta)$. The following
conditions are met:

  * •

$\hat{v}$ is the exact solution of the exact linear system and $\tilde{v}$ is
the exact solution of the approximate linear system

$\displaystyle\hat{A}^{\top}\hat{v}$
$\displaystyle\triangleq\nabla_{\hat{S}}\mathcal{C}(\hat{\beta})$
$\displaystyle A^{\top}\tilde{v}$
$\displaystyle\triangleq\nabla_{\hat{S}}\mathcal{C}(\beta)\enspace.$

  * •

One can control the difference between the exact matrix in the linear system
$\hat{A}$ and the approximate matrix $A$.

$\displaystyle\lVert
A-\hat{A}\rVert_{2}\leq\delta\lVert\beta-\hat{\beta}\rVert\enspace,$

for a certain $\delta>0$ (Equation 32).

  * •

One can control the difference between the two right-hand side of the linear
systems

$\displaystyle\lVert\nabla_{\hat{S}}\mathcal{C}(\beta)-\nabla_{\hat{S}}\mathcal{C}(\hat{\beta})\rVert\leq
L_{\mathcal{C}}\lVert\beta-\hat{\beta}\rVert\enspace,$

since $\beta\mapsto\nabla\mathcal{C}(\beta)$ is $L_{\mathcal{C}}$-Lipschitz
continuous (H6).

  * •

One can control the product of the perturbations

$\displaystyle\delta\cdot\lVert\beta-\hat{\beta}\rVert\cdot\lVert\hat{A}^{-1}\rVert_{2}\leq\rho<1\enspace.$

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

$\displaystyle\lVert\tilde{v}-\hat{v}\rVert$
$\displaystyle\leq\frac{\epsilon}{1-\epsilon\lVert\hat{A}^{-1}\rVert\delta}\left(L_{\mathcal{C}}\lVert\hat{A}^{-1}\rVert+\lVert\hat{v}\rVert\cdot\lVert\hat{A}^{-1}\rVert\delta\right)$
$\displaystyle\leq\frac{\epsilon}{1-\rho}\left(L_{\mathcal{C}}\lVert\hat{A}^{-1}\rVert+\lVert\hat{v}\rVert\cdot\lVert\hat{A}^{-1}\rVert\delta\right)$
$\displaystyle=\mathcal{O}(\epsilon)\enspace.$ (33)

The bound on $\lVert\tilde{v}-\hat{v}\rVert$ finally yields a bound on the
first quantity in Equation 3, $\lVert v-\hat{v}\rVert$:

$\displaystyle\lVert v-\hat{v}\rVert$ $\displaystyle=\lVert
v-\tilde{v}+\tilde{v}-\hat{v}\rVert$ $\displaystyle\leq\lVert
v-\tilde{v}\rVert+\lVert\tilde{v}-\hat{v}\rVert$ $\displaystyle\leq\lVert
A^{-1}A(v-\tilde{v})\rVert+\lVert\tilde{v}-\hat{v}\rVert$
$\displaystyle\leq\lVert A^{-1}\rVert_{2}\times\underbrace{\left\lVert
A(v-\tilde{v})\right\rVert}_{\leq\epsilon\text{
\ref{hyp:espilon_sol}}}+\underbrace{\lVert\tilde{v}-\hat{v}\rVert}_{\mathcal{O}(\epsilon)\text{
(\lx@cref{creftypecap~refnum}{eq:bound_v_v_tilde})}}$
$\displaystyle=\mathcal{O}(\epsilon)\enspace.$ (34)

#### Bound on $\lVert B-\hat{B}\rVert_{2}$.

We now bound the second quantity in Equation 3 $\lVert B-\hat{B}\rVert_{2}$:

$\displaystyle\lVert B-\hat{B}\rVert_{2}$
$\displaystyle\leq\lVert\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\beta-\gamma\nabla
f(\beta))_{\hat{S}:}-\partial_{\lambda}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{\beta}-\gamma\nabla
f(\hat{\beta}))_{\hat{S}:}\rVert_{2}$
$\displaystyle+\gamma\lVert\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\hat{\beta}-\gamma\nabla
f(\hat{\beta}))_{\hat{S}}\nabla^{2}_{\hat{S},\hat{S}^{c}}f(\hat{\beta})\hat{\mathcal{J}}_{\hat{S}^{c}:}-\partial_{z}\operatorname{prox}_{\gamma
g(\cdot,\lambda)}(\beta-\gamma\nabla
f(\beta))_{\hat{S}}\nabla^{2}_{\hat{S},\hat{S}^{c}}f(\beta)\hat{\mathcal{J}}_{\hat{S}^{c}:}\rVert_{2}$
$\displaystyle\leq L_{1}\lVert\beta-\gamma\nabla
f(\beta)_{\hat{S}:}-\hat{\beta}+\gamma\nabla f(\hat{\beta})\rVert\text{ using
\ref{hyp:prox_lip}}$
$\displaystyle+L_{2}\lVert\hat{\beta}-\beta\rVert\cdot\lVert\hat{\mathcal{J}}_{\hat{S}^{c}:}\rVert\text{
using \ref{hyp:prox_lip} and
\lx@cref{creftypecap~refnum}{ass:taylor_expansion}}$
$\displaystyle=\mathcal{O}(\lVert\hat{\beta}-\beta\rVert)\enspace.$ (35)

Plugging Equations 34 and 35 into Equation 3 yields the desired result:
$\lVert\nabla\mathcal{L}(\lambda)-h\rVert=\mathcal{O}(\epsilon)$.

##  C Additional experiments

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 (Massias et al., 2018, 2020) for the Lasso and Blitz (Johnson and
Guestrin, 2015) for the sparse logistic regression. Table 6 summarizes the
values of the hyperparameters $\lambda$ used in Figures 3, 10 and 11.

Figure 10:  Local linear convergence of the Jacobian for the Lasso. Distance
to optimum for the coefficients $\beta$ (top) and the Jacobian $\mathcal{J}$
(bottom) of the forward-mode differentiation of proximal coordinate descent
(Algorithm 3) on multiple datasets.

Figure 11:  Local linear convergence of the Jacobian for sparse logistic regression. Distance to optimum for the coefficients $\beta$ (top) and the Jacobian $\mathcal{J}$ (bottom) of the forward-mode differentiation of proximal coordinate descent (Algorithm 3) on multiple datasets. Table 6: Dataset characteristics and regularization parameters used in Figures 10, 11 and 3. Datasets | _leukemia_ | _rcv1_ | _news20_ | _real-sim_
---|---|---|---|---
# samples | $n=38$ | $n=$20242$$ | $n=$19996$$ | $n=$72309$$
# features | $p=$7129$$ | $p=$19959$$ | $p=$632982$$ | $p=$20958$$
Lasso | $e^{\lambda}=0.01\,e^{\lambda_{\text{max}}}$ | $e^{\lambda}=0.075\,e^{\lambda_{\text{max}}}$ | $e^{\lambda}=0.3\,e^{\lambda_{\text{max}}}$ | $e^{\lambda}=0.1\,e^{\lambda_{\text{max}}}$
Logistic regression | $e^{\lambda}=0.1\,e^{\lambda_{\text{max}}}$ | $e^{\lambda}=0.25\,e^{\lambda_{\text{max}}}$ | $e^{\lambda}=0.8\,e^{\lambda_{\text{max}}}$ | $e^{\lambda}=0.15\,e^{\lambda_{\text{max}}}$
SVM | $e^{\lambda}=10^{-5}$ | $e^{\lambda}=3\times 10^{-2}$ | $e^{\lambda}=10^{-3}$ | $e^{\lambda}=5\times 10^{-2}$