Smooth Bilevel Programming for Sparse Regularization
Year: 2021

Abstract:
Iteratively reweighted least square (IRLS) is a popular approach to solve sparsityenforcing regression problems in machine learning. State of the art approaches are more efficient but typically rely on specific coordinate pruning schemes. In this work, we show how a surprisingly simple re-parametrization of IRLS, coupled with a bilevel resolution (instead of an alternating scheme) is able to achieve top performances on a wide range of sparsity (such as Lasso, group Lasso and trace norm regularizations), regularization strength (including hard constraints), and design matrices (ranging from correlated designs to differential operators). Similarly to IRLS, our method only involves linear systems resolutions, but in sharp contrast, corresponds to the minimization of a smooth function. Despite being non-convex, we show that there are no spurious minima and that saddle points are “ridable”, so that there always exists a descent direction. We thus advocate for the use of a BFGS quasi-Newton solver, which makes our approach simple, robust and efficient. We perform a numerical benchmark of the convergence speed of our algorithm against state of the art solvers for Lasso, group Lasso, trace norm and linearly constrained problems. These results highlight the versatility of our approach, removing the need to use different solvers depending on the specificity of the ML problem under study.

heading: 1 Introduction
text: Regularized empirical risk minimization is a workhorse of supervised learning, and for a linear model, it reads
min β∈Rn
R(β) + 1
λ L(Xβ, y) (Pλ)
whereX ∈ Rm×n is the design matrix (n being the number of samples andm the number of features), L : Rm × Rm → [0,∞) is the loss function, and R : Rn → [0,∞) the regularizer. Here λ > 0 is the regularisation parameter which is typically tuned by cross-validation, and in the limit case λ = 0, (P0) is a constraint problem minβ R(β) under the constraint L(Xβ, y) = 0. In this work, we focus our attention to sparsity enforcing penalties, which induce some form of structure on the solution of (Pλ), the most celebrated examples (reviewed in Section 2) being the Lasso, group-Lasso and trace norm regularizers. All these regularizers, and much more (as detailed in Section), can be conveniently re-written as an infimum of quadratic functions. While Section 2 reviews more general formulations, this so-called “quadratic variational form” is especially simple in the case of block-separable functionals (such as Lasso and group-Lasso), where one has
R(β) = min η∈Rk+
1
2 ∑ g∈G ||βg||22 ηg + 1 2 h(η), (1)
∗Department of mathematical sciences, University of Bath, Bath BA2 7AY, UK cmshp20@bath.ac.uk †CNRS and DMA, Ecole Normale Supérieure, PSL University, 45 rue d’Ulm, F-75230 PARIS cedex 05,
FRANCE, gabriel.peyre@ens.fr
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
where G is a partition of {1, . . . , n}, k = |G| is the number of groups and h : Rk+ → [0,∞). An important example is the group-Lasso, where R(β) = ∑ g ||βg||2 is a group-`1 norm, in which case
h(η) = ∑ i ηi. The special case of the Lasso, corresponding to the `1 norm is obtained when g = {i} for i = 1, . . . , n and k = n. This quadratic variational form (1) is at the heart of the Iterative Reweighted Least Squares (IRLS) approach, reviewed in Section 1.1. We refer to Section 2 for an in-depth exposition of these formulations.
Sparsity regularized problems (Pλ) are notoriously difficult to solve, especially for small λ, because R is a non-smooth function. It is the non-smoothness of R which forces the solutions of (Pλ) to belong to low-dimensional spaces (or more generally manifolds), the canonical example being spaces of sparse vectors when solving a Lasso problem. We refer to Bach et al. [2011] for an overview of sparsity-enforcing regularization methods. The core idea of our algorithm is that a simple reparameterization of (1) combined with a bi-level programming (i.e. solving two nested optimization problems) can turn (Pλ) into a smooth program which is much better conditioned, and can be tackled using standard but highly efficient optimization techniques such as quasi-Newton (L-BFGS). Indeed, by doing the change of variable (vg, ug) , ( √ ηg, βg/ √ ηg) in (1), (Pλ) is equivalent to
min v∈Rk f(v) where f(v) , min u∈Rn G(u, v) (2)
G(u, v) , 1 2 h(v v) + 1 2 ||u||2 + 1 λ L(X(v G u), y). (3)
Throughout, we define to be the standard Hadamard product and for v ∈ Rk and u ∈ Rn, we define v G u ∈ Rn to be such that (v G u)g = vgug. Provided that v 7→ h(v v) is differentiable and L(·, y) is a convex, proper, lower semicontinuous function, the inner minimisation problem has a unique solution and f is differentiable. Moreover, in the case of the quadratic loss L(z, y) , 12 ||z − y|| 2 2, the gradient of f can be computed in closed form, by solving a linear system of dimension m or n. This paper is thus devoted to study the theoretical and algorithmic implications of this simple twist on the celebrated IRLS approach.
Comparison with proximal gradient To provide some intuition about the proposed approach, the figure on the right contrasts the iterations of a gradient descent on f and of the iterative soft thresholding algorithm (ISTA) on the Lasso. We consider a random Gaussian matrix X ∈ R10×20 with λ = ||X>y||∞/10. The ISTA trajectory is non-smooth when some feature crosses 0. In particular, if a coefficient (such as the red one on the figure) is initialized with the wrong sign, it takes many iteration for ISTA to flip sign. In sharp contrast, the gradient flow of f does not exhibit such a singularity and exhibits a smooth geometric convergence. We refer to the appendix for an analysis of this phenomenon.
Contributions Our main contribution is a new versatile algorithm for sparse regularization, which applies standard smooth optimization methods to minimize the function f in (2). We first propose in Section 2 a generic class of regularizers R that enjoy a quadratic variational form. This section recaps existing results under a common umbrella and shows the generality of our approach. Section 3.1 then gathers the theoretical analysis of the method, and in particular the proof that while being non-convex, the function f has no local minimum and only “ridable” saddle points. As a result, one can guarantee convergence to a global minimum for many optimisation schemes, such as gradient descent with random perturbations [Lee et al., 2017, Jin et al., 2017] or trust region type methods [Pascanu et al., 2014]. Furthermore, for the case of the group Lasso, we show that f is an infinitely differentiable function with uniformly bounded Hessian. Consequently, standard solvers such as Newton’s method/BFGS can be applied and with a superlinear convergence guarantee. Section 4 performs a detailed numerical study of the method and benchmarks it against several popular competing algorithms for Lasso, group-Lasso and trace norm regularization. Our method is consistently amongst the best performers, and is in particular very efficient for small values of λ, and can even cope with the constrained case λ = 0.

heading: 1.1 Related works
text: State of the art solvers for sparse optimisation Popular approaches to nonsmooth sparse optimisation include proximal based methods. The simplest instance is the Forward-Backward algorithm
[Lions and Mercier, 1979, Daubechies et al., 2004] which handles the case where L is a smooth function. There are many related inertial based acceleration schemes, such as FISTA [Beck and Teboulle, 2009] and particularly effective techniques that leads to substantial speedups are the adaptive use of stepsizes and restarting strategies [O’donoghue and Candes, 2015]. Other related approaches are proximal quasi-Newton and variable metric methods [Combettes and Vũ, 2014, Becker et al., 2019], which incorporate variable metrics into the quadratic term of the proximal operator. Another popular approach, particularly for the Lasso-like problems are coordinate descent schemes [Friedman et al., 2010]. These schemes are typically combined with support pruning schemes [Ghaoui et al., 2010, Ndiaye et al., 2017, Massias et al., 2018]. In the case where both the regularisation and loss terms are nonsmooth (e.g. in the basis pursuit setting), typical solvers are the primal-dual [Chambolle and Pock, 2011], ADMM [Boyd et al., 2011] and Douglas-Rachford algorithms [Douglas and Rachford, 1956]. Although these schemes are very popular due to their relatively low per iteration complexity, these methods have sublinear convergence rates in general, with linear convergence under strong convexity.
The quadratic variational formulation and IRLS Quadratic variational formulations such as (1) have been exploited in many early computer vision works, such as Geiger and Yuille [1991] and Geman and Reynolds [1992]. One of the first theoretical results can be found in Geman and Reynolds [1992], see also Black and Rangarajan [1996] which lists many examples. Our result in Theorem 1 can be seen as a generalisation of these early works. Norm regularizers of this form were introduced in Micchelli et al. [2013], and further studied in the monograph Bach et al. [2011] under the name of subquadratic norms. The main algorithmic consequence of the quadratic variational formulation in the literature is the iterative reweighted least squares (IRLS) algorithm. When L(z, y) = 12 ||z − y|| 2 2 is the quadratic loss, a natural optimisation strategy is alternating minimising. However, due to the 1/ηg term, one needs to introduce some regularisation to ensure convergence. One popular approach is to add ε2 ∑ g η −1 g to the formulation (1) which leads to the IRLS algorithm [Daubechies et al., 2010]. The ε-term can be seen as a barrier to keep ηg positive which is reminiscent of interior point methods. The idea of IRLS is to do alternating minimisation over β and η, where the minimisation with respect to β is a least squares problem, and the minimisation with respect to η admits a closed form solution. A nuclear norm version of IRLS has been used in Argyriou et al. [2008] where an alternating minimisation algorithm was introduced. Finally, we remark that although nonconvex formulations for various low complexity regularizers have appeared in the literature, see for instance Rennie and Srebro [2005], Hastie et al. [2015], Mardani and Giannakis [2015], Hoff [2017] for the case of the `1 and nuclear norms, they are typically associated with alternating minimisation algorithms.
Variable projection/reduced gradient approaches IRLS methods are quite slow because the resulting minimization problem is poorly conditioned. Adding the smoothing term ε2 ∑ g η −1 g only partly alleviates this, and also breaks the sparsity enforcing property of the regularizer R. We avoid both issues in (2) by solving a “reduced” problem which is much better conditioned and smooth. This idea of solving a bi-variate problem by re-casting it as a bilevel program is classical, we refer in particular to [Rockafellar and Wets, 2009, Chap. 10] for some general theoretical results on reduced gradients, and also Danskin’s theorem (although this is restricted to the convex setting) in Bertsekas [1997]. Our formulation falls directly into the framework of variable projection [Ruhe and Wedin, 1980, Golub and Pereyra, 2003], introduced initially for solving nonlinear least squares problems. Properties and advantages of variable projection have been studied in Ruhe and Wedin [1980], we refer also to Hong et al. [2017], Zach and Bourmaud [2018] for more recent studies. Nonsmooth variable projection is studied in van Leeuwen and Aravkin [2016], although the present work is in the classical setting of variable projection due to our smooth reparametrization. Reduced gradients have also been associated with the quadratic variational formulation in several works [Bach et al., 2011, Pong et al., 2010, Rakotomamonjy et al., 2008]. The idea is to apply descent methods over g(η) = minβ R0(η, β) + 1 2 ||Xβ − y|| 2 2. Although the function over η and β is discontinuous, the function g over η is smooth and one can apply first order methods, such as proximal gradient descent to minimise g under positivity constraints. While quasi-Newton methods can be applied in this setting with bound constraints, we show in Section 4.1 that this approach is typically less effective than our nonconvex bilevel approach. In the setting of the trace norm, the optimisation problem is constrained on the set of positive semidefinite matrices, so one is restricted to using first order methods [Pong et al., 2010].

heading: 2 Quadratic variational formulations
text: We describe in this section some general results about when a regulariser has a quadratic variational form. Our first result brings together results which are scattered in the literature: it is closely related to Theorem 1 in Geman and Reynolds [1992], but their proof was only for strictly concave differentiable functions and did not explicitly connect to convex conjugates, while the setting for norms have been characterized in the monograph Bach et al. [2011] under the name of subquadratic norms.
Theorem 1. Let R : Rn → R. The following are equivalent: (i) R(β) = ϕ(β β) where ϕ is proper, concave and upper semi-continuous, with domain Rd+. (ii) There exists a convex function ψ for which R(β) = infz∈Rn+ 1 2 ∑n i=1 ziβ 2 i + ψ(z). Furthermore, ψ(z) = (−ϕ)∗(−z/2) is defined via the convex conjugate (−ϕ)∗ of −ϕ, leading to (1) using the change of variable η ← 1/z and h(η) = 2ψ(1/η). When R is a norm, the function h can be written in terms of the dual norm R∗ as h(η) = maxR∗(w)61 ∑ i w 2 i ηi. Moreover,
R(β)2 = infη∈Rn+ {∑ i β 2 i η −1 i \ h(η) 6 1 } .
See Appendix B for the proof to this Theorem. Some additional properties of ψ are derived in Lemma 1 of Appendix B, one property is that if R is coercive, then lim||z||→0 ψ(z) = +∞, so the function f is coercive (see also remark 2).

heading: 2.1 Examples
text: Let us first give some simple examples of both convex and non-convex norms: – Euclidean norms: for R = || · ||2, making use of R∗ as stated in Theorem 1, one has h = || · ||∞. – Group norms: for G is a partition of {1, . . . , n}, the group norm is R(β) = ∑ g∈G ||βg||. Using
the previous result for the Euclidean norm, one has h(η) = ∑ g∈G ||(ηi)i∈g||∞. This expression can be further simplified to obtain (1) with a reduced vector η in R|G|+ in place of Rn by noticing that the optimal η is constant in each group g.
– `q (quasi) norms: For R(β) = |β|q where q ∈ (0, 2), one has ϕ(u) = uq/2 and one verifies that h(η) = Cqη q 2−q where Cq = (2− q)qq/(2−q). Note that for q > 23 , v 7→ h(v
2) = vγ for γ > 1 is differentiable. Analysis and numerics for this nonconvex setting can be found in Appendix F.
Matrix regularizer The extension of Theorem 1 to the case where β = B is a matrix can be found in the appendix. When R = ϕ(BB>) is a function on matrices, the analogous quadratic variational formulation is
R(B) = min Z∈Sn+ min B∈Rn×r
1
2 ∑ g tr(B>Z−1B) + 1 2 h(Z), (4)
where Sn+ denotes the set of symmetric positive semidefinite matrices and h(Z) = 2(−ϕ)∗(−Z−1/2). Letting U = Z−1/2B and V = Z1/2, we have B = V U and the equivalence
min B
R(B) + 1
2λ ||A(B)− y||22 (5)
= min V ∈Rn×n f(V ) , min U∈Rn×r
1 2 ||U ||2F + 1 2 h(V >V ) + 1 2λ ||A(V U)− y||22. (6)
where A : Rn×r → Rm is a linear operator. Again, provided that V 7→ h(V >V ) is differentiable, f is a differentiable function with ∇f(V ) = V ∂h(V >V ) + 1λA
∗(A(V U)− y))U> and U such that λU + V >A∗(A(V U)− y) = 0. For the trace norm, R(B) = tr( √ B>B), we have h(Z) = tr(Z) and ∂h(Z) = Id. Note that, just in the vectorial case, one could write the inner minimisation problem over the dual variable α ∈ Rm and handle the case of λ = 0.

heading: 3 Theoretical analysis
text: Our first result shows the equivalence between (Pλ) and a smooth bilevel problem.
Theorem 2. Denote Ly , L(·, y) and let L∗y denote the convex conjugate of Ly. Assume that Ly is a convex, proper, lower semicontinuous function and R takes the form (1). The problem (Pλ) is equivalent to
min v∈Rk f(v) , min u∈Rn
1 2 h(v v) + 1 2 ||u||2 + 1 λ Ly (X(v G u)) . (7)
= max α∈Rm
1 2 h(v v)− 1 λ L∗y (λα)− 1 2 ||v G (X>α)||22. (8)
where v G u , (ugvg)g∈G . The minimiser β to (Pλ) and the minimiser v to (7) are related by β = v G u = −v2 G X>α. Provided that v 7→ h(v v) is differentiable, the function f is differentiable with gradient
∇f(v) = v ∂h(v2)−v G ( ||X>g α||2 ) g where α ∈ argmaxα̃−L∗y(α̃)− 1 2 ||v G(X>α̃)||22. (9) Note that ∇f is uniquely defined even if α is not unique. If λ > 0 and Ly is differentiable, then we have the additional formula, with u ∈ argminũ 12 ||ũ|| 2 2 + 1 λLy(X(v G ũ)),
∇f(v) = v ∂h(v2) + 1 λ
( 〈ug, X>g ∂Ly(X(v G u)〉 ) g∈G . (10)
Proof. The equivalence between (Pλ) and (7) is simply due to the quadratic variational form of R, and the change of variable vg = √ ηg , and ug = βg/ √ ηg . The equivalence to (8) follows by convex duality on the inner minimisation problem, that is
f(v) = min u G1(v, u) ,
1 2 h(v2) + 1 2 ||u||2 + 1 λ Ly(X(v G u))
= min u,z
1 2 h(v2) + 1 2 ||u||2 + 1 λ Ly(z) where z = X(v G u)
= min u,z max α
1 2 h(v2) + 1 2 ||u||2 + 1 λ Ly(z)− 〈α, z〉+ 〈α, X(v G u)〉
= max α min u
1 2 h(v2) + 1 2 ||u||2 + 〈α, X(v G u)〉 − 1 λ L∗y(λα).
Using the optimality condition over u, we obtain u = −v G X>α and hence,
f(v) = max α
G2(v, α) , 1 2 h(v2)− 1 2 ||v G X>α||22 − 1 λ L∗y(λα).
By [Rockafellar and Wets, 2009, Theorem 10.58], if the following set S(v) , { ∂vG2(v, α) = v − v G (||X>g α||22)g∈G \ α ∈ argminαG2(v, α) } is singled valued, then f is differentiable with ∇f(v) = ∂vG2(v, α) for α ∈ argminαG2(v, α). Observe that even if argminαG2(v, α) is not single valued, sinceG2 is strongly convex for v GX>α, S(v) is single-valued and hence, f is a differentiable function. In the case where Ly is differentiable, we can again apply [Rockafellar and Wets, 2009, Theorem 10.58], to obtain ∇f(v) = v ∂h(v2) + ∂vG1(v, y) with u = argminuG1(v, u) (noting that G1(v, ·) is strongly convex and has a unique minimiser) which is precisely the gradient formula (10).
For the Lasso and basis pursuit setting, the gradient of f can be computed in closed form: Corollary 1. If R has a quadratic variational form and Ly(z) = 12 ||z − y|| 2 2, then ∂Ly(z) = z − y, L∗y(α) = 〈y, α〉+ 12 ||α|| 2 2 and the gradient of f can be written as in (9) and additionally (10) when λ > 0. Furthermore, α ∈ Rm in (9) and u ∈ Rn in (10) solves (X diag(v̄2)X> + λId)α = −y, and (diag(v̄)X>X diag(v̄) + λId)u = v G (X>y), (11)
where v̄ ∈ Rn is defined as v̄ u = (vgug)g∈G for all u ∈ Rn and Id denotes the identity matrix.
This shows that our method caters for the case λ = 0 with the same algorithm in a seamless manner. This is unlike most existing approach which work well for λ > 0 (and typically do not require matrix inversion) but fails when λ is small, whereas solvers dedicated for λ = 0 might require inverting a linear system, see Section 4.4 for an illustrative example.

heading: 3.1 Properties of the projected function f
text: In this section, we analyse the case of the group Lasso. The following theorem ensures that the projected function f has only strict saddle points or global minima. We say that v is a second order stationary point if ∇f(v) = 0 and ∇2f(v) 0. We say that v is a strict saddle point (often called “ridable”) if it is a stationary point but not a second order stationary point. One can thus always find a direction of descent outside the set of global minimum. This can be exploited to derive convergence guarantees to second order stationary points for trust region methods [Pascanu et al., 2014] and gradient descent methods [Lee et al., 2017, Jin et al., 2017]. Theorem 3. In the case h(z) = ∑ i zi and L(z, y) = 1 2 ||z−y|| 2 2, the projected function f is infinitely
continuously differentiable and for v ∈ Rk,∇f(v) = v ( 1− |ξ|2 ) where ξg = 1λX > g (X(u v)−y) and u solves the inner least squares problem for v. Let J denote the support of v, by rearranging the columns and rows, the Hessian of f can be written as the following block diagonal matrix
∇2f(v) = (
diag(1− ||ξg||22)g∈J + 4U>WU 0 0 diag(1− ||ξg||22)g∈Jc
) (12)
where W , Id− λ ( (vgX > g Xhvh)g,h∈J + λIdJ )−1 and U is the block diagonal matrix with blocks (ξg)g∈J , with maxg∈G ||ξg||2 6 C and ||∇2f(v)|| 6 1 + 3C2 where C , ||y||2 maxg∈G ||Xg||/λ. Moreover, all stationary points of f are either global minima or strict saddles. At stationary points, the eigenvalues of the Hessian of f are at most 4 and is at least
min ( 4(1− λ/(λ+ σ̂)),min
g 6∈J (1− ||ξg||2) ) where σ̂ is the smallest eigenvalue of (vgX>g Xhvh)g,h∈J .
The proof can be found in Appendix C. We simply mention here that by examining the first order condition of (Pλ), we see that β is a minimizer if and only if ξ satisfies −ξg = βg||βg||2 for all g ∈ Supp(β) and ||ξg||2 6 1 for all g ∈ G. The first condition on the support of β is always satisfied at stationary points of the nonconvex function (2), and by examining (12), the second condition is also satisfied unless the stationary point is strict. Remark 1 (Example of strict saddle point for our f ). One can observe that v = 0 is a strict saddle point, as the solution to the associated linear system yields u = 0 and hence ∇f(v) = 0. If λ > ||X>y||∞, then u = v = 0 corresponds to a global minimum, otherwise, it is clear to see that there exists g such that 1− ||ξg||2 < 0 and v = 0 is a strict saddle point. Remark 2. Since f ∈ C∞, it is Lipschitz smooth on any bounded domain. As mentioned, f is coercive when R is coercive, and hence, its sublevel sets are bounded. So, for any descent algorithm, we can apply results based on ∇kf being Lipschitz smooth for all k. Remark 3. The nondegeneracy condition that ||ξg||2 < 1 outside the support of v and invertibility of (X>g Xh)g,h∈J is often used to derive consistency results [Bach, 2008, Zhao and Yu, 2006]. By Proposition 1, we see that this condition guarantees that the Hessian of f is positive definite at the minimum, and hence, combining with the smoothness properties of f explained in the previous remark, BFGS is guaranteed to converge superlinearly for starting points sufficiently close to the optimum [Nocedal and Wright, 2006, Theorem 6.6].

heading: 4 Numerical experiments
text: In this section, we use L-BFGS [Byrd et al., 1995] to optimise our bilevel function f and we denote the resulting algorithm “Noncvx-Pro”. Throughout, the inner problem is solved exactly using either a full or a sparse Cholesky solver. One observation from our numerics below is that although NoncvxPro is not always the best performing, unlike other solvers, it is robust to a wide range of settings: for example, our solver is mostly unaffected by the choice of λ while one can observe in Figures 1 and 2 that this has a large impact on the proximal based methods and coordinate descent. Moreover, Noncvx-Pro is simple to code and rely on existing robust numerical routines (Cholesky/ conjugate gradient + BFGS) which naturally handle sparse/implicit operators, and we thus inherit their nice convergence properties. All numerics are conducted on 2.4 GHz Quad-Core Intel Core i5 processor with 16GB RAM. The code to reproduce the results of this article is available online3. 3 https://github.com/gpeyre/2021-NonCvxPro

heading: 4.1 Lasso
text: We first consider the Lasso problem where R(β) = ∑n i=1 |βi|.
Datasets. We tested on 8 datasets from the Libsvm repository4. These datasets are mean subtracted and normalised by m.
Solvers. We compare against the following 10 methods: 1. 0-mem SR1: a proximal quasi newton method [Becker et al., 2019]. 2. FISTA w/ BB: FISTA with Barzilai–Borwein stepsize [Barzilai and Borwein, 1988] and restarts
[O’donoghue and Candes, 2015]. 3. SPG/SpaRSA: spectral projected gradient [Wright et al., 2009]. 4. CGIST: an active set method with conjugate gradient iterative shrinkage/thresholding [Goldstein
and Setzer, 2010]. 5. Interior point method: from Koh et al. [2007]. 6. CELER: a coordinate descent method with support pruning [Massias et al., 2018]. 7. Non-cvx-Alternating-min: alternating minimisation of u and v in G from (3) [Hoff, 2017]. 8. Non-cvx-LBFGS: Apply L-BFGS to minimise the function (u, v) 7→ G(u, v) in (3). 9. L-BFGS-B [Byrd et al., 1995]: apply L-BFGS-B under positivity constraints to minu,v∈Rn+ ∑ i ui+∑
i vi+ 1 2λ ||X(u−v)−y|| 2 2. This is the standard approach for applying L-BFGS to `1 minimisation and corresponds to splitting β into its positive and negative parts. 10. Quad-variational: Based on our idea of Noncvx-Pro, another natural (and to our knowledge
novel) approach is to apply L-BFGS-B to the bilevel formulation of (1) without nonconvex reparametrization. Indeed, by applying (1) and using convex duality, the Lasso can solved by minimizing g(η) , maxα∈Rm 12 ∑ i ηi − λ 2 ||α|| 2 − 12 ∑ i ηi|〈xi, α〉|2 + 〈α, y〉. The gradient of g is g(η) = λ2 − 1 λ |X
>α|2 where | · |2 is in a pointwise sense and α maximises the inner problem, and we apply L-BFGS-B with positivity constraints to minimise g.
Experiments. The results are shown in Figure 1 (with further experiments in the appendix). We show comparisons at different regularisation strengths, with λ∗ being the regularisation parameter found by 10 fold cross validation on the mean squared error, and λmax = ||X>y||∞ is the smallest parameter at which the Lasso solution is guaranteed to be trivial.

heading: 4.2 Group Lasso
text: The multi-task Lasso [Gramfort et al., 2012] is the problem (7) where one minimises over β ∈ Rn×q , the observed data is y ∈ Rm×q and R(β) = ∑n j=1 ||βj ||2 with βj ∈ Rq denotes the jth row of the matrix β and the loss function is 12 ||y −Xβ|| 2 F in terms of the the Frobenius norm.
4 https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/
Datasets. We benchmark our proposed scheme on a joint MEG/EEG data from Ndiaye et al. [2015]. X denotes the forward operator with n = 22494 source locations, and is constructed from 301 MEG sensors and 59 EEG sensors, so that m = 360. The observations y ∈ Rm×q represent time series measurements at m sensor with q = 181 timepoints each. The β corresponds to the source locations which are assumed to remain constant across time.
Solvers. We perform a comparison against FISTA with restarts and BB step, the spectral projected gradient method SPG/SpaRSA, and CELER. Since n is much larger than q andm, we use for NonCvxPro the saddle point formulation (8) where the maximisation is over α ∈ Rm×q . Computation of α in ∇f(v) involves solving (Idm + 1λX diag(v 2)X>)α = y., that is, q linear systems each of size m.
Experiments. Figure 2 displays the objective convergence plots against running time for different λ: λ = λmax/r with λmax = maxi ||X>i y||2 for r = 10, 20, 50, 100. We observe substantial performance gains over the benchmarked methods: In MEG/EEG problems, the design matrix tends to exhibit high correlation of columns and proximal-based algorithms tend to perform poorly here. Coordinate descent with pruning is known to perform well here when the regularisation parameter large [Massias et al., 2018], but its performance deteriorates as λ increases.

heading: 4.3 Trace norm
text: In multi-task feature learning [Argyriou et al., 2008], for each task t = 1, . . . , T , we aim to find ft : Rn → R given training examples (xt,i, yt,i) ∈ Rn × R for i = 1, . . . ,mt. One approach is to jointly solve these T regression problems by minimising (5) where R is the trace norm (also called “nuclear norm”), r = T , y = (yt) T t=1, A(B) = (XtBt) T t=1 with Xt being the matrix with ith row as xt,i and Bt being the tth column of the matrix B. Note that letting ut ∈ Rn denote the tth column of U , the computation of U in ∇f(V ) involves solving T linear systems (λIdn + V >X>t XtV )ut = V >X>t yt. Here, the trace norm encourages the tasks to share a small number of linear features.
Datasets. We consider the three datasets commonly considered in previous works. The Schools dataset 5 from Argyriou et al. [2008] consists of the scores of 15362 students from 139 schools. There are therefore T = 139 tasks with 15362 = ∑ tmt data points in total, and the goal is to map n = 27 student attributes to exam performance. The SARCOS dataset 6 [Zhang and Yang, 2021, Zhang and Yeung, 2012] has 7 tasks, each corresponding to learning the dynamics of a SARCOS anthropomorphic robot arms. There are n = 21 features and m = 48, 933 data points, which are shared across all tasks. The Parkinsons dataset [Tsanas et al., 2009] 7 which is made up of m = 5875 datapoints from T = 42 patients. The goal is to map n = 19 biomarkers to Parkinson’s disease symptom scores for each patient.
Solvers. Figure 3 reports a comparison against FISTA with restarts and IRLS. The IRLS algorithm for (5) is introduced in Argyriou et al. [2008] (see also Bach et al. [2011]), and applies alternate minimisation after adding the regularisation term ελ tr(Z−1)/2 to (4). The update of B is a simple least squares problem while the update for Z is Z ← (BB> + εId) 12 . Our nonconvex approach has

heading: 5 https://home.ttic.edu/~argyriou/code/ 6 http://www.gaussianprocess.org/gpml/data/
text: 7 http://archive.ics.uci.edu/ml/datasets/Parkinsons+Telemonitoring
the same per-iteration complexity as IRLS, but one advantage is that we directly deal with (5) without any ε regularisation. Remark 4. Quad-variational mentioned in Section 4.1 does not extend to the trace norm case, since the function g would be minimised over Sn+, for which the application of L-BFGS-B is unclear. For this reason, bilevel formulations for the trace norm [Pong et al., 2010] have been restricted to the use of first order methods.
Experiments. For each dataset, we compute the regularisation parameter λ∗ by 10-fold crossvalidation on the RMSE averaged across 10 random splits. Then, with this regularisation parameter, we compare Non-convex-pro, FISTA and IRLS with different choices of ε. The convergence plots are show in Figure 3. We observe substantial computational gains for the Schools and Parkinson’s dataset. For the SARCOS dataset, IRLS performed the best, and even though Nonconvex-pro is comparatively less effective here, although we remark that the number of tasks is much smaller (T = 7) and the recorded times are much shorter (less than 0.2s). Further numerical illustrations with synthetic data are shown in the appendix – our method is typically less sensitive to problem variations.

heading: 4.4 Constraint Group Lasso and Optimal Transport
text: A salient feature of our method is that it can handle arbitrary small regularization parameter λ and can even cope with the constrained formulation, when λ = 0, which cannot be tackled by most state of the art Lasso solvers. To illustrate this, we consider the computation of an Optimal Transport (OT) map, which has recently gained a lot of attention in ML [Peyré et al., 2019]. We focus on the Monge problem, where the ground cost is the geodesic distance on either a graph or a surface (which extends original Monge’s problem where the cost is the Euclidean distance). This type of OT problems has been used for instance to analyze and process brain signals in M/EEG [Gramfort et al., 2015], for application in computer graphics [Solomon et al., 2014] and is now being applied to genomic datasets [Schiebinger et al., 2019]. As explained for instance in [Santambrogio, 2015, Sec.4.2], the optimal transport between two probability measures a and b on a surface can be computed by advecting the mass along a vector field v(x) ∈ R3 (tangent to the surface) with minimum vectorial L1 norm ∫ ||v(x)||dx (where dx is the surface area measure) subject to the conservation of mass div(v) = a − b. Once discretized on a 3-D mesh, this boils down to solving a constrained group Lasso problem (P0) where βg ∈ R3 is a discretization of v(zg) at some vertex zg of the mesh, X is a finite element discretization of the divergence using finite elements and y = a− b. The same applies on a graph, in which case the vector field is aligned with the edge of the graph and the divergence is the discrete divergence associated to the graph adjacency matrix, see Peyré et al. [2019]. This formulation is often called “Beckmann problem”.
Datasets. We consider two experiments: (i) following Gramfort et al. [2015] on a 3-D mesh of the brain with n = 20000 vertices with localized Gaussian-like distributions a and b (in blue and red), (ii) a 5-nearest neighbors graph in a 7-dimensional space of gene expression (corresponding to 7 different time steps) of baker’s yeast, which is the dataset from DeRisi et al. [1997]. The n = 614 nodes correspond to the most active genes (maximum variance across time) and this results in a graph with 1927 edges. The distributions are synthetic data, where a is a localized source whereas b is more delocalized.
Solvers. We test our method against two popular first order schemes: the Douglas-Rachford (DR) algorithm [Lions and Mercier, 1979, Combettes and Pesquet, 2007] (DR) and the primal-dual scheme of Chambolle and Pock [2011]. DR is used in its dual formulation (the Alternating Direction Method of Multipliers – ADMM) but on a re-parameterized problem in Solomon et al. [2014]. The details of these algorithms are given in Appendix E.
Experiments. Figure 4 shows the solution of this Beckmann problem in these two settings. While DR has the same complexity per iteration as the computation of the gradient of f (resolution of a linear system), a chief advantage of PD is that it only involves the application of X and X> et each iterations. Both DR and PD have stepping size parameters (denoted µ and σ) which have been tuned manually (the rightmost figure shows two other sub-optimal choices of parameters). In contrast, our algorithm has no parameter to tune and is faster than both DR and PD on these two problems.

heading: 5 Conclusion
text: Most existing approaches to sparse regularisation involve careful smoothing of the nonsmooth term, either by proximal operators or explicit regularisation as in IRLS. We propose a different direction: a simple reparameterization leads to a smooth optimisation problem, and allows for the use of standard numerical tools, such as BFGS. Our numerical results demonstrate that this approach is versatile, effective and can handle a wide range of problems.
We end by making some remarks on possible future research directions. The application of our method to other loss functions requires the use of an inexact solver for the inner problems, and controlling the impact of its approximation is an interesting avenue for future work. Furthermore, it is possible that one can obtain further acceleration by combining with screening rules or active set techniques.

heading: Acknowledgments
text: The work of G. Peyré was supported by the French government under management of Agence Nationale de la Recherche as part of the “Investissements d’avenir” program, reference ANR19-P3IA0001 (PRAIRIE 3IA Insti- tute) and by the European Research Council (ERC project NORIA).

Review summaries:
The authors study regularized regression with a sparsity inducing constraint. (In the appendix also the extension to the low-rankness inducing case is considered.) The authors propose a new non-convex formulation, which is optimized via BFGS. The non-convex formulation is motivated by the well-known IRLS (Iteratively Reweighted Least Squares) algorithm. Extensive numerical simulations are conducted on a wide variety of problems and the proposed algorithm is compared with a number of state-of-the-art approaches. While the proposed algorithm does not always perform best, it is in all cases at least competitive to the best performing approach. Moreover, the authors prove that the induced non-convex function, which is optimized in their approach, is ridable -- meaning that all local minima are global minima and the Hessians of all saddle points have at least one negative eigenvalue. 

This paper proposes a reparametrization of common objective functions in the sparse estimation literature (for linear models). This allows rewriting LASSO, group LASSO, low-rank estimation as a smooth bilevel problem, and then applying fast methods for smooth problems such as L-BFGS. Critically, the gradient to the outer problem has a simple form to compute. The smooth problem is shown to have no spurious saddles, although it is nonconvex. The paper then shows impressive empirical results on multiple problems of interest stemming from this simple reparametrization.

The authors propose a bilevel optimization approach for solving regularized empirical risk minimization. The considered minimization problem ($\cal{P}_{\lambda}$) refers to a linear model with parameters $\beta$ that are regularized in terms of a sparsity enforcing term $R(\beta)$. More specifically, the authors treat regularizers $R$ that can be recast in a so-called quadratic variational form. Therewith, a reparametrization of the original problem can be regarded as a bilevel program with smooth inner problem. Subsequent to the introduction, the authors provide a characterization of regularizers that have a quadratic variational form, followed by a theoretical analysis and numerical experiments.

Authors proposed a re-parametrization of the variational form of usually sparsity enforcing penalties (including the Lasso and its variations). The proposed variational re-parametrization leads to a non-convex but smooth bilevel optimization problem, that authors propose to solve using usual L-BFGS. In addition to an extensive benchmark, authors show that the potential saddle points of the new minimized function are not "pathological".