Alain Rakotomamonjy

Gilles Gasso

Joseph Salmon

Screening Rules for Lasso with Non-Convex Sparse Regularizers

  L'archive ouverte pluridisciplinaire

Introduction

Sparsity-inducing penalties are classical tools in statistical machine learning, especially in settings where the data available for learning is scarce and of high-dimension. In addition, when the solution of the learning problem is known to be sparse, using those penalties yield to models that can leverage this prior knowledge. The Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and the Basis pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Chen | Basis pursuit[END_REF] where the first approaches that have employed 1 -norm penalty for inducing sparsity.

While the Lasso has had a great impact on machine learning and signal processing communities given some success stories [START_REF] Shevade | A simple and efficient algorithm for gene selection using sparse logistic regression[END_REF][START_REF] Donoho | Compressed sensing[END_REF][START_REF] Lustig | Compressed sensing MRI[END_REF][START_REF] Ye | Sparse methods for biomedical data[END_REF] it also comes with some From an optimization point of view, the Lasso can benefit from a large palette of algorithms ranging from blockcoordinate descent [START_REF] Friedman | Pathwise coordinate optimization[END_REF][START_REF] Fu | Penalized regressions: the bridge versus the lasso[END_REF] to (accelerated) proximal gradient approaches [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF]. In addition, by leveraging convex optimization theory, some of these algorithms can be further accelerated by combining them with sequential or dynamic safe screening rules [START_REF] El Ghaoui | Safe feature elimination in sparse supervised learning[END_REF][START_REF] Bonnefoy | Dynamic screening: Accelerating first-order algorithms for the lasso and group-lasso[END_REF][START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF], which allow to safely set some useless variables to 0 before terminating (or even sometimes before starting) the algorithm.

While learning sparse models with non-convex penalties are seemingly more challenging to address, blockcoordinate descent algorithms [START_REF] Breheny | Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection[END_REF][START_REF] Mazumder | Sparsenet: Coordinate descent with nonconvex penalties[END_REF] or iterative shrinkage thresholding [START_REF] Gong | A general iterative shrinkage and thresholding algorithm for non-convex regularized optimization problems[END_REF] algorithms can be extended to those learning problems. Another popular way of handling non-convexity is to consider the majorization-minimization (MM) [START_REF] Hunter | A tutorial on mm algorithms[END_REF] principle which consists in iteratively minimizing a majorization of a (non-convex) objective function. When applied to non-convex sparsity enforcing penalty, the MM scheme leads to solving a sequence of weighted Lasso problems [START_REF] Candès | Enhancing sparsity by reweighted l 1 minimization[END_REF][START_REF] Gasso | Recovering sparse signals with a certain family of nonconvex penalties and dc programming[END_REF][START_REF] Mairal | Stochastic majorization-minimization algorithms for large-scale optimization[END_REF].

In this paper, we propose a screening strategy that can be applied when dealing with non-convex penalties. As far as we know, this work is the first attempt in that direction. For that, we consider a MM framework and in this context, our contributions are

• the definition of a MM algorithm that produces a se-quence of iterates known to converge towards a critical point of our non-convex learning problem,

• the proposition of a duality gap based screening rule for weighted Lasso, which is the core algorithm of our MM framework,

• the introduction of conditions allowing to propagate screened variables from one MM iteration to the next,

• we also empirically show that our screening strategy indeed improves the running time of MM algorithms with respect to block-coordinate descent or proximal gradient descent methods.

Global non-convex framework

We introduce the problem we are interesting in, its first order optimality conditions and propose an MM approach for its resolution.

The optimization problem

We consider solving the problem of least-squares regression with a generic penalty of the form

min w∈R d 1 2 y -Xw 2 2 + d j=1 r λ (|w j |) , (1)
where y ∈ R n is a target vector, X = [x 1 , . . . , x d ] ∈ R n×d is the design matrix with column-wise features x j , w is the coefficient vector of the model and the map r λ : R + → R + is concave and differentiable on [0, +∞) with a regularization parameter λ > 0. In addition, we assume that r λ (|w|) is lower semi-continuous function. Note that most penalty functions such as SCAD, MCP or log sum (see their definitions in Table 1) admit such a property.

We consider tools such as Fréchet subdifferentials and limiting-subdifferentials [START_REF] Kruger | On Fréchet subdifferentials[END_REF][START_REF] Rockafellar | Variational analysis[END_REF][START_REF] Mordukhovich | Fréchet subdifferential calculus and optimality conditions in nondifferentiable programming[END_REF] well suited for non-smooth and non-convex optimization, so that a vector w belongs to the set of minimizers (not necessarily global) of Problem (1) if the following Fermat condition holds [START_REF] Clarke | Method of Dynamic and Nonsmooth Optimization[END_REF][START_REF] Kruger | On Fréchet subdifferentials[END_REF]:

X (y -Xw ) ∈ d j=1 ∂r λ (|w j |) , (2)
with ∂r λ (•) being the Fréchet subdifferential of r λ , assuming it exists at w . In particular this is the case for the MCP, log sum and SCAD penalties presented in Table 1. For an illustration, we present the optimality conditions for MCP and log sum.

Example 1. For the MCP penalty (see Table 1 for the definition and subdifferential), it is easy to show that the ∂r λ (0) = [-λ, λ]. Hence, the Fermat condition becomes

-x j (y -Xw ) = 0, if |w j | > λθ |x j (y -Xw )| ≤ λ, if w j = 0 -x j (y -Xw ) + λ sign(w j ) - w j θ = 0, otherwise . (3)
Example 2. For the log sum penalty one can explicitly compute ∂r λ (0) = [-λ θ , λ θ ] and leverage on smoothness of r λ (|w|) when |w| > 0 for computing ∂r λ (|w|). Then the above necessary condition can be translated as

-x j (y -Xw ) + λ sign(w j ) θ + |w j | = 0, if w j = 0 , |x j (y -Xw )| ≤ λ θ , if w j = 0 . (4)
As we can see, first order optimality conditions lead to simple equations and inclusion.

Remark 1. There exists a critical parameter λ max such that 0 is a critical point for the primal problem for all λ ≥ λ max . This parameter depends on the subdifferential of r λ at 0. For the MCP penalty we have λ max max j |x j y|, and for the log sum penalty λ max θ max j |x j y|. From now on, we assume that λ ≤ λ max to avoid such irrelevant local solutions.

Majorization-Minimization approach

There exists several majorization-minimization algorithms for solving non-smooth and non-convex problems involving sparsity-inducing penalties [START_REF] Gasso | Recovering sparse signals with a certain family of nonconvex penalties and dc programming[END_REF][START_REF] Gong | A general iterative shrinkage and thresholding algorithm for non-convex regularized optimization problems[END_REF][START_REF] Mairal | Stochastic majorization-minimization algorithms for large-scale optimization[END_REF].

In this work, we focus on MM algorithms provably convergent to a critical point of Problem 1 such as those described by [START_REF] Kang | On the global convergence of majorization minimization algorithms for nonconvex optimization problems[END_REF]. Their main mechanism is to iteratively build a majorizing surrogate objective function that is easier to solve than the original learning problem. In the case of non-convex penalties that are either fully concave or that can be written as the sum of a convex and concave functions, the idea is to linearize the concave part, and the next iterate is obtained by optimizing the resulting surrogate function. Since r λ (| • |) is a concave and differentiable function on [0, +∞), at any iterate k, we have :

r λ (|w j |) ≤ r λ (|w k j |) + r λ (|w k j |)(|w j | -|w k j |) .
To take advantage of MM convergence properties [START_REF] Kang | On the global convergence of majorization minimization algorithms for nonconvex optimization problems[END_REF], we also majorize the objective function and

Log sum λ log(1 + |w|/θ) -λ θ , λ θ if w = 0 λ sign(w) θ+|w| if w > 0 MCP λ|w| -w 2 2θ if |w| ≤ λθ θλ 2 /2 if |w| > θλ    [-λ, λ] if w = 0 {λ sign(w) -w θ } if 0 < |w| ≤ λθ {0} if |w| > θλ SCAD    λ|w| if |w| ≤ λ 1 2(θ-1) (-w 2 + 2θλ|w| -λ 2 ) if λ < |w| ≤ λθ λ 2 (1+θ) 2 if |w| > θλ        [-λ, λ] if w = 0 {λ sign(w)} if 0 < |w| ≤ λ 1 θ-1 (-w + θλ sign(w)) if 0 < |w| ≤ λθ {0} if |w| > θλ
our algorithm boils down to the following iterative process

w k+1 = arg min w∈R d 1 2 y -Xw 2 2 + 1 2α w -w k 2 2 (5) + d j=1 r λ (|w k j |)|w j | ,
where α > 0 is some user-defined parameter controlling the proximal regularization strength. Note that as α → ∞, the above problem recovers the reweighted 1 iterations investigated by [START_REF] Candès | Enhancing sparsity by reweighted l 1 minimization[END_REF]; [START_REF] Gasso | Recovering sparse signals with a certain family of nonconvex penalties and dc programming[END_REF]. Moreover, when using w k = (e.g., when evaluating the first λ in a path-wise fashion, for k = 0) this recovers the Elastic-net penalty [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF].

Proximal Weighted Lasso: coordinate descent and screening

As we have stated, the screening rule we have developed for regression with non-convex sparsity enforcing penalties is based on iterative minimization of Proximal Weighted Lasso problem. Hereafter, we briefly show that such subproblems can be solved by iteratively optimizing coordinate wise. Then, we derive a duality gap based screening rule to screen coordinate-wise. In what follows, we assume that 0 0 = 0 1 .

Primal and dual problems

Solving the following primal problem (as it encompass Problem 5 as a special case) would prove useful in the application of our MM framework:

min w∈R d P Λ (w) 1 2 y-Xw 2 2 + 1 2α w-w 2 2 + d j=1 λ j |w j | .
(6) with w is some pre-defined vector of R d , and Λ = (λ 1 . . . λ d ) with λ 1 ≥ 0, . . . , λ d ≥ 0 some regulariza-1 This will be of interest cases where λj = 0 in Problem 6. tion parameters. In the sequel we denote Problem 6 as the Proximal Weighted Lasso (PWL) problem.

Problem ( 6) can be solved through proximal gradient descent algorithm [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF] but in order to benefit from screening rules, coordinate descent algorithms are to be privileged. [START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF] have proposed a coordinate-wise update for the Elastic-net penalty and similarly, it can be shown that for the above, the following update holds for any coordinate w j and λ j ≥ 0:

w j ← 1 x j 2 2 + 1 α sign(t j ) max(0, |t j | -λ j ) , (7)
with t j = x j (y -Xw + x j w j ) + 1 α w j . Typically, coordinate descent algorithms visit all the variables in a cyclic way (another popular choice is sampling uniformly at random among the coordinates) unless some screening rules prevent them from unnecessary updates.

Recent efficient screening methods rely on producing primal-dual approximate solutions and on defining tests based on these approximate solutions [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF][START_REF] Shibagaki | Simultaneous safe screening of features and samples in doubly sparse modeling[END_REF][START_REF] Ndiaye | Gap safe screening rules for sparsity enforcing penalties[END_REF]. While our approach follows the same road, it does not result from a straightforward extension of the works of [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF] and [START_REF] Shibagaki | Simultaneous safe screening of features and samples in doubly sparse modeling[END_REF] due to the proximal regularizer.

The primal objective of Proximal Weighted Lasso (given in Equation ( 6)) is convex (and lower bounded) and admits at least one global solution. To derive our screening tests, we need to investigate the dual formulation associated to this problem, that reads:

max s∈R n v∈R d D(s, v) - 1 2 s 2 2 - α 2 v 2 2 + s y -v w (8) s.t. |X s -v| Λ , (9)
with the inequality operator being applied in a coordinate-wise manner. As a side result of the dual derivation (see Appendix for the details), key conditions for screening any primal variable w j are obtained as:

|x j s -v j | -λ j < 0 =⇒ w j = 0 ,
with s and v being solutions of the dual formulation from Equation (8).

Screening test

Screening test on the Proximal Weighted Lasso problem can thus be derived if we are able to provide an upper bound on |x j s -v j | that is guaranteed to be strictly smaller than λ j . Suppose that we have an intermediate triplet of primaldual solution ( ŵ, ŝ, v) with ŝ and v being dual feasible2 , then we can derive the following bound

|x j s -v j | = |x j ŝ -vj + x j (s -ŝ) -(v j -vj )| ≤ |x j ŝ -vj | + x j s -ŝ + |v j -vj | .
Now we need an upper bound on the distance between the approximated and the optimal dual solution in order to make screening condition exploitable. By exploiting the property that the objective function D(s, v) of the dual problem given in Equation ( 8) is quadratic and strongly concave3 , the following inequality holds

D(ŝ, v) ≤D(s , v ) -∇ s D(s , v ) (ŝ -s ) -∇ v D(s , v ) (v -v ) - 1 2 ŝ -s 2 2 - α 2 v -v 2 2 , with ∇ s D = ∂D ∂s1 , . . . , ∂D ∂sn , ∇ v D = ∂D ∂v1 , . . . , ∂D ∂v d .
As the dual problem is a constrained optimization problem, the first-order optimality condition for (s , v ) reads

∇ s D(s , v ) (s -s ) + ∇ v D(s , v ) (v -v ) ≥ 0, ∀s ∈ R n , v ∈ R d ; thus we have 2(D(s , v ) -D(ŝ, v)) ≥ ŝ -s 2 2 + α v -v 2 2 . By strong duality, we have P Λ ( ŵ) ≥ D(s , v ), hence 2(P Λ ( ŵ) -D(ŝ, v)) ≥ ŝ -s 2 2 + α v -v 2 2 .
We can now use the duality gap for bounding ŝs and

|v j -v j |. Hence, given a primal-dual intermediate solu- tion ( ŵ, ŝ, v), with duality gap G Λ ( ŵ, ŝ, v) P Λ ( ŵ) - D(ŝ, v), the screening test for a variable j is |x j ŝ -vj | + 2G Λ ( ŵ, ŝ, v) x j + 1 α T (λ j ) j ( ŵ,ŝ,v) < λ j , (10)
and we can safely state that the j-th coordinate of w is zero when this happens.

Finding approximate primal-dual solutions: In our case of interest, the Proximal Weighted Lasso problem is solved in its primal form using a coordinate descent algorithm that optimizes one coordinate at a time. Hence, an approximate primal solution ŵ is easy to obtain by considering the current solution at a given iteration of the algorithm. From this primal solution, we show how to obtain a dual feasible solution ŝ that can be considered for the screening test.

One can check the following primal/dual link (for instance by deriving the first order conditions of the maximization in Equation 16or equivalent formulation 21, see Appendix)

y -Xw = s and w -w = αv with the constraints that |x j s -v j | ≤ λ j , ∀j ∈ [d].
Hence, a good approximation of the dual solution can be obtained by scaling the residual vector y -X ŵ such that it becomes dual feasible. Indeed, the condition

|x j s -v j | ≤ λ j , ∀j ∈ [d]
is guaranteed only at optimality. To avoid issues with dividing by vanishing λ j 's, let us consider the set S = {j ∈ [d] : λ j > 0} of associated indexes and assume that this set is non-empty. Then, one can define

j † = arg max j∈S 1 λj x j (y -X ŵ) -1 α ( ŵj -w j ) ρ(j)
.

For all j ∈ S, if ρ(j † ) ≤ 1 then , ŝ y -X ŵ and v = 1 α ( ŵ -w ) are dual feasible, no scaling is needed. If ρ(j † ) > 1, we define the approximate dual solution as ŝ = y-X ŵ ρ(j † ) and v = ŵ-w αρ(j † ) which are dual feasible. Hence, in practice, we compute our screening test using the triplet { ŵ, y-X ŵ max(1,ρ(j † )) , ŵ-w α max(1,ρ(j † )) }. Remark 2. As the above dual approximation is valid only for components of λ j > 0, we have added a special treatment for components j ∈ [d] \ S setting vj = x j ŝ for such indexes, guaranteeing dual feasibility. Note that the special case λ j = 0 is not an innocuous case. For some non-convex penalties like MCP or SCAD, their gradient r λ is equal to 0 for large values and this is one of their key statistical property. Hence, the situation in which λ j = 0 is very likely to occur in practice for these penalties.

Comparing with other screening tests: Our screening test exploits duality gap and strong concavity of the dual function, and the resulting test is similar to the one derived by [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF]. Indeed, it can be shown that Equation 10 boils down to be a test based on a sphere centered on ŝ with radius 2G Λ ( ŵ, ŝ, v). The main difference relies on the threshold of the test λ j , which differs on a coordinate basis instead of being uniform.

Similarly to the GAP test of [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF], our screening rule come with theoretical properties. It can be Algorithm 1 Proximal Weighted Lasso (PWL) Input: X, y, w 0 ,Λ, α, LstScreen Output: w, r = y -Xw 1: k ← 0 2: repeat 3:

for variable j ∈ LstScreen do Proof. Because ŵk , ŝk and vk are convergent, owing to the strong duality, the duality gap also converges towards zero. Now, for any given , define k 0 such that ∀k ≥ k 0 , we have

ŝk -s 2 ≤ , vk -v ∞ ≤ and 2G Λ ≤ .
For j ∈ E , we have

|x j (ŝ k -s ) -(v k j -v j )| ≤ |x j (ŝ k -s )| + |(v k j -v j )| ≤ (max j ∈E x j + 1) .
From triangle inequality we have:

|x j ŝk -vk j | ≤ |x j (ŝ k -s ) -(v k j -v j )| + |x j s -v j | ≤ (max j ∈E x j + 1) + |x j s -v j | If we add √ 2G Λ x j + 1 α on both sides, we get T (λj ) j ≤ 2 max j ∈E x j + 1 + 1 α + |x j s -v j | .
4 following a terminology introduced byTibshirani (2013)

Now define the constant C min j ∈E [λ j -|x j s -v j |] and because j ∈ E , C > 0. Hence, if we choose

< C 2 max j ∈E x j + 1 + 1 α ,
we have T (λj ) j < λ j which means that the variable j has been screened, hence j ∈ E k . To conclude, we have j ∈ E implies that j ∈ E k , which also translates in E k ⊂ E . This property thus tells that all the zero variables of w are correctly detected and screened by our screening rule in a finite number of iterations of the algorithm.

Screening rule for non-convex regularizers

Now that we have described the inner solver and its screening rule, we are going to analyze how this rule can be improved into a majorization-minimization (MM) approach.

Majorization-minimization approaches and screening

At first, let us discuss some properties of our MM algorithm. According to the first order condition of the Proximal Weighted Lasso related to the MM problem at iteration k, the following inequality holds for any j ∈ [d]

|x j s k, -v k, j | ≤ λ k j ,
where the superscript denotes the optimal solution at iteration k. Owing to the convergence properties of the MM algorithm [START_REF] Kang | On the global convergence of majorization minimization algorithms for nonconvex optimization problems[END_REF], we know that the sequence {w k } converges towards a vector satisfying Equation (2).

Owing to the continuity of r λ (|w|), we deduce that the sequence {λ k j } also converges towards a λ j . Thus, by taking the limits of the above inequality, the following condition holds for w j = 0:

|x j s -v j | ≤ λ j ,
with λ j = r λ (|w j |). This inequality basically tells us about the relation between vanishing primal component and the optimal dual variable at each iteration k. While this suggests that screening within each inner problem by defining λ j in Equation ( 6) as λ k j = r λ (|w k j |) should improve efficiency of the global MM solver, it does not tell whether screened variables at iteration k are going to be screened at the next iteration as the λ k j 's are also expected to vary between iterations.

Remark 3. In general, the behavior of a λ j across MM iterations strongly depends on an initial w 0 and the optimal w . For instance, if one variable w 0 j is initialized at 0 and w j is large, λ j will tend to be decreasing across iterations. Conversely, λ j will tend to increase.

Propagating screening conditions

In what follows, we derive conditions on allowing to propagate screened coefficients from one iteration to another in the MM framework. Property 2. Consider a Proximal Weighted Lasso problem with weights {λ j } and its primal-dual approximate solutions ŵ, ŝ and v allowing to evaluate a screening test in Equation 10. Suppose that we have a new set of weight Λ ν = {λ ν j } j=1,...,d defining a new Proximal Weighted Lasso problem. Given a primal-dual approximate solution defined by the triplet ( ŵν ,ŝ ν , vν ) for the latter problem, a screening test for variable j reads

T (λj ) j ( ŵ, ŝ, v)+ x j (a+ √ 2b)+c+ 1 α √ 2b < λ ν j , ( 12
)
where T

(λj ) j ( ŵ, ŝ, v) is the screening test for j at λ j , a, b and c are constants such that ŝν -

ŝ 2 ≤ a, |G Λ ( ŵ, ŝ, v) -G Λ ν ( ŵν , ŝν , vν )| ≤ b and |v ν j -vj | ≤ c.
Proof. The screening test for the novel Proximal Weighted Lasso problem for parameter Λ ν j can be written as

|x j ŝν -vν | + 2G Λ ν ( ŵν , ŝν , v ν ) x j + 1 α < λ ν j .
Let us bound the terms in the left-hand side of this inequality. At first, we have:

|x j ŝν -vν | ≤ |x j ŝ -vj | + |x j (ŝ ν -ŝ)| + |v ν j -vj | ≤ |x j ŝ -vj | + x j ŝν -ŝ + |v ν j -vj | , and
G Λ ν ≤ |G Λ ν -G Λ | + |G Λ | (13) ≤ |G Λ ν -G Λ | + |G Λ | .
where the last inequality holds by applying the norm property

x 2 ≤ x 1 to the 2-dimensional vector of compo- nent [ |G Λ ν -G Λ |, |G Λ |],
where we have drop the dependence on ( ŵ, ŝ, v) and ( ŵν , ŝν , vν ) for simplicity. By gathering the pieces together, the left-hand side of Equation 4.2 can be bounded by

|x j ŝ -vj | + x j ŝν -ŝ + 2G Λ x j + 1 α (14) + 2|G Λ ν -G Λ | x j + 1 α + |v ν j -vj | .
which leads us to the novel screening test

T (λj ) j ( ŵ, ŝ, v)+ x j (a+ √ 2b)+c+ 1 α √ 2b ≤ λ ν j . ( 15
)
with T

(λj ) j ( ŵ, ŝ, v) defined in Equation ( 10).

In order to make this screening test tractable, we need at first an approximation ŝν of the dual solution, then an upper bound on the norm of ŝŝν and a bound on the Algorithm 2 MM algorithm for Lasso with penalty r λ (|w|)

Input: X, y, λ, w 0 , α 1: k ← 0 2: repeat 3:

Λ k = {λ k j } j=1,...,
k ← k + 1 14: until convergence difference in duality gap |G Λ ν -G Λ |.
In practice, we will define ŝν = y-X ŵ ρ ν (j † ) and then compute exactly ŝŝν and |G Λ ν -G Λ |. Interestingly, since we consider the primal solution ŵ as our approximate primal solution for the new problem, computing ρ ν (j) is only costs an element-wise division since y -X ŵ and |x j (y -X ŵ)| have already been pre-computed at the previous MM iteration. Another interesting point to highlight is that given pre-computed screening values {T (λj ) j ( ŵ, ŝ, v)}, the screening test given in Equation ( 12) does not involve any additional dot product and thus is cheaper to compute.

Algorithm

When considering MM algorithms for handling nonconvex penalties, we advocate the use of weighted Lasso with screening as a solver for Equation 5 and a screening variable propagation condition as given in Equation 12. This last property allows us to cheaply evaluating whether a variable can be screened before entering in the inner solver. However, Equation ( 12) also needs the screening score computed at some previous {λ j } and a trade-off has thus to be sought between relevance of the test when λ ν j is not too far from λ j and the computational time needed for evaluating T (λj ) j

. In practice, we compute the exact score T (λj ) j ( ŵ, ŝ, v) every 10 iterations and apply Equation ( 12) for the rest of the iterations. This results in Algorithm 2 where the inner solvers, denoted as PWL and solved according to Algorithm 1, are warm-started with previous outputs and provided with a list of already screened variables. In practice, α > 0 helps us guaranteeing theoretical convergence of the sequence {w k } and we have set

Numerical experiments

The goal of screening rules is to improve the efficiency of solvers by focusing only on variables that are nonvanishing. In this section, we thus report the computational gains we obtain owing to our screening strategy.

Experimental set-up

Our goal is to compare the computational running time of different algorithms for computing a regularization path on problem (1). For most relevant non-convex regularizers, coordinatewise update [START_REF] Breheny | Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection[END_REF] and proximal operator [START_REF] Gong | A general iterative shrinkage and thresholding algorithm for non-convex regularized optimization problems[END_REF] can be derived in a closed-form. For our experiments, we have used the logsum penalty which has an hyperparameter θ. Hence, our regularization path involves 2 parameters (λ t ,θ t ). The set of {λ t } Nt-1 t=0 has been defined as λ t λ max 10 -3t N t -1 , N t being dependent of the problems and θ ∈ {0.01, 0.1, 1}.

Our baseline method should have been an MM algorithm, in which each subproblem as given in Equation 5is solved using a coordinate descent algorithm without screening. However, due to its very poor running time, we have omitted its performance. Two other competitors that directly address the non-convex learning problems have instead been investigated: the first one, denoted as GIST, uses a majorization of the loss function and iterative shrinkage-thresholding [START_REF] Gong | A general iterative shrinkage and thresholding algorithm for non-convex regularized optimization problems[END_REF] while the second one, named ncxCD is a coordinate descent algorithm that directly handles non-convex penalties [START_REF] Breheny | Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection[END_REF][START_REF] Mazumder | Sparsenet: Coordinate descent with nonconvex penalties[END_REF]. We have named MM-screen our method that screens within each Proximal Weighted Lasso and propagates screening scores while the genuine version drops the screening propagation (but screening inside inner solvers is kept), is denoted MMgenuine. All algorithms have been stopped when optimality conditions as described in Equation 4 are satisfied up to the same tolerance τ . Note that all algorithms have been implemented in Python/Numpy. Hence, this may give GIST a slight computational advantage over coordinate descent approaches that heavily benefit on loop efficiency of low level languages. For our MM approaches, we stop the inner solver of the Proximal Weighted Lasso when the duality gap is smaller than 10 -4 and screening is computed every 5 iterations. In the outer iterations, we perform the screening every 10 iterations.

Toy problem

Our toy regression problem has been built as follows. The entries of the regression design matrix X ∈ R n×d are drawn uniformly from a Gaussian distribution of zero mean and variance 4. For a given n and d and a number p of active variables, the true coefficient vector w is obtained as follows. The p non-zero positions are chosen randomly, and their values are drawn from a zero-mean unit variance Gaussian distribution, to which we added ±0.1 according to the sign of w j . Finally, the target vector is obtained as y = Xw + e where e is a random noise vector drawn from a Gaussian distribution with zero-mean and standard deviation σ. For the toy problem, we have set N t = 50.

Figure 1 presents the running time needed for the different algorithms to reach convergence under different settings. We note that indifferently to the settings our screening rules help in reducing the computational time by a fac- tor of about 5 compared to a plain non-convex coordinate descent. The gain compared to GIST is mostly notable for high-precision and noisy problems.

We have also analyzed the benefit brought by our screening propagation strategy. Figure 2 presents the gain in computation time when comparing MM screening and MM genuine. As we have kept the number of relevant features to 5, one can note that for a fixed amount of noise the gain is around 1.8 for a wide range of dimensionality. Similarly, the gain is almost constant for a large range of noise and the best benefit occurs in a low-noise setting. More interestingly, we have compared this gain for increasing precision and for a low-noise situation σ = 0.01, which is classical noise level in the screening literature [START_REF] Ndiaye | Gap safe screening rules for sparse-group lasso[END_REF][START_REF] Tibshirani | Strong rules for discarding predictors in lasso-type problems[END_REF]. We can note from the mostright panel of Figure 2 that the more precision we require on the resolution of the learning problem, the more we benefit from screening propagation. For a tolerance of 10 -8 , the gain peaks at about 6.5 whereas in most regimes of number of features, the gain is about 4. Even for a lower tolerance, the genuine screening approach is 4 times slower than the full approach we propose.

Real-world datasets

We have also run the comparison on different real datasets. Figure 3 presents the results obtained on the leukemia dataset, which is a dense data with n = 72 examples in dimension d = 7129. For the path computation, the set of λ t has been fixed to N t = 20 elements. Remark that the gain compared to ncvxCD varies from 3 to 5 depending on the tolerance and is about 4 compared to GIST at high tolerance. On the right panel of Figure 3, we compare the gain in running time brought by screening propagation rule in a real world sparse dataset newsgroups in which we have kept only 2 categories (religion and graphics) resulting in n = 961 and d = 21319. We can note that the gain is similar to what we have observed on the toy problem ranging in between 1.3 and 1.8.

Conclusion

We have presented the first screening rule strategy that handles sparsity-inducing non-convex regularizers. The approach we propose is based on a majorizationminimization framework in which each inner iteration solves a Proximal Weighted Lasso problem. We introduced a screening rule for this learning problem and a rule for propagating screened variables within MM iterations. In-terestingly, our screening rule for the weighted Lasso is able to identify all the variables to be screened in a finite amount of time. We have carried out several numerical experiments showing the benefits of the proposed approach compared to methods directly handling the non-convexity of the regularizers and illustrating the situation in which our propagating-screening rule helps in accelerating efficiency of the solver.

Supplementary Material of "Screening rules for Lasso with non-convex sparse regularizers" Dual problem of the weighted Lasso minimization Let recall the weighted lasso problem Let ãi = X i,: and φ i (z i ) = 1 2 (ỹ i -z i ) 2 being the quadratic loss function. Let its convex conjugate [START_REF] Boyd | Convex Optimization[END_REF] being φ * i (η i ) = max zi η i z i -φ i (z i ) for a scalar η i , which results in φ * i (η i ) = 1 2 η 2 i + η i ỹi . Note also that φ i = φ * * i as φ i is convex.

Following [START_REF] Johnson | Blitz: A principled metaalgorithm for scaling sparse optimization[END_REF] we derive the dual of the weighted problem through these steps

The dual objective function is obtained by substituting the expression φ * i and using the optimality condition of the problem min

w∈R d w X η + d i=1 λ j |w j | . ( 18
)
This problem is separable and optimality condition with respect to any w j is as follows, provided g j = X η j g j + λ j sign(w j ) = 0 if w j = 0

|g j | ≤ λ j if w j = 0. ( 19
)
The latter condition implies the coordinate-wise inequality constraint | X η| Λ with Λ = λ 1 , . . . , λ d . Also we can easily establish that g j w j + λ j |w j | = 0. Hence the objective function in Equation ( 18) vanishes.

Finally let us decompose the dual vector as η = -s √ αv where s ∈ R n and v ∈ R d . Recalling the form of ỹ and X, it is easy to see that the dual problem (17) becomes max s,v:|X s-v| Λ -1 2 s 2 2 -α 2 v 2 2 + s yv w .

Also, from (19), it holds the screening conditions

|x j s -v j | < λ j =⇒ w j = 0 , ∀j ∈ [d] , (20)
remind that x j = X :,j is the j th covariate. In addition, the maximisation in Equation ( 16) takes the form

max s,v -1 2 s 2 2 -α 2 v 2 2 + s (y -Xw) + v (w -w ) . (21)
Thus given an optimal solution w , we may have s = y -Xw and αv = w -w , by deriving the first order optimality conditions of this maximization problem.

  convergenceshown that if we have a converging sequence { ŵk } of primal coefficients lim k→∞ ŵk = w , then ŝk and vk defined as above converge to s and v . Moreover, we can state a property showing the ability of our screening rule to remove irrelevant variables after a finite number of iterations.Property 1. Define the Equicorrelation set 4 of the Proximal Weighted Lasso asE {j ∈ [d] : |x j s -v j | = λ j } and E k {j ∈ [d] : |x j ŝk -vk j | ≥ λ j }obtained at iteration k of an algorithm solving the Proximal Weighted Lasso. Then, there exists an iteration k 0 ∈ N s.t. ∀k ≥ k 0 , E k ⊂ E .

Figure 1 .

 1 Figure 1. Comparing running time for coordinate descent (CD) and proximal gradient algorithms as well as screening-based CD and MM algorithms under different tolerances on the stopping condition and under different samples, features settings. For all experiments, we have 5 active variables. (left) n = 50, d = 100, σ = 2. (right) n = 500, d = 5000, σ = 2

Figure 3 .

 3 Figure 3. Running time for computing the regularization path on (left) the dense Leukemia dataset n = 72, d = 7129. (right) sparse dataset Newsgroup with n = 961 and d = 21319.

  i w)η i -φ * i (η i )] +

Table 1 .

 1 Common non-convex penalties with their sub-differentials. Here λ > 0, θ > 0 (θ > 1 for MCP, θ > 2 for SCAD).

	Penalty	r λ (|w|)	∂r λ (|w|)

  d ← {r λ (|w k j |)} j=1,...,d

	4:	compute approximate dual s k and duality gap G Λ k
		given w k and Λ k
	5: 6:	if needed then compute screening scores T j (λj )	(w k , s k , v k ) ac-
		cording to Eq. (10)
	7:	store reference duality gap and approximate dual
	8:	else
	9:	estimate screening scores according to Eq. (12)
	10:	end if
	11:	LstScreen ← updated screened variables list based
		on results of Line 6 or 9
	12:	w

k+1 , y -Xw k+1 ← PWL(X,y,w k ,Λ k ,α,Lstcreen) 13:

  Figure 2. Computational gain of propagating screening within MM iterations (left) w.r.t. number of features at a KKT tolerance of 10 -4 . (middle) w.r.t. to noise level at a KKT tolerance of 10 -4 . (right) w.r.t. number of features at different tolerance in a low-noise regime.

		2.0 2.2			2.75 3.00								7 8		tol=1e-4 tol=1e-6 tol=1e-8
		1.8			2.50								6
	Gain	1.6		Gain	2.25							Gain	5
		1.4			1.75 2.00								4
		1.2			1.50								3
		0 2000 4000 6000 8000 10000 12000 14000 Number of features		0	1	2	3	4	5	6	7	0	2000	4000 Number of features 6000	8000	10000
		1.00e-03	1.00e-04 Tolerance	1.00e-05

Université Rouen Normandie, Criteo AI Lab

INSA Rouen Normandie

Université de Montpellier. Correspondence to: alain.rakotomamonjy <first.last@insa-rouen.fr>.February 13, 2019 theoretical drawbacks (e.g., biased estimates of large coefficient of the model). Hence, several authors have proposed non-convex penalties that approximate better the 0 -(pseudo)norm, the later being the original measure of sparsity though it leads to NP-hard learning problem. The most commonly used non-convex penalties are the Smoothly Clipped Absolute Deviation (SCAD)[START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF], the Log Sum penalty (LSP)[START_REF] Candès | Enhancing sparsity by reweighted l 1 minimization[END_REF], the capped-1 penalty(Zhang, 2010), the Minimax Concave Penalty (MCP)(Zhang et al., 2010). We refer the interested reader to[START_REF] Soubies | A unified view of exact continuous penalties for 2 -0 minimization[END_REF] for a discussion on the pros and cons of such non-convex formulations.

meaning |X ŝ -v| Λ.

see[START_REF] Nesterov | Introductory lectures on convex optimization[END_REF] for a precise definition of strong convexity/concavity.