C Elvira

C Herzet

Safe Squeezing for Antisparse Coding Clément Elvira and Cédric Herzet

Abstract-Spreading the information over all coefficients of a representation is a desirable property in many applications such as digital communication or machine learning. This socalled antisparse representation can be obtained by solving a convex program involving an ∞-norm penalty combined with a quadratic discrepancy. In this paper, we propose a new methodology, dubbed safe squeezing, to accelerate the computation of antisparse representation. We describe a test that allows to detect saturated entries in the solution of the optimization problem. The contribution of these entries is compacted into a single vector, thus operating a form of dimensionality reduction. We propose two algorithms to solve the resulting lower dimensional problem. Numerical experiments show the effectiveness of the proposed method to detect the saturated components of the solution and illustrates the induced computational gains in the resolution of the antisparse problem.

Index Terms-antisparse coding, safe screening, scaled projected-gradient algorithm, Frank-Wolfe algorithm.

I. INTRODUCTION

I N the last decades, convex optimization has become a central element in the resolution of inverse problems. This success revolves around two main ingredients. First, it has been shown that proper "regularizing" functions enforce desirable properties on the solutions of convex problems. One of the first (and maybe most striking) example of such a behavior is the use of " 1 -norm" regularization, promoting sparsity of the solutions. The second ingredient which sparked the success of convex optimization is the advent of numerical procedures able to solve efficiently (up to some accuracy) problems involving thousands to billions of variables. To name a few, let us mention the Augmented Lagrangian methods [START_REF] Rockafellar | A dual approach to solving nonlinear programming problems by unconstrained optimization[END_REF], the forward-backward splitting [START_REF] Bruck | An iterative solution of a variational inequality for certain monotone operators in hilbert space[END_REF], or the Alternating direction method of multipliers [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF].

Of particular interest in this paper is an acceleration method first proposed by El Ghaoui et al. in [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF] in the context of sparsity-promoting convex problems, namely "safe screening". This procedure leverages two main elements. First, the solutions of convex problems involving 1 regularization are typically sparse, i.e., contains a large number of zeros. For example, the solutions of the well-known "Lasso" problem x ∈ arg min

x∈R n 1 2 y -Ax 2 2 + λ x 1 , (1)
where y ∈ R m is an observation vector, A ∈ R m×n a representation matrix and λ a penalization parameter, are known to contain at most m nonzero components, see [START_REF] Foucart | A Mathematical Introduction to Compressive Sensing[END_REF]Theorem 3.1] or [START_REF] Rosset | Boosting as a regularized path to a maximum margin classifier[END_REF]Theorem 3]. Second, if the position of (some of) the zeros in x λ are known, (1) can be transformed into a problem of smaller dimension where the columns of A associated to the position of known zeros are simply discarded. Solving this reduced problem may then potentially result in huge memory and computational savings. As far as sparsitypromoting problems are concerned, the idea of safe screening thus consists in designing simple tests allowing to detect zeros of the solution.

Since the seminal work by El Ghaoui et al., safe screening has been developed and improved in many contributions of the literature, see e.g., [START_REF] Fercoq | Mind the duality gap: safer rules for the Lasso[END_REF][START_REF] Bonnefoy | Dynamic Screening: Accelerating First-Order Algorithms for the Lasso and Group-Lasso[END_REF][START_REF] Malti | Safe screening tests for Lasso based on firmly non-expansiveness[END_REF][START_REF] Herzet | Gather and conquer: Region-based strategies to accelerate safe screening tests[END_REF]. Recently, this strategy has also been extended to other families of regularizers [START_REF] Ndiaye | Gap safe screening rules for sparse multi-task and multi-class models[END_REF][START_REF]Gap safe screening rules for sparsity enforcing penalties[END_REF][START_REF] Rakotomamonjy | Screening rules for Lasso with non-convex sparse regularizers[END_REF]. However, to the best of our knowledge, all these contributions have focused so far on the resolution of "sparsity-promoting" problems, where the solutions contain many zeros. In this paper, we show that the principles ruling safe screening may be extended to another family of convex optimization problems. More specifically, we address the following optimization problem

x ∈ arg min

x∈R n 1 2 y -Ax 2 2 + λ x ∞ , ( 2
-P ∞ λ )
where the admissible range of the coefficients x is penalized through an ∞ -norm. Here, the ∞ -norm is advocated for spreading the information over all representation coefficients in the most uniform way. For this reason, the solutions of (2 -P ∞ λ ) are sometimes referred to as "antisparse" or "spread" since, contrary to the Lasso problem, they are known to be dense with many "saturated" entries satisfying

|x (i)| = x ∞ . (3)
Although less popular than its "sparsity-promoting" 1 counterpart, the "spreading" property of the ∞ norm has proved to be of practical interest in several applicative domains, e.g., to design robust analog-to-digital conversion schemes [START_REF] Cvetković | Resilience properties of redundant expansions under additive noise and quantization[END_REF][START_REF] Calderbank | The pros and cons of democracy[END_REF], reduce the peak-to-average power ratio in multi-carrier transmissions [START_REF] Farrell | A Kashin approach to the capacity of the discrete amplitude constrained Gaussian channel[END_REF][START_REF] Ilic | Papr reductioni in ofdm using kashin's representation[END_REF][START_REF] Studer | PAR-aware large-scale multi-user MIMO-OFDM downlink[END_REF], perform approximate nearest neighbor search [START_REF] Jegou | Anti-sparse coding for approximate nearest neighbor search[END_REF], outlier detection [START_REF] Vural | A new outlier detection method based on anti-sparse representations[END_REF] or robust beamforming [START_REF] Jiang | Large-scale robust beamforming via ∞ -minimization[END_REF].

Besides, when combined with a set of linear equality constraints, minimizing an ∞ -norm is referred to as the minimumeffort control problem in the optimal-control framework [START_REF] Neustadt | Minimum effort control systems[END_REF][START_REF] Cadzow | Algorithm for the minimum-effort problem[END_REF].

From a numerical point of view, although many generic convex optimization tools may be straightforwardly applied to problem (2 -P ∞ λ ), the design of algorithms specifically dedicated to this kind of problems has also focused less attention than its 1 counterpart. In this contribution, we make one step in that direction by showing that the main principles underlying "safe screening" may be extended to problem (2 -P ∞ λ ). More specifically, we emphasize that (2 -P ∞ λ ) can be reduced to a problem of smaller dimension by identifying "saturated" components of the solution x and propose simple tests to identify them. We refer to the proposed methodology as "safe squeezing" since it allows to compact several components of x into one single variable. We show numerically that the proposed safe squeezing procedure can significantly reduce the computational complexity needed to solve the antisparse problem (2 -P ∞ λ ) to some accuracy. The rest of this paper is organized as follows. Section II defines the main notations used in this work. In Section III, we introduce the proposed safe squeezing test and discuss the problem's dimensionality reduction it allows. In Section IV we describe two algorithms to solve the equivalent "reduced" optimization problem. In Section V, we illustrate the performance of the proposed methodology in numerical experiments. Concluding remarks are finally made in Section VI.

II. NOTATIONS

Unless otherwise specified, we will use the following convention throughout the paper. The vectors are denoted by lowercase bold letters (e.g., x) and matrices by uppercase bold letters (e.g., A); x(i) refers to the ith component of x and a i to the ith column of A. Given q ∈ R q and w ∈ R, the shorthand notations "q ≤ w" should be understood as "q(i) ≤ w ∀i". The notations ker(A) and kruskal(A) will be used to denote the null space and the Kruskal rank [5, p. 56] of A, respectively.

Caligraphic letters are used to denote sets (e.g., I). Ī stands for the complementary set of I and card(I) refers to the cardinality of I.

Given a vector x ∈ R n and a set of indices I ⊆ {1, . . . , n}, we let x I be the vector of components of x with indices in I. Similary, A I denotes the submatrix of A whose columns have indices in I. Finally, σ min (A) and σ max (A) refer to the lowest and highest eigenvalue of A, respectively.

III. DETECTING AND EXPLOITING SATURATION

A. Working hypotheses

In the rest of this paper, we will assume that the following working hypothesis is verified:

0 < λ < A T y 1 . (4)
The left inequality is natural since letting λ = 0 is tantamount to removing the penalization of problem (2 -P ∞ λ ). On the other hand, the right inequality ensures that the solution of problem (2 -P ∞ λ ) is not the all-zero vector 0 n . More specifically, it can be shown that

x ∞ > 0 ⇐⇒ λ verifies (4). (5)
A proof of this result is available in Appendix A. Hypothesis (4) implies in particular that y / ∈ ker(A T ) since the two inequalities cannot be satisfied simultaneously in the opposite case. In the sequel, we will always assume that (4) holds although not explicitly mentioned in our statements.

We will also assume that (2 -P ∞ λ ) admits a unique minimizer x . Although our subsequent derivations may be quite easily extended to the general case, this working hypothesis greatly simplifies the exposition of our procedure. We will add different comments all along our presentation to point out how the non-uniqueness of the solutions of (2 -P ∞ λ ) would modify the quantities and results at stake.

B. Dimensionality reduction via saturation detection

In this section, we illustrate how the knowledge of the positions of saturated entries in x can lead to memory and complexity savings in the resolution of (2 -P ∞ λ ). Let1 I + {i :

x (i) = + x ∞ }, (6a)
I - {i : x (i) = -x ∞ }, (6b)
be the sets of positive and negative saturated components of x , and I = I -∪ I + . Then, for any I = I + ∪ I -with I + ⊆ I + and I -⊆ I -, problem (2 -P ∞ λ ) can be equivalently rewritten as

(w , q ) ∈ arg min (w,q)∈R×R n-card(I) 1 2 y -AĪq -s w 2 2 + λw subject to q ≤ w, -q ≤ w (7 -P sq λ ) where s i∈I+ a i - i∈I- a i . (8)
In particular, the following bijection holds between the solutions of (2 -P ∞ λ ) and (7 -P sq λ ):

x Ī = q (9) x (i) = -w if i ∈ I - +w if i ∈ I + . (10)
In other words, any minimizer of (2 -P ∞ λ ) defines a minimizer of (7 -P sq λ ), and vice-versa. The rationale behind problem (7 -P sq λ ) is as follows. Since any minimizer of (2 -P ∞ λ ) weighs (up to a sign) the elements in I by the same value (namely x ∞ ), the signed summation of the columns of A I can be done once for all in advance. This is the meaning of vector s in [START_REF] Bonnefoy | Dynamic Screening: Accelerating First-Order Algorithms for the Lasso and Group-Lasso[END_REF]. Variable w plays the role of the largest absolute value of the elements of x. In particular, the constraints "q ≤ w" and "-q ≤ w" ensure that the absolute value of the coefficients weighting the elements of AĪ is no larger than w.

We note that, although defining the same set of solutions, the dimensionality of problems (2 -P ∞ λ ) and (7 -P sq λ ) can be quite different: whereas the former manipulates n-dimensional variables, the latter only involves an optimization space of dimension n -card(I) + 1. As a limit case, when I = I , the size of problem (7 -P sq λ ) drops down to n-card(I )+1. The next lemma shows that the dimensionality reduction can be drastic in this case:

2 Lemma 1. If kruskal(A) = m then card(I ) ≥ n -m + 1.
Lemma 1 is a direct consequence of Theorem 1 below. The hypothesis "kruskal(A) = m" holds as soon as any group of m columns of A ∈ R m×n is linearly independent [5, p. 56]. This is, for example, a typical setup in machine learning applications where the columns of matrix A contain randomly generated features. 3 In this case, we see from Lemma 1 that the dimension of problem (7 -P sq λ ) is no larger than m if I = I . In the overcomplete setting, when m n, solving (7 -P sq λ ) instead of (2 -P ∞ λ ) then leads to dramatic dimensionality reduction.

A decrease of the size of the problem leads de facto to a reduction of the storage and computational costs necessary to evaluate the solution. Tackling problem (7 -P sq λ ) instead of (2 -P ∞ λ ) may thus be obviously of interest to save computational resources. The relevance of (7 -P sq λ ) is however conditioned to the identification of a (large) subset of I . In the rest of this section, we propose a procedure, dubbed "safe squeezing", to perform efficiently this task. The term "safe" refers to the fact that all the indices identified by our procedure necessarily belong to I . The term "squeezing" relates to the dimensionality reduction allowed by the identification of (some of) the elements of I .

C. Safe squeezing test

In this section, we provide the main ingredient of our safe squeezing procedure. It takes the form of a simple test ensuring that an index i belongs to I , and is grounded on the following result:

Theorem 1. Let I ⊆ I and u = arg max u∈U I 1 2 y 2 2 -1 2 y -u 2 2 (11-D ∞ λ )
where

U I u : A T I u 1 + s T u ≤ λ ( 12
)
and s is defined as in [START_REF] Bonnefoy | Dynamic Screening: Accelerating First-Order Algorithms for the Lasso and Group-Lasso[END_REF]. Then,

∀i : a T i u > 0 ⇒ i ∈ I sign(a T i u ) . (13)
Moreover, if kruskal(A) = m, a T i u = 0 for at least n-m+1 a i 's.

A proof of this result is available in Appendix B. Theorem 1 provides an appealing way of identifying saturated components in x : we have from (13) that i ∈ I as soon as a i is not orthogonal to some vector u defined in (11-D ∞ λ ). Moreover, if kruskal(A) = m, then at least n -m + 1 saturated components of x can be identified by verifying that a T i u > 0.

Unfortunately, in practice, finding a maximizer of (11-D ∞ λ ) turns out to be as difficult as solving our target optimization problem (7 -P sq λ ). In particular, if I = ∅, addressing (11-D ∞ λ ) requires the same order of complexity as solving (2 -P ∞ λ ). Hence, the direct use of ( 13) is of poor interest to save computational resources. Nevertheless, we emphasize below that relaxed versions of ( 13) can be devised to identify subsets of I with a low computational burden.

Before giving a detailed description of the proposed method, we make two important remarks about u . First, we see from (11-D ∞ λ ) that u is defined as a maximizer of a feasible4 and strictly-concave problem with a continuous and coercive cost function. Hence, u always exists [START_REF] Bertsekas | Nonlinear Programming[END_REF]Propositions A.8] and is unique [START_REF] Bertsekas | Nonlinear Programming[END_REF]Propositions B.10]. This justifies the equality sign used in (11-D ∞ λ ). In fact, as shown in Appendix B-A, (11-D ∞ λ ) corresponds to the Lagrangian dual of problem (7 -P sq λ ). Second, the value of u does not depend on the particular subset I ⊆ I considered in (11-D ∞ λ ). This fact is only stated here for the sake of keeping the discussion at a reasonable level of complexity, but is proved in Appendix B. We note that, although the value of u is independent of I, there may be some computational advantages in considering I's with large numbers of elements. In particular, we emphasize in Section III-D below that the computational cost of the "dual scaling" operation, needed to implement the proposed safe squeezing test, evolves linearly with card( Ī).

We are now ready to expose the main building block of our safe squeezing procedure. Let S be some subset of R m such that

u ∈ S. (14)
For now, we only assume that some S verifying ( 14) is available. We will see in the next section how to construct such a region. A region satisfying ( 14) is usually referred to as "safe region" in the screening literature, and the same convention will be used hereafter. We note that since u ∈ S, we have

min u∈S a T i u < a T i u max u∈S a T i u > a T i u , (15)
so that

min u∈S a T i u > 0 ⇒ a T i u > 0 max u∈S a T i u < 0 ⇒ a T i u < 0. (16)
Using ( 13), we thus obtain

min u∈S a T i u > 0 ⇒ i ∈ I + max u∈S a T i u < 0 ⇒ i ∈ I -. (17)
Hence, if the maximum and minimum values of a T i u over S are easy to evaluate, the left-hand side of ( 17) provides a simple test to determine whether i belongs to I + or I -. In particular, considering a safe region S with a spherical geometry, that is

S = B(c, r) {u ∈ R m : u -c 2 ≤ r}, (18)
for some c ∈ R m and r > 0, leads to min u∈B(c,r)

a T i u = a T i c -a i 2 r max u∈B(c,r) a T i u = a T i c + a i 2 r. (19)
Plugging these expressions into (17) yields the following result:

Theorem 2 (Safe sphere squeezing test). If u ∈ B(c, r), then

a T i c > r a i 2 ⇒ i ∈ I sign(a T i c) . (20)
Squeezing test [START_REF] Vural | A new outlier detection method based on anti-sparse representations[END_REF] provides a very practical way of identifying some elements of I . In particular, testing whether i ∈ I only requires to evaluate one inner product between a i and the center c of a safe sphere. We note however that passing ( 20) is only a sufficient condition for i ∈ I . In practice, depending on the choice of the center and the radius of the safe region, some elements of I may not be identified by the proposed safe squeezing test. As a general rule of thumb, "small" safe regions will lead to squeezing tests able to identify more saturated components. We will illustrate this behavior in our numerical results in Section V. We elaborate on the construction of "good" safe sphere regions in the next subsection. For now, let us just notice that if c = u and r = 0, S = B(c, r) is safe and its volume is equal to zero. In this case, the safe squeezing test (20) boils down to [START_REF] Rakotomamonjy | Screening rules for Lasso with non-convex sparse regularizers[END_REF]. In particular, provided that kruskal(A) = m, at least n -m + 1 saturated components of x can be identified (see last part of Theorem 1). As a final remark, let us mention that Theorem 2 corresponds to one specific particularization of (17) (i.e., to spherical regions S). As with procedures derived for Lasso screening, it is also possible to devise squeezing tests based on safe regions having more refined geometries (e.g., dome [START_REF] Xiang | Fast Lasso screening tests based on correlations[END_REF] or truncated dome [START_REF] Xiang | Screening tests for Lasso problems[END_REF]). The extension of our safe squeezing procedure to these geometries is however left for future work and will not be considered hereafter.

D. Construction of safe spheres

A cornerstone of the safe squeezing methodology presented in Section III-C is the identification of a "safe" region, that is a region S verifying [START_REF] Cvetković | Resilience properties of redundant expansions under additive noise and quantization[END_REF]. We elaborate on this problem hereafter: we emphasize that several safe spheres derived in the context of safe screening for Lasso can be reused for safe squeezing with minor modifications. We focus here on the spheres "ST1" and "GAP" respectively proposed in [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF] and [START_REF] Fercoq | Mind the duality gap: safer rules for the Lasso[END_REF].

Our reasoning is based on two key observations. First, the Lasso dual problem shares the same cost function as problem (11-D ∞ λ ) but with a definition of the dual feasible set different from [START_REF]Gap safe screening rules for sparsity enforcing penalties[END_REF], see e.g., [9, Eq. ( 2)]. Second, the safe spheres proposed in [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF][START_REF] Fercoq | Mind the duality gap: safer rules for the Lasso[END_REF] are based on convex optimality conditions and the knowledge of some dual feasible point u, but never on the specific definition of the dual feasible set. Hence, the center and radius of the "ST1" and "GAP" spheres proposed in [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF] and [START_REF] Fercoq | Mind the duality gap: safer rules for the Lasso[END_REF] can be reused here with the provisio that the dual variable u appearing in these expressions is feasible for (11-D ∞ λ ) (that is u ∈ U I ). For example, the ST1 safe sphere derives from the fact that, by definition of u :

∀u ∈ U I : y -u 2 ≥ y -u 2 . (21)
Hence, the sphere B(c, r) with

c = y r = y -u 2 (22)
is safe for problem (11-D ∞ λ ). We see that expression of the center and the radius in [START_REF] Neustadt | Minimum effort control systems[END_REF] are the same as those proposed in [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF] with the difference that u must here be feasible for (11-D ∞ λ ) rather than for the Lasso dual problem. Similarly, the expression of the GAP safe sphere can be extended as follows in the context of safe squeezing: for any I ⊆ I and (w, q, u) verifying

q ≤ w, -q ≤ w, u ∈ U I , (23)
the sphere B(c, r) with

c = u r = 2 gap(w, q, u) (24)
and

gap(w, q, u) = 1 2 y -AĪq -s w 2 2 + λw -( 1 2 y 2 2 -1 2 y -u 2 2 ) (25)
is safe for problem (11-D ∞ λ ). The proof of this result is similar to that provided in [START_REF] Fercoq | Mind the duality gap: safer rules for the Lasso[END_REF] and is omitted here. As for the sphere ST1, the expressions of the center and radius remain the same as for the Lasso problem but the dual (resp. primal) variable(s) coming into play must now be feasible for problem (11-D ∞ λ ) (resp. (7 -P sq λ )). The only computational difficulty in evaluating ( 22) and ( 24) stands in the identification of a dual feasible point. We may nevertheless identify such a point by extending the "dual scaling" procedure proposed in [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF]Section 3.3] to the present setup. More precisely, given any z ∈ R m , we clearly have that

u = dualscal(z) (26)
where

dualscal(z) = z if A T Ī z 1 + s T z ≤ 0, λ A T Ī z 1 +s T z z otherwise.
is dual feasible. We note that the implementation of ( 26) requires the computation of card( Ī)+1 inner products in R m .

E. Static versus dynamic squeezing

Similarly to screening methods, the proposed safe squeezing procedure can be used either in a "static" or a "dynamic" way.

Static squeezing refers to the case where ( 20) is applied once for all (with I = ∅) on the columns of matrix A before the application of a numerical optimization procedure. In such a case, the dual feasible point u used to build the safe sphere is commonly defined as follows:

u = dualscal(y). ( 27
)
Static screening may for example be of interest to reduce the problem's dimensionality so that its size fits the computer's memory requirements. Dynamic squeezing consists in interleaving test [START_REF] Vural | A new outlier detection method based on anti-sparse representations[END_REF] with the iterations of an optimization procedure. Here, the goal is to refine the quality of the safe region and increase the cardinality of I all along the optimization process. The general principle of dynamic squeezing (applied to the resolution of (2

-P ∞ λ )) is described in Algorithm 1. At each iteration, Algorithm 1: Principle of dynamic squeezing // Initialization 1 x (0) = 0 n , I (0) = ∅ 2 u (0) = dualscal(y) 3 t = 1 // Iteration index 4 repeat
// Squeezing test 5 (c (t) , r (t) ) = sphere param(x (t-1) , u (t-1) , I (t-1) )

6

I (t-½) = squeezing test(c (t) , r (t) ) 7 I (t) = I (t-½) ∪ I (t-1)
// Iterations of the optimization procedure

8 (x (t) , u (t) ) = optim update(x (t-1) , u (t-1) , I (t) )
// Update iteration index

9 t = t + 1 10 until convergence criterion is met Output: x (t) , I (t)
the parameters of a safe sphere are evaluated by using a couple a primal-dual feasible points (x (t-1) , u (t-1) ), see line 5. The safe sphere is then exploited in the squeezing test [START_REF] Vural | A new outlier detection method based on anti-sparse representations[END_REF] to identify some set I (t-½) ⊆ I , see line 6. In line 7, the indices identified at the current iteration are merged to those previously identified. Finally, the primal-dual feasible couple is updated in line 8. This operation is typically carried out by running a few iterations of a numerical procedure addressing problem (7 -P sq λ ) with I = I (t) and initialized with the previous iterate (x (t-1) , u (t-1) ). In the common case where the optimization method only returns a primal iterate x (t) , a feasible dual point can for example be computed by dual scaling of z = y -Ax (t) . By construction, we have

I (0) ⊆ I (1) ⊆ . . . ⊆ I (t) (28)
so that the dimension of the reduced problem (7 -P sq λ ) considered in line 8 does not increase through the iterations of Algorithm 1. In particular, if (x (t) , u (t) ) converges to (x , u ), the parameters of the GAP safe sphere converge to (c = u , r = 0) and, from Theorem 1, card(I (t) ) ≥ n-m+1 after a sufficient number of iterations when kruskal(A) = m.

IV. ALGORITHMIC SOLUTIONS FOR (7 -P sq λ ) A crucial difference between "screening" and "squeezing" is the nature of the optimization problems obtained after reduction. Whereas the reduced problem has the same mathematical form as the original one in the context of screening, the reduced problem (7 -P sq λ ) obtained after squeezing is structurally different from (2 -P ∞ λ ). As a consequence, numerical procedures devised to solve (2 -P ∞ λ ) (e.g., [START_REF] Studer | PAR-aware large-scale multi-user MIMO-OFDM downlink[END_REF]) cannot be used straightforwardly to address (7 -P sq λ ). In this section, we describe two algorithmic solutions for problem (7 -P sq λ ): in Section IV-A a numerical procedure based on the Frank-Wolfe algorithm [25, Section 2.2] is proposed; in Section IV-B, we present a strategy based on a rescaled projected gradient method [25, Section 2.3]. These Algorithm 2:

Frank-Wolfe algorithm for problem [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF].

Input: y, AĪ, s, w, w 0 , q 0 1 w (0) = w 0 , q (0) = q 0 // Initialization 2 t = 1 // Iteration index 3 repeat // Search of a feasible descent direction 4 w (t-½) = w 5 q (t-½) = w sign(A T Ī (y -AĪq (t-1) -s w (t-1) ))

// Convex update 6 w (t) = γw (t-1) + (1 -γ)w (t-½) 7 q (t) = γq (t-1) + (1 -γ)q (t-½)
// Update iteration index

8 t = t + 1 9 until convergence criterion is met Output: w (t) , q (t)
procedures will be used in the next section to assess computational gain allowed by the proposed safe squeezing procedure.

A. Frank-Wolfe algorithm

Let us consider the following equivalent formulation of (7 -P sq λ ): min

(w,q)∈R×R n-card(I) 1 2 y -AĪq -s w 2 2 + λw subject to q ≤ w, -q ≤ w, w ≤ w, ( 29
)
where w is some constant such that w ≤ w. We will discuss the identification of w at the end of this section. For now, let us notice that the additional constraint "w ≤ w" does not change the minimizers of the problem but makes the feasible set compact. As consequence, the conditional gradient method, also known as Frank-Wolfe algorithm [START_REF] Frank | An algorithm for quadratic programming[END_REF], can be applied to find a numerical solution of [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF].

The recursions of the Frank-Wolfe algorithm particularized to problem (29) are described in Algorithm 2. One iteration of the procedure consists of two main steps: the identification of a feasible descent direction and the update of the current iterate.

The feasible descent direction is identified as the solution of the following optimization problem:

(w (t-½) , q (t-½) ) = arg min (w,q) w q T g (t-1) (30)
subject to q ≤ w, -q ≤ w w ≤ w, where g (t-1) is the gradient of the cost function evaluated at the current iterate (w (t-1) , q (t-1) ). Problem ( 30) is linear and admits the closed-form solution given in lines 4-5 of Algorithm 2.

The update of the current iterate is performed via a convex combination of (w (t-½) , q (t-½) ) and (w (t-1) , q (t-1) ), see lines 6-7 in Algorithm 2. Parameter γ is chosen to minimize the cost function over the line segment joining ( w (t-½) , q (t-½) ) and (w (t-1) , q (t-1) ). In the particular setup considered here, this amounts to finding the value of γ minimizing a quadratic function over [0, 1], which admits a closed-form solution.

We conclude this section by showing that parameter w can be chosen as follows:

w = 1 2λ y -AĪq -s w 2 2 + w (31)
where (w, q) is any feasible point of problem (7 -P sq λ ). This follows from the following observation:

λw ≤ 1 2 y -AĪq -s w 2 2 + λw ≤ 1 2 y -AĪq -s w 2 2
+ λw, where the last inequality is valid for any couple (w, q) feasible for (7 -P sq λ ). Hence, the definition of w in (31) verifies "w ≤ w" as required at the beginning of this section. In particular, (w, q) = (0, 0 card( Ī) ) is feasible for (7 -P sq λ ) and leads to w = y 2 2 2λ . We will consider the latter choice in our simulations in Section V.

B. Rescaled projected gradient method

We consider a rescaled version of (7 -P sq λ ): min

( w,q)∈R×R n-card(I) 1 2 y -AĪq -s α w 2 2 + λ α w (32)
subject to αq ≤ w, -αq ≤ w

where α > 0. The choice of α will be discussed later on in this section and will be made based on convergence rate arguments.

Particularizing the gradient projection algorithm described in [START_REF] Bertsekas | Nonlinear Programming[END_REF]Section 2.3] to problem (32) leads to Algorithm 3. Each iteration is divided into three main stages: a "gradient descent" step (lines 4-6), a "projection" operation (line 8) and a "convex update" (line 9-10).

Starting from a current estimate (q (t-1) , w (t-1) ), the gradient step consists in updating the current iterate in the direction of the negative gradient of the cost function (by an amount η > 0). The step η is chosen to minimize the cost function along the direction of the (negative) gradient. We note that since the cost function is quadratic, the value of η admits a simple closed-form expression.

The projection step consists in solving the optimization problem specified in line 8 of Algorithm 3. Although conceptually simpler than (32), the latter does not admit any closed-form solution. We show nevertheless in Appendix C that the procedure described in Algorithm 4 can compute the unique minimizer of this problem in a finite number of steps (upper-bounded by card( Ī)) with a complexity scaling as O(card( Ī) log(card( Ī))). The interested reader is referred to Appendix C-A for more details on the computational complexity.

Finally, similarly to the Frank-Wolfe algorithm, the update of the current iterate is performed via a convex combination, see lines 9-10 in Algorithm 3. The choice of γ is made to minimize the value of the cost function over the line segment joining the two points.

Algorithm 3: Projected gradient method for problem (32)

Input: y, AĪ, s, α, η, w 0 , q 0 1 w (0) = w 0 , q (0) = q 0 // Initialization 2 t = 1 // Iteration index 3 repeat

// Gradient step 4 q (t-2 ⁄3) = q (t-1) + η A T Ī z (t-1) 5 w (t-2 ⁄3) = w (t-1) + η α s T z (t-1) -λ 6
where z (t-1) = y -AĪq (t-1) -s α w (t-1)

// Projection step (solved via Algorithm 4)

7 ( w (t-1 ⁄3) , q (t-1 ⁄3) ) 8 = arg min q, w q -q (t-2 ⁄3) 2 2 + ( w -w (t-2 ⁄3) ) 2 subject to αq ≤ w, -αq ≤ w // Convex update 9 w (t) = γw (t-1) + (1 -γ)w (t-1 ⁄3) 10 q (t) = γq (t-1) + (1 -γ)q (t-1 ⁄3) // Update iteration index 11 t = t + 1 12 until convergence criterion is met Output: q (t) , w (t)
Before concluding this section, we elaborate on the choice of the parameters α. It can be seen (see [START_REF] Bertsekas | Nonlinear Programming[END_REF]Section 2.3.1]), that the speed of convergence of Algorithm 3 to the global minimizer of (32) is a direct function of the conditioning of the Hessian of the cost function5

H = α -2 s 2 2 α -1 s T AĪ α -1 A T Ī s A T Ī AĪ . (33)
A proper choice of α may thus help improving the conditioning of this matrix and, consequently, enhance the speed of convergence of the overall procedure. In this paper, we suggest the following rule of thumb:

α = s 2 if s = 0 m 1 otherwise. ( 34
)
This choice is motivated by the ideal case where the columns of AĪ and s are orthogonal. In this case, the Hessian matrix in (33) is diagonal and (because the columns of AĪ have unitnorm) we have:

σ min (H) σ max (H) = min α -1 s 2 , α s -1 2 . ( 35
)
Setting α as in (34) ensures that the conditioning of H is equal to 1 in this case. Although in practice AĪ and s are commonly not orthogonal, we noticed in our numerical experiments that (34) leads to much better speeds of convergence than the trivial choice α = 1.

Algorithm 4: Projection of ( w, q) onto the feasible set of (32). Input: q, w, α

// Initialization 1 w (0) = w 2 Q (0) {i : α|q(i)| ≥ w (0) } 3 t = 0 // Iteration index // Main recursion 4 repeat 5 t = t + 1 6 w (t) = α 2 α 2 + card(Q (t-1) )   w + i∈Q (t-1) |q(i)| α   7 Q (t) i : α|q(i)| ≥ w (t) 8 until Q (t) = Q (t-1)
// Final assignment 9 if w (t) ≤ 0 then 10 w proj = 0 and q proj = 0 11 else 12

w proj = w (t) 13 q proj (i) = q(i) if i / ∈ Q (t) α -1 sign(q(i)) w proj otherwise 14 end
Output: w proj , q proj

V. NUMERICAL EXPERIMENTS

In this section, we report several simulation results illustrating the relevance of the proposed squeezing methodology. In Section V-A, we study the effectiveness of the proposed strategy as a function of the "quality" of the safe sphere used in its construction. In Sections V-B and V-C, we demonstrate that safe squeezing can be advocated to accelerate the resolution of problem (2 -P ∞ λ ): in Section V-B, we investigate the computational savings allowed by "dynamic squeezing"; in Section V-C, we show that for a given computational budget, safe squeezing enables to reach better convergence properties.

In our simulations, new realizations of the dictionary A and observation vector y are drawn for each trial according to the following distributions. The observation vector is generated according to a standard normal distribution. The dictionary obeys one of the following distribution: i) the entries are i.i.d. realizations of a centered Gaussian; ii) the entries are i.i.d. realizations of a uniform law on [0, 1]; iii) the rows are randomly-sampled rows of a DCT matrix; iv) the rows are shifted versions of a Gaussian curve. In the following, these four options will be respectively referred to as "Gaussian", "Uniform", "DCT" and "Toeplitz".

A. Effectiveness of the safe sphere squeezing test

In this section, we assess the effectiveness of the proposed squeezing test. More specifically, we evaluate the proportion of saturated entries of x which can be identified by the proposed procedure as a function of the "quality" of the safe sphere. Fig. 1: Percentage of saturated entries identified with the safe sphere squeezing test as a function of the radius for (a) the Gaussian dictionary, 3 values of λ/λ max using the GAP sphere (thick line) and the ST1 sphere (dashed line) (b) the four dictionaries and λ/λ max = 0.2.

Fig. 1 presents our results. They have been obtained by repeating the experiment described below 50 times. For each simulation trial, we draw a realization of y ∈ R 200 and A ∈ R 200×300 according to the distributions described at the beginning of Section V. We then compute a "high-accuracy" primal-dual solution (x a , u a ) of ( 2

-P ∞ λ )-(11-D ∞ λ ).
x a is evaluated numerically by solving (2 -P ∞ λ ) with Algorithm 3 and a dual gap of 10 -14 as stopping criterion. u a is obtained by dual scaling of y -Ax a (see [START_REF] Xiang | Fast Lasso screening tests based on correlations[END_REF]). This high-accuracy couple is used to construct the following safe spheres:

c ST1 = y (36) r ST1 = r 0 + y -u a 2 (37)
and

c GAP = u a (38) r GAP = r 0 + 2 gap( x a ∞ , x a , u a ) ( 39
)
where gap is defined in [START_REF] Bertsekas | Nonlinear Programming[END_REF] and r 0 ≥ 0 is some parameter. If r 0 = 0, (36)-( 37) and ( 38)-( 39) respectively correspond to the parameters 6 of the ST1 and GAP safe spheres presented in Section III-D. By construction, these expressions thus lead to safe spheres for any value of r 0 ≥ 0. In our experiments, we consider (36)-( 39) with different r 0 ≥ 0 to simulate different qualities of safe spheres. More specifically, if r 0 = 0 and (x a , u a ) = (x , u ), r ST1 (resp. r GAP ) is the smallest possible radius for an ST1 (resp. GAP) safe sphere. For each value of r 0 and each simulation trial, we apply the proposed safe squeezing test [START_REF] Vural | A new outlier detection method based on anti-sparse representations[END_REF] with the spheres defined in (36)-(39). Fig. 1 represents the proportion of saturated entries of x detected by our squeezing test as a function of parameter r 0 . The results are averaged over 50 realizations for each value of r 0 . Fig. 1a shows the percentage of saturated entries detected with a Gaussian dictionary and three different values of λ/λ max . Although we restrict our attention to the Gaussian dictionary, the following comments are valid for the other dictionaries as well. When r 0 = 0 (i.e., the radius of the safe sphere is equal to r GAP or r ST1 ), one observes that all the saturated entries are detected with the GAP sphere. However, only a subset of these entries are detected with the ST1 sphere (up to 60% in the best case). This result was expected since the high accuracy estimate of u leads to u a u and r GAP is therefore close to zero when r 0 = 0 (see Section III-E). On the other side, the radius of the ST1 sphere is usually bounded away from zero even when u a = u and r 0 = 0.

For the two tests, the number of saturated entries detected decreases as the volume of the safe sphere (parametrized by r 0 here) increases. The performance also seems to depend on λ: higher values lead to better detection results.

Fig. 1b shows the percentage of saturated entries detected for the four dictionaries when λ/λ max = 0.2. We observe that some dictionaries seem to be more prone to squeezing. For example, up to 30% additional saturated entries are detected with the DCT dictionary as compared to the Toeplitz dictionary when r 0 0.3.

B. Complexity savings

In this section, we evaluate the overall computational gain induced by the proposed method in the resolution of problem (2 -P ∞ λ ). More precisely, the total number of operations 7 carried out by several algorithms (with and without squeezing) to solve (2 -P ∞ λ ) to some accuracy are compared. In the case where a sequential implementation of the optimization procedures is considered, this figure of merit can be directly related to CPU time needed to solve (2 -P ∞ λ ); if parallel implementation is envisaged, the number of operations can be interpreted in terms of "energy consumption".

Four numerical procedures are assessed: (i) FITRA, an accelerated proximal gradient algorithm proposed in [START_REF] Studer | PAR-aware large-scale multi-user MIMO-OFDM downlink[END_REF]; (ii) the Frank-Wolfe procedure described in Algorithm 2 with no squeezing (I = ∅); (iii) the projected gradient presented in Algorithm 3 with dynamic squeezing (see Algorithm 1); (iv) the Frank-Wolfe in Algorithm 2 with dynamic screening. Procedures (ii), (iii) and (iv) are respectively denoted as "FW", "PGs", "FWs" in the sequel. Finally, empirical evidence suggests that the GAP sphere test is more effective than ST1, and we concentrate therefore on the former in our experiments.

Since this experiment aims at demonstrating the advantage of resorting to safe squeezing, we only compare procedures belonging to the same family of algorithms and therefore with similar rate of convergence. In such a setting, the induced reduction in terms of number of operations is most likely only due to safe squeezing. In particular, we compare FW against FWs and FITRA against PGs since the two couples of procedures belong to the families of Frank-Wolfe and proximal/projection gradient methods, respectively.

The numbers of operations are computed for a decreasing sequence of penalization parameters {λ j } p j=1 . Moreover, we consider a "warm-start" initialization: for each λ j , the optimization procedure is initialized with the solution obtained for λ j-1 . For λ 1 , all algorithms are initialized with the zero vector 0 ∈ R n . For each run, a stopping criterion in terms of dual gap 7 We restrict our attention to multiplications since they entail a much higher computational burden than additions in floating-point arithmetic. is used. To obtain results with the same order of magnitude 8 , two values of the dual gap are used: 10 -7 for FITRA and PGs and 10 -4 for the two Frank-Wolfe procedures (FW and FWs). All results are averaged over 50 simulation trials. For each simulation trial, we draw a realization of y ∈ R 100 and A ∈ R 100×150 according to the distributions described at the beginning of Section V. Fig. 2 represents the averaged number of operations needed by FITRA, PGs, FW and FWs to achieve the required convergence accuracy. Each figure corresponds to a different dictionary (Gaussian, Uniform, DCT and Toeplitz). The curves are plotted as a function of -log 10 (λ/λ max ). We note that, to allow a fair comparison between the different algorithms, we took particular care in the counting of the operations. In particular, each operation which can be reused latter on in an optimization procedure is not counted twice in our simulation.

In all cases, one observes that PGs (straight blue line) achieves convergence with fewer operations than FITRA (dashed blue line). As an example, up to 10 4 times more operations are required by FITRA when -log 10 (λ/λ max ) = 0.1 and A is a DCT dictionary (third row). However, such a "computational" saving tends to decrease as -log 10 (λ/λ max ) increases. Such a behavior is in good accordance with typical results of the screening literature where the performance degrades for low values of the penalization parameter.

As far as our simulation setup is concerned, these results show the advantage of considering safe squeezing to reduce the computational burden of computing antisparse representations. The observed reduction of complexity may be explained by analyzing the cost of one iteration of the four algorithms. Indeed, since the squeezing tests are done at few additional costs, the implementations of the four considered algorithms results in a cost per iteration of O(mn (t) ) where n (t) = n -card(I (t) ) and I (t) is the number of saturated entries detected at iteration t. While n (t) is decreasing along iterations for PGs and FWs, it remains constant and equal to n for FITRA and FW. The two mains consequences are (i) the computational cost of one iteration is smaller and (ii) the induced dimensionality reduction may help improving the conditioning of the optimization problem, resulting potentially in fewer iterations needed to reach convergence.

C. Benchmark profiles

As a final assessment of the performance of the proposed squeezing test, we make use of the Dolan-Moré performance profiles [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF]. For each family of dictionaries, we generate a test set consisting of p = 50 random dictionaries A ∈ R 100×150 and observations y ∈ R 100 according to the distributions described at the beginning of Section V. Each problem is then solved by using the four algorithms mentioned in the previous section, with a budget of 10 8 operations. 9 For each run j ∈ {1 . . . p} and solvers s ∈ {FITRA, PGs, FW, FWs}, we define d j,s dual gap achieved by solver s at run j.

(40)

In this experiment, we compare the performance profile

ρ s (τ ) 100 p card({j : d j,s ≤ τ }) ∀τ (41)
of each algorithm s. In other words, for all possible values of the dual gap τ , ρ s (τ ) is the (empirical) probability that solver s reaches a dual gap no greater than τ .

Figure 3 presents the performance profiles of the four algorithms with different dictionaries and values of the ratio λ/λ max . As far as our simulation setup is concerned, these results show that for a given budget of operations and any targeted dual gap smaller that 10 -4 , resorting to PGs (resp. FWs) allows to solve 10 more problems than FITRA (resp. FW). Moreover, PGs outperforms the other procedures: in all scenarios, at least 45% of the runs reach a dual gap with order of magnitude equal to machine precision (τ = 10 -16 ).

As detailed in the last paragraph of Section V-B, resorting to safe squeezing may reduce the computational burden of one iteration of an optimization procedure and help improving the conditioning of the optimization problem through an easy rescaling. With this conclusion in mind, the performance improvements observed in this experiment can be explained by the combination of two following factors. First, improving the conditioning of the problem results potentially in fewer iterations needed to reach a given dual gap. Second, reducing the computational cost of one iteration allows to perform more iterations for the same budget of operations, reaching possibly a solution with a lower dual gap.

VI. CONCLUSIONS

In this paper, we proposed a new methodology, dubbed safe squeezing, to accelerate algorithms that compute a so-called antisparse representation. The proposed procedure consists in a set of tests that detect saturated entries in the solution of the optimization problem. These tests rely on the notion of safe region, that is a set that contains the solution of the dual problem. Benefiting from the recent developments in the screening literature, we showed that two existing safe spheres can be extended in the context of safe squeezing. The contribution of the saturated entries is then compacted into one single vector, allowing for a potentially dramatic dimensionality reduction. Computing an antisparse representation becomes equivalent to solving a lower dimensional quadratic problem. When interleaved with the iterations of an optimization procedure, the approach permits to dynamically detect new saturated entries all along the optimization process. Numerical simulations showed that the proposed safe squeezing methodology leads to significant computational gains in several numerical setups. [START_REF] Foucart | A Mathematical Introduction to Compressive Sensing[END_REF] In this section, we prove that "λ < A T y 1 " is indeed a necessary and sufficient condition for the solution of (2 -P ∞ λ ) to be nonzero. A direct application of the Fermat's rule [START_REF] Bauschke | Subdifferentiability[END_REF]Th. 16.2] ensures that x is a minimizer of (2

APPENDIX A PROOF OF THE RESULT IN

-P ∞ λ ) if and only if λ -1 A T (Ax -y) ∈ ∂ x ∞ (42)
where ∂ x ∞ denotes the sub-differential of the ∞ -norm evaluated at x . Danskin's theorem [25, Th. B.25] leads to

∂ x ∞ d ∈ R n : d 1 ≤ 1 and x T d = x ∞ (43)
for all x ∈ R n . Result (5) follows by particularizing (42)-(43) to x = 0 n .

APPENDIX B PROOF OF THE RESULTS IN SECTION III

In this section, we provide a proof of Theorem 1. We also show that the quantity u defined in (11-D ∞ λ ) is independent of the set I.

In our derivations, we focus on the following generic optimization problem:

min z∈R m ,q∈R p w≥0 1 2 y -z 2 2 + λw s.t.      q ≤ w -q ≤ w z = Bq + sw. (44)
Conclusions about problems (2 -P ∞ λ ) and (7 -P sq λ ) are then drawn by considering the following substitutions:

Problem (2 -P ∞ λ ): B = A, s = 0 m (45) Problem (7 -P sq λ ): B = AĪ, s = i∈I+ a i - i∈I- a . ( 46
)
We organize the presentation of the results as follows. In Section B-A we derive the Lagrangian dual problem of (44). In Section B-B, we emphasize some properties of the solutions of the dual problem. The latter are then exploited in Section B-C to derive (the first part of) the result stated in Theorem 1. In Section B-D, we elaborate on the uniqueness of the maximizer of the dual problem and its independence with respect to the set I. Finally, in Section B-E we prove the last part of Theorem 1.

A. Dual problem of (44) Let v + ∈ R m and v -∈ R m (resp. u ∈ R n ) denote the dual variables associated to the inequality (resp. equality) constraint(s) in (44). The Lagrangian function associated to this problem can then be rewritten as

L (q, z, w, u, v + , v -) = 1 2 y -z 2 2 + u T z + w λ -s T u -1 T n v + -1 T n v - + q T -B T u + v + -v -.
Moreover, the corresponding dual function is defined as

D(u, v + , v -) min q,z,w≥0 L (q, z, w, u, v + , v -).
We note that since v + ≥ 0 and v -≥ 0, we have D(u, v + , v -) = -∞ as soon as the following conditions are not satisfied:

λ -s T u -1 T n v + -1 T n v -≥ 0, -B T u + v + -v -= 0 p . (47)
Hence, defining

D {(u, v + , v -) : v + ≥ 0, v -≥ 0
and (47) holds}, the dual function D(u, v + , v -) can be written as

D(u, v + , v -) = 1 2 y 2 2 -1 2 y -u 2 2 if (u, v + , v -) ∈ D -∞ otherwise.
The dual problem associated to (44) thus reads

max (u,v+,v-)∈D 1 2 y 2 2 -1 2 y -u 2 2 . (48)
From [25, Proposition 5.2.1], 11 there is no duality gap between (44) and (48); moreover the dual problem (48) has (at least) one maximizer. In fact, we will show in Section B-D below that (48) admits one unique maximizer, denoted (u , v + , v -) hereafter.

Since card(Q (t-1) ) must be nonnegative and card(Q (0) ) ≤ q, we thus have necessarily from (73) that t ≤ q. This proves the convergence of Algorithm 4 in at most q steps. In practice, Algorithm 4 is implemented as follows. Vector q is first sorted by absolute value and its cumulative sum is precomputed. This requires a complexity of O(q log q). Then, each of the (at most) q repetitions of lines 5-6-7 scales as O(log(card(Q (t) )) since the complexity is concentrated in the identification of the components of Q (t) verifying "α|q(i)| ≥ w (t+1) ". Inasmuch as card(Q (t) ) ≤ q, the overall complexity of the proposed procedure is therefore O(q log q).

B. Optimality conditions for (70)

In this section, we derive necessary and sufficient optimality conditions for problem (70). These conditions will be exploited in the next section to show the correctness of Algorithm 4.

We first note that sign(q ) = sign(q). (74) Indeed, q is also the minimizer of min q qq 2 2 subject to ± αq ≤ w ,

which admits the following simple analytical expression:

∀i : q (i) = sign(q(i)) min(|q(i)|, w ).

This shows (74).

In the rest of the appendix, we will thus focus, without loss of generality, on the case where q(i) ≥ 0 ∀i.

(77)

We note that this setup can always be obtained by some simple changes of variables. In this case, we have from (74) that q ≥ 0 q and problem (70) can equivalently be expressed as min q ≥0q, w ≥0

( w -w) 2 + qq 2 2 s.t. αq ≤ w . (78)

Since ( 78) is convex with linear inequality constraints, any minimizer ( w , q ) must satisfy the following set of necessary and sufficient optimality conditions [25, Propositions 5.2.1 and 5.1.5]:

αq ≤ w , (79)
and ∃v ≥ 0 q such that ∀i : 0 = v (i)(αq (i) -w ),

( w , q ) = arg min w ≥0,q ≥0q

( w -w) 2 + qq 2 2 (81)

+ q i=1
v (i)(αq (i) -w ) .

In the rest of this section, we derive other optimality conditions, equivalent to (79)-( 81), which are more amenable to prove the correctness of Algorithm 4.

We first note that (81) is a separable nonnegative leastsquare problem and its minimizer admits the following closedform expression:

w = max 0, w + q i=1 v (i) (82)
q (i) = max 0, q(i) -αv (i) .

We then distinguish between the cases " w = 0" and " w > 0". From (79), we see that the case " w = 0" may occur only if q = 0 q . Using (82)-(83), this leads to ( w , q ) = (0, 0 q ) ⇐⇒ ∃v ≥ 0 q s.t.

     w + q i=1 v (i) ≤ 0, α -1 q(i) ≤ v (i),
or more compactly, ( w , q ) = (0, 0 q ) ⇐⇒ w + α -1 q i=1 q(i) ≤ 0.

(84)

We next consider the case " w > 0". First, this assumption together with (82) leads to

w = w + q i=1 v (i). (85)
Moreover, from (80) we have that

v (i) = 0 ∀i / ∈ Q , (86)
where

Q {i : αq (i) = w }. (87)
(85) thus further simplifies to

w = w + i∈Q v (i). (88)
On the other hand, we have (by definition of Q ) that

q (i) = α -1 w ∀i ∈ Q (89)
so that, invoking again our hypothesis " w > 0", we obtain from (83)

q (i) = q(i) -αv (i) ∀i ∈ Q . ( 90
)
Finally, combining (88), ( 89) and (90) leads to

w = α 2 α 2 + card(Q ) w + α -1 i∈Q q(i) . (91)
As a consequence ( w , q ) with w > 0 is a solution of (78) if and only if ∃v ≥ 0 q such that (88)-(91) are satisfied.

C. Optimality of ( w proj , q proj )

We finally prove the correctness of Algorithm 4. We assume that the algorithm terminates its main recursion (steps 4-8) after t iterations (that is Q (t) = Q (t-1) ) and show that the assignments defined in steps 9-14 are such that ( w proj , q proj ) = ( w , q ). First, if w (t) ≤ 0, we must have (by definition of Q (t) ) that

Q (t) = {1, . . . , q} (92)
since q ≥ 0 q . Combining this equality with the definition of w (t) in line 6 of Algorithm 4 leads to

w (t) = α 2 α 2 + q w + α -1 q i=1 q(i) . (93)
Since w (t) ≤ 0, the last equality implies

w + α -1 q i=1 q(i) ≤ 0. (94)
In view of (84), we finally have w (t) ≤ 0 ⇒ ( w , q ) = (0, 0 q ). (95)

This corresponds to the assignment in line 9-10 of Algorithm 4. If w (t) > 0, it can easily be verified that the optimality conditions (88)-(91) are verified for q = q proj w = w proj v (i) = 0 if αq proj (i) < w proj α -2 (αq(i) -w proj ) otherwise, where w proj and q proj are defined in lines 12-13 of Algorithm 4. This concludes the proof.

Fig. 2 :

 2 Fig. 2: Number of operations required to achieve convergence for different values of λ/λ max and four dictionaries: Uniform (first row), Gaussian (second row), DCT (third row) and Toeplitz (fourth row).

Fig. 3 :

 3 Fig. 3: Performance profiles obtained with the Uniform, Gaussian, DCT and Toeplitz dictionaries. First row: λ/λ max = 0.3; second row: λ/λ max = 0.8.

If x is not unique, I must be understood as the set of indices of the components saturating for all the solutions of (2 -P ∞ λ ).

A related result is stated in [24, Sec. 4.1].

More particularly, this is the case with probability 1 when the columns of A are drawn according to any distribution that admits a density with respect to the Lebesgue measure.

We always have u = 0m ∈ U I since λ > 0 by hypothesis.

That is the ratio between σmax(H) and σ min (H).

Computed with the primal-dual feasible couple (xa, ua).

Frank-Wolfe procedures are indeed known to converge much more slowly than proximal/projection gradient methods.

The algorithm stops as soon as

[START_REF] Bonnefoy | Dynamic Screening: Accelerating First-Order Algorithms for the Lasso and Group-Lasso[END_REF] multiplications have been carried out.[START_REF] Herzet | Gather and conquer: Region-based strategies to accelerate safe screening tests[END_REF] that is, reach a solution with the desired dual gap.

More specifically, (44) is a feasible convex problem with linear constraints and the minimum value of the cost function is upper bounded.

This property holds true in the more generally setup where x is not unique although the proof is slightly more involved.

In particular, q should not be confused with the minimizer of problem (7 -P sq λ ).

Part of this work has been funded thanks to the Becose ANR project no. ANR-15-CE23-0021 and the Labex CominLabs.

The research presented in this paper is reproducible. Code and data are available at https://gitlab.inria.fr/celvira/safe-squeezing.

B. Connection between u ,v + and v - In this section, we emphasize that v + and v -are strongly connected to u . More specifically, we have:

This connection can be shown as follows. First, since there is no duality gap (see Section B-A), any couple of primaldual optimal variables, say (w , q , z ) and (u , v + , v -), must verify the well-known Karush-Kuhn-Tucker optimality conditions [START_REF] Bertsekas | Nonlinear Programming[END_REF]Proposition 5.1.5]. In particular, we have

More precisely, (50a) corresponds to first-order optimality condition while (50b) and (50c) partially describe the set of dual feasible variables (see (47)). One also recognizes in (50d) and (50e) the complementary slackness conditions of the problem. Using these conditions, we have

> 0 (52

and

Therefore, combining (50b), (50c) and (51), we obtain the following implications:

Expression (49) corresponds to a compact reformulation of (56), (57) and (58).

C. Necessary optimality condition

We now exploit (49) together with conditions (50d)-(50e) to prove the necessary optimality condition (13) stated in Theorem 1. Let us first show that u can be expressed as

where

To this end, let us define

Clearly, from ( 48) and (49) we have that

Moreover, combining (47) and (49) leads to

Plugging this definition into (62), we obtain that u is defined by ( 59)-(60). Moreover, v + and v -are univocally specified from u via (49). We finally show that the following condition holds between the optimal primal and dual variables:

Condition ( 13) then corresponds to a particular instantiation of (63) when B and s are defined as in (46). Relation (63) can be obtained as follows. Taking the difference between conditions (50d)-(50e), we first have

Now, using (49) we obtain

Using (65a)-(65b) and the fact that b T i u = b T i u sign(b T i u ), we finally end up with the desired result.

D. Uniqueness of u and independence with respect to I

Problem (59) admits a unique minimizer since the cost function is continuous, coercive and strictly concave, see [START_REF] Bertsekas | Nonlinear Programming[END_REF]Propositions A.8 & B.10]. This shows the uniqueness of u . We note that the uniqueness of (u , v + , v -) claimed in Section B-A then straightforwardly follows from property (49).

Setting B and s as in (46), we now show the independence of u with respect to the definition of I ⊆ I . We note that in this case, problem (44) is equivalent to (7 -P sq λ ). From optimality condition (50a), we have

Now, we have by definition that

Moreover, because problems (2 -P ∞ λ ) and (7 -P sq λ ) are equivalent (see Section III), we have from ( 9)-( 10) that

where x is the unique (see Section III-A) minimizer of (2 -P ∞ λ ). Combining the last three equalities, we obtain that u does not depend on the choice of I as long as I ⊆ I . 12

E. Minimum number of saturated components

In this section we prove the last part of Theorem 1. Since we showed in Section B-D that u is independent of the choice of I ⊆ I , we focus on the case where I = ∅ and B = A, s = 0 m .

We first note that our working hypothesis (4) prevents us from having u = 0 m . Indeed, considering the optimality conditions of problem (59), it can be seen that

Now, as pointed out in Section III-A, hypothesis (4) implies that y / ∈ ker(A T ) and thus necessarily u = 0 m . Since u = 0 m , we observe that "a T i u = 0" can only occur for at most m -1 a i 's when kruskal(A) = m. Indeed, kruskal(A) = m implies that any subset of m columns of A are linearly independent. On the other hand, the vector space P {a ∈ R m : a T u = 0} is a (m-1)-dimensional subspace of R m . Hence, if more than m -1 columns of A belong to P, they must be linearly dependent, which is in contradiction with kruskal(A) = m. This proves the last part of Theorem 1.

APPENDIX C CORRECTNESS OF ALGORITHM 4

In this appendix, we prove that the procedure described in Algorithm 4 solves the following projection problem min q ,w (w -w) 2 + qq 2 2 subject to αq ≤ w , -αq ≤ w (70) in a finite number of steps. With a slight abuse of notations, 13 we denote hereafter the unique minimizer of this problem as ( w , q ). The dimension of q will be denoted q hereafter. The appendix is organized as follows. In Section C-A, we prove the convergence of Algorithm 4 in a number of steps no greater than q and discuss its computational complexity. In Section C-B we elaborate on necessary and sufficient conditions for ( w , q ) to be a minimizer of (70). Finally, we exploit these optimality conditions in Section C-C to show that the output of Algorithm 4 corresponds to the unique minimizer of (70).

A. Convergence of Algorithm 4 in a finite number of steps

We first note that, by construction, the sequence w (t) is non-decreasing. As a matter of fact, by definition of Q (t) we must have

Moreover, if equality occurs in (71), then the stopping criterion of Algorithm 4 is satisfied and the procedure terminates. Algorithm 4 carries out t iterations provided that

and, in particular, card(Q (t-1) ) < card(Q (t-2) ) < . . . < card(Q (0) ). (73)