Student Member, IEEE Cássio F Dantas
email: cassio.fraga-dantas@inria.fr

Fellow, IEEE Rémi Gribonval

Stable safe screening and structured dictionaries for faster ℓ 1 regularization

In this paper, we propose a way to combine two acceleration techniques for the ℓ1-regularized least squares problem: safe screening tests, which allow to eliminate useless dictionary atoms; and the use of fast structured approximations of the dictionary matrix. To do so, we introduce a new family of screening tests, termed stable screening, which can cope with approximation errors on the dictionary atoms while keeping the safety of the test (i.e. zero risk of rejecting atoms belonging to the solution support). Some of the main existing screening tests are extended to this new framework. The proposed algorithm consists in using a coarser (but faster) approximation of the dictionary at the initial iterations and then switching to better approximations until eventually adopting the original dictionary. A systematic switching criterion based on the duality gap saturation and the screening ratio is derived. Simulation results show significant reductions in both computational complexity and execution times for a wide range of tested scenarios.

I. INTRODUCTION

Sparsity-constrained linear regression has found numerous applications in signal processing and machine learning, tackling under-determined inverse problems. These appear in many forms such as image or audio inpainting [START_REF] Adler | Audio inpainting[END_REF], source localization [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF] or spike deconvolution [START_REF] Dossal | Sparse spike deconvolution with minimum scale[END_REF], to cite only a few.

There are many computational approaches to estimate the sparse coefficient vector in such settings with two main families: greedy algorithms [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF], and convex optimization approaches relying on ℓ 1 -norm minimization, where the regression problem is expressed as an ℓ 1 regularized optimization problem -known as Basis Pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] or LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]which solution is computed using iterative convex optimization algorithms such as ISTA or its variants [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]- [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

For very high-dimensional problems, however, iterative algorithms to solve ℓ 1 minimization problems can become computationally prohibitive, which is why accelerating techniques are still an intense research topic. This paper demonstrates how to combine two such techniques:

1) Fast structured operators [START_REF] Rubinstein | Double sparsity: Learning sparse dictionaries for sparse signal approximation[END_REF]- [START_REF] Dantas | Learning dictionaries as a sum of Kronecker products[END_REF] which provide faster matrix-vector products (see Section III-C); 2) Safe screening tests [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF]- [START_REF] Xiang | Screening tests for lasso problems[END_REF], which safely eliminate unused explanatory variables (see Section II). This paper extends the results in [START_REF] Dantas | Dynamic screening with approximate dictionaries[END_REF] [START_REF]Faster and still safe: combining screening techniques and structured dictionaries to accelerate the lasso[END_REF]. In [START_REF] Dantas | Dynamic screening with approximate dictionaries[END_REF] we introduced a safe screening sphere test that manipulates an approximate dictionary and in [START_REF]Faster and still safe: combining screening techniques and structured dictionaries to accelerate the lasso[END_REF] the GAP Safe sphere test [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF] was extended to this new setting. Here, we revisit the previous results under a broader formalism and propose a fast algorithm for ℓ 1 -minimization which combines safe screening and (potentially multiple) fast approximate dictionaries.

A. Sparsity constrained regularization

Let N and K be respectively the dimension of the observed vector and that of the unknown coefficient vector. The observed vector is denoted y ∈ R N and modeled as y ≈ Ax 0 where x 0 ∈ R K is sparse and A = [a 1 , . . . , a K ] ∈ R N ×K . In the context of linear inverse problems, A is the so-called measurement matrix [START_REF] Foucart | A Mathematical Introduction to Compressive Sensing[END_REF], while for sparse signal representations A would be the dictionary matrix [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], and in statistics the design matrix [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. In this paper, we adopt the terminology of sparse signal representations hence A is called a dictionary and its columns, denoted a j ∈ R N , are called atoms.

The ℓ 1 -regularized least squares, referred to as Lasso or Basis Pursuit, consists in finding a sparse coefficient vector x ∈ R K , solution of the following optimization problem:

L(λ, A, y) : x ⋆ = argmin x 1 2
Axy 2 2 + λ x 1

P (x|A) (1)
where P (x|A) is called the primal objective and the parameter λ > 0 controls the trade-off between data fidelity and sparsity of the solution. We suppose that

λ ≤ λ max := A T y ∞ (2)
since otherwise x ⋆ = 0 ∈ R K is the unique solution.

B. Iterative algorithms

First-order iterative algorithms, especially proximal-gradient methods (ISTA [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], TwIST [START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF], SpaRSA [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF], Chambolle-Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]), are popular approaches for solving [START_REF] Dantas | Stable screening -Python code[END_REF].

One can use an abstract notations to represent the update step for a generic iterative algorithm for the problem L(λ, A, y): {x t+1 , α t+1 } ← p(x t , A, α t ) where x t is the current estimate of the primal variable and α t is a list of updated auxiliary scalars (e.g. the gradient step-size, and possibly a few previous estimates of the primal variable). For instance, the ISTA algorithm is given by

x t+1 ← p(x t , A, L t ) = ST λ L t x t + 1 L t A T (y -Ax t ) (3)
where ST u (z) = sign(z)(|z| -u) + denotes the softthresholding operation, which is the proximal operator associated to the ℓ 1 -norm.

C. Computational bottlenecks

The main bottleneck of existing iterative algorithms in terms of computational complexity is the cost of the required matrixvector products involving the dictionary matrix, which dominates the overall iteration cost. For example, ISTA requires two matrix-vector multiplications at each iteration (or at least one if the Gram matrix A T A is pre-computed and stored).

D. Addressing the computational bottlenecks

A popular way to address this limitation is to constrain the dictionary matrix to a certain type of structure which would allow for fast matrix-vector products. Another approach is to use screening rules to identify a set indexing zero components of the solution x ⋆ before even computing it.

1) Acceleration with structured dictionaries: Different types of structure, more or less suited to a given application, can be imagined. In this kind of work, there is always a compromise between the flexibility (generality) of the structure and the provided speedup.

In [START_REF] Rubinstein | Double sparsity: Learning sparse dictionaries for sparse signal approximation[END_REF], the authors constrain the dictionary to be the product of a fixed base dictionary that has a fast implementation (a DCT, for instance) and a column sparse matrix. This model is extended in [START_REF] Sulam | Trainlets: Dictionary learning in high dimensions[END_REF] by replacing the base dictionary with cropped wavelets, which have a certain degree of adaptability. A more general model, which has the two previous ones as particular cases, was proposed in [START_REF] Magoarou | Flexible multilayer sparse approximations of matrices and applications[END_REF]. The dictionary is constrained to be the product of several sparse matrices.

Another interesting line of research, particularly suited to multi-dimensional signals, consists in finding dictionaries formed as the Kronecker (or tensorial) product of subdictionaries of smaller size. It was proposed in [START_REF] Hawe | Separable dictionary learning[END_REF] for the Kronecker product of two sub-dictionaries aiming 2D signals (like images). The model was then extended to a sum of Kronecker products in [START_REF] Dantas | Learning dictionaries as a sum of Kronecker products[END_REF] and [START_REF] Dantas | Learning fast dictionaries for sparse representations using low-rank tensor decompositions[END_REF] (respectively for 2dimensional and higher-order tensorial data), increasing its flexibility in exchange of some additional complexity.

In practical applications, the dictionary matrix A is not necessarily structured. A possible strategy is to replace certain iterations x t+1 = p(x t , A, α t ) with x t+1 = p(x t , Ã, α t ) where à is a structured approximation of A, i.e.,

A = Ã + E (4) with approximation error E = [e 1 , . . . , e K ] ∈ R N ×K .
2) Acceleration with screening rules: For a given instance L(λ, A, y) of (1), characterized by a regularization parameter λ, a dictionary A, and an input vector y, the safe screening approach, introduced by [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF], consists in identifying and removing from the dictionary a subset of atoms which are guaranteed to have zero weight in a solution x ⋆ , before solving the problem. By removing these so-called inactive atoms a more compact and readily solved problem is obtained, with decreased matrix-vector multiplication cost, while not affecting at all the original solution which can be obtained by simply zero-padding its restricted version.

Put differently, safe screening is a feature selection technique for a given instance of (1). However, unlike other previous feature selection heuristics [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF], [START_REF] Tibshirani | Strong rules for discarding predictors in lasso-type problems[END_REF] based on correlation measures between the atoms a j and the input signal y, safe screening has zero risk of false rejections. Moreover, screening techniques are transparent to the underlying ℓ 1 solver and can be readily combined with almost any existing solvers.

Basically, three categories of screening rules can be distinguished: 1) Static; 2) Dynamic; 3) Sequential. The two first categories being mutually exclusive, but not the third one.

The earliest safe screening techniques [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF] [20] [START_REF] Xiang | Fast lasso screening tests based on correlations[END_REF], currently classified as static rules, where designed to be applied once and for all before starting the optimization. In contrast, the more recent dynamic rules, introduced in [21] and followed by [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF], are repeatedly applied during the iterative optimization algorithm leveraging its current solution estimate and gradually increasing the set of rejected atoms.

Sequential rules [START_REF] Xiang | Learning sparse representations of high dimensional data on large scale dictionaries[END_REF], [START_REF] Wang | Lasso screening rules via dual polytope projection[END_REF]- [START_REF] Liu | Safe screening with variational inequalities and its application to lasso[END_REF] exploit the fact that ( 1) is commonly solved over a grid of regularization parameters λ i and reuse the results from a previous configuration λ i-1 to improve the screening performance.

E. Contributions

First, we introduce a formalism for defining safe screening tests which are robust to approximation errors on the dictionary matrix given that an error bound is provided. The resulting tests are called stable screening tests. The stable tests proposed here are general and apply to any context in which the atoms are known up to a certain error margin. The source of this error can be manifold, but we will focus on the case where it is a side effect of manipulating structured approximations of the true dictionary matrix. The proposed framework is also general with respect to the form of the safe region (be it a sphere, dome, etc) as well as to that of the error bound. Here, again, we stick to ℓ p -balls, which allows us to exemplify the stable tests on some existing ℓ 2 -sphere tests -one of which being the state-of-the-art GAP Safe [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF]for a general ℓ q -ball error bound. Extension to dome tests is not particularly difficult within the proposed formalism.

In a second part, we exploit the proposed stable screening in a fast algorithm for ℓ 1 -minimization problems, by making use of fast structured approximations of the problem's dictionary matrix. It consists in starting the iterative optimization by using a coarser (but faster) version of the true dictionary and, as approaching convergence, finer approximations are progressively adopted until eventually the original dictionary takes over. Choosing an appropriate moment to switch to a more precise dictionary is crucial in the proposed algorithm. A robust switching criterion based on both the duality gap saturation and the screening ratio is proposed.

F. Related Work

Apart from the aforementioned structured dictionaries and safe screening tests -as well as other preceding correlationbased feature selection heuristics [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF], [START_REF] Tibshirani | Strong rules for discarding predictors in lasso-type problems[END_REF]-some related acceleration strategies for sparsity-inducing optimization problems (or more specifically for problem (1)) are worth citing.

Instead of starting from the full problem and pruning the feature set, working set techniques [START_REF] Kim | Fast active-set-type algorithms for l1-regularized linear regression[END_REF], [START_REF] Johnson | Blitz: A principled meta-algorithm for scaling sparse optimization[END_REF] start with small restricted problems and progressively include more promising features. In [START_REF] Massias | From safe screening rules to working sets for faster lasso-type solvers[END_REF], the authors combine a working set strategy with safe screening and in [START_REF] Massias | Celer: a fast solver for the lasso with dual extrapolation[END_REF] they incorporate a dual extrapolation technique to further enhance the screening performance and accelerate convergence. This idea is conceivably complementary to the techniques proposed here.

Joint screening [START_REF] Herzet | Joint screening tests for lasso[END_REF] allows to screen many atoms which lie close together in one single test, reducing the number of required tests for a given dictionary. Interestingly, the resulting tests share many similarities and mathematical connections to the stable screening tests introduced here, despite arising from an essentially different premise.

G. Outline of the Paper

This paper is organized as follows. In Section II we briefly review standard (static and dynamic) safe screening rules. Stable safe screening rules are introduced in Section III. A complete description of the proposed algorithm is given in Section IV and experimental results in Section V.

II. REMINDERS ON SAFE SCREENING

Given a set A of non-repeating integers, A [A] := [a i ] i∈A denotes a sub-matrix of A composed of the columns indexed by the elements in A. The notation extends to vectors:

x [A] := [x(i)] i∈A ,
where x(i) denotes the i-th component of x. Screening rules allow to identify such sets of reduced cardinality |A t | < K, so that the complexity of applying A [A] or its transpose is reduced compared to that of applying A, while ensuring that A t contains the support of the solution {i : x ⋆ (i) = 0}. This is achieved using the dual of (1) [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF]:

Since the optimal solution θ ⋆ of the dual problem (5) is not known, the inner products a T j θ ⋆ cannot be evaluated. Fortunately, given only A and y, it is possible to identify at a moderate computational cost a region R ⊂ R N , called safe region, which is guaranteed to contain the optimal θ ⋆ . Definition 1 (Safe region). A region R ⊂ R N is safe (with respect to the dual problem (5)) if and only if it contains the dual solution θ ⋆ , i.e. θ ⋆ ∈ R.

B. Screening an atom given a safe region

Consider an atom a j . If the inequality |a T j θ| < 1 holds for all θ ∈ R where R is a safe region, then the above analysis ensures that a j is inactive. This gives rise to a screening test.

Definition 2 (Screening test for an atom). Given a region R, a screening test for the atom a ∈ R N is given by:

µ(a|R) = sup θ∈R |a T θ| (7)
A sufficient condition for an atom a j to be inactive can be expressed as follows: if R is safe then

µ(a j |R) < 1 =⇒ x ⋆ (j) = 0.
In practice, for each atom a j , computing the test µ(a j |R) allows to eliminate or not the atom. Formally, given a safe region, the atoms can be partitioned into a preserved set A = {j ∈ {1, . . . , K} : µ(a j |R) ≥ 1} and its complement, the rejection set A c , that gathers the indices of the identified inactive atoms.

In practice, safe regions have simple parameterized shapes which parameters need to be identified with moderate computations from the only knowledge of A, y, and possibly the current iterate x t of an iterative algorithm addressing [START_REF] Dantas | Stable screening -Python code[END_REF].

The two most common shapes of safe regions in the literature are spheres and domes (i.e. intersection between a sphere and one or more half-planes).

This paper focuses on sphere tests. Assume we are given a safe region R which is a closed ℓ p -ball with center c and radius R, denoted B p (c, R) = {z : zc p ≤ R}. The screening test for this region has a closed form (a proof is given in Appendix A)

µ(a|B p (c, R)) = |a T c| + R a p * (8)
where • p * denotes the dual norm associated to the p-norm, with 1 p + 1 p * = 1. For simplicity, we omit subscripts for the ℓ 2 -ball (p = p * = 2) denoted B(c, R).

C. Construction of a safe spherical region

A safe region should be as small as possible (to maximize the screening effect) while requiring as little computational overhead as possible. In light of (8), as a p * can be precomputed for each atom, the computational overhead of a screening test with a spherical region is governed by the cost of computing the radius R and the inner products a T c for all atoms that have not been screened out yet. This calls for techniques where the choice of R and c allows to reuse either quantities that have already been computed along previous iterations of the optimization algorithm.

The construction of the first (static) safe region, obtained in [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF], follows from the simple observation that the solution θ ⋆ of the dual problem [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] is the Euclidean projection of y/λ on the feasible set ∆ A . As a result, if some feasible point θ F ∈ ∆ A is known, then θ ⋆ can't be further away from y/λ than θ F in the ℓ 2 sense, i.e. θ ⋆y/λ 2 ≤ θ Fy/λ 2 ∀θ F ∈ ∆ A .

This leads to an ℓ 2 -spherical safe region (p = p * = 2),

R = B(c = y/λ, R = θ F -y/λ 2 ), (9)
which full specification requires a feasible point θ F ∈ ∆ A .

To generate a feasible point θ F ∈ ∆ A , one could compute the Euclidean projection of any given point z ∈ R N over the closed convex set ∆ A . As this is too computationally demanding, simpler scaling approaches are preferred: given z ∈ R N we determine a scaling factor α so that Θ(z|A) := αz ∈ ∆ A .

By definition of ∆ A , if A T z ∞ ≤ 1 no scaling is needed, otherwise one only needs to divide z by A T z ∞ . This yields

Θ(z|A) := z max (1, A T z ∞ ) . ( 10
)
While the scaling in 10 selects a feasible point to build a safe sphere [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] given an arbitrary z, it neglects the fact that one would like to minimize the radius of this sphere to maximize the effect of screening. This leads to the dual scaling [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF] formalized as follows (a proof is given in [START_REF] Bonnefoy | Dynamic screening: Accelerating first-order algorithms for the lasso and grouplasso[END_REF]Lemma 8])

Proposition 1. Denoting proj z y λ = y T z λ z 2 2
z the Euclidean projection of y/λ onto the direction of z, we have

Θ proj z y λ | A = argmin z ′ =αz z ′ -y/λ 2 s.t. A T z ′ ∞ ≤ 1
i.e. among all points which are both dual feasible and proportional to z, Θ proj z y λ | A is the one closest to y/λ.

To obtain dynamic spherical regions, we need dual feasible points that change over iterations (θ F = θ t at iteration t) and leverage the current state of the optimization algorithm. Following [START_REF] Bonnefoy | Dynamic screening: Accelerating first-order algorithms for the lasso and grouplasso[END_REF], this can be achieved by defining θ t proportional to the current residual ρ t := y -Ax t , since the dual solution θ ⋆ is proportional to y -Ax ⋆ . Taking z = ρ t and

θ t = Θ proj ρ t y λ | A yields θ t = y T ρ t λ ρ t 2 2 αt -αt ρ t , with α t = ( A T ρ t ∞ ) -1 (11)
where The above principles yield existing safe spherical regions. 1) Static Safe sphere [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF]: With z = y/λ in (10) one gets θ F = Θ(y/λ|A) = y/λ max , and (9

) reads R = B(c, R) with Static Safe: c = y/λ, R = | 1 λmax -1 λ | • y 2 (12)
where we recall that λ max is defined in (2).

2) Dynamic Safe sphere [START_REF] Bonnefoy | Dynamic screening: Accelerating first-order algorithms for the lasso and grouplasso[END_REF]: With the dual feasible point defined in [START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF], [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]

yields R = B(c, R) with Dynamic Safe: c = y/λ, R = θ t -y/λ 2 . ( 13
)
3) GAP Safe sphere [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF]: Given any primal-dual feasible pair (x t , θ t ) ∈ R K × ∆ A (say the solution estimations at iteration t) and denoting

G(x t , θ t |A) := P (x t |A) -D(θ t ) (14)
the corresponding duality gap, the GAP Safe sphere is

R = B(c, R) with GAP Safe: c = θ t , R = 1 λ 2G(x t , θ t |A). ( 15
)
This region is provably safe for any dual feasible point θ t ∈ ∆ A [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF]Theorem 2]. The authors, following [START_REF] Bonnefoy | Dynamic screening: Accelerating first-order algorithms for the lasso and grouplasso[END_REF], propose to use θ t defined by [START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF]. Note that the cited safe regions reuse quantities computed along the iterations of classical proximal algorithms such as the residual ρ t , the inner products A T ρ t , or the duality gap (which typically serves as a stopping criterion).

III. STABLE SAFE SCREENING

When iterating with à instead of A, the computational overhead of accessing inner products a T c is no longer moderate, hence the screening techniques reviewed in Section II cannot be directly applied in this context. A possibility would of course be to recompute residuals / inner products / duality gaps associated to A (instead of Ã) to implement the above screening tests, but the whole purpose of re-using quantities already computed with a structured (computationally efficient) matrix à would be lost.

In this section, we propose an alternative solution by deriving screening rules which remain safe w.r.t the original problem (1) even when manipulating an approximate version of the dictionary. These stable screening rules only requires some knowledge on the magnitude of certain approximation errors between à and A.

A. Screening a zone given a safe region

In conventional screening tests (Definition 2) the dual solution θ ⋆ is assumed to belong to a known safe region R. This allows to screen (or not) an atom a that is perfectly known, (re)using computations of inner products with this atom.

When iterating with Ã, efficient screening tests are only entitled to reuse inner products with approximated atoms, ÃT c. Yet, they should be able to screen certain atoms a, using the only knowledge that these are "close" to the corresponding approximated atoms ã.

One way to capture this knowledge is to consider S ⊂ R N some neighborhood of a, assumed to be known. We will call S a zone. This gives rise to screening tests for zones. Definition 3 (Screening test for a zone). Given a safe region R, a screening test for the zone S ⊂ R N is given by:

µ(S|R) := sup a∈S µ(a|R) = sup a∈S sup θ∈R |a T θ| (16)
A sufficient condition for an atom a j to be inactive can be expressed as follows: if R is safe and if a j ∈ S, then

µ(S|R) < 1 =⇒ x ⋆ (j) = 0.

B. Stable screening tests

Given an approximate dictionary Ã, error bounds ǫ j such that ãja j q ≤ ǫ j , 1 ≤ j ≤ K, can be pre-computed for some choice of q ∈ [1, ∞], yielding spherical zones S j := B q (ã j , ǫ j ) known to contain the atoms: a j ∈ S j . Using (8), a spherical zone

S := B q (ã, ǫ) (17)
gives rise to the following estimate:

µ(S|B p (c, R)) ≤ |ã T c| + ǫ c q * + R ã p * + CRǫ. ( 18
)
where

C = C p,q := N (1-1 p -1 q ) + and q * ∈ [1, ∞] is such that 1 q + 1 q * = 1. See proof in Appendix B.
In practice, the norms a j = a j p * can be pre-computed and stored with negligible overhead, yielding restricted spherical zones

S ′ := S ∩ {a : a p * = a}. ( 19
)
For such zones we get

µ(S ′ |B p (c, R)) = |ã T c| + ǫ c q * + Ra. ( 20
)
Lemma 1. Consider a dictionary A and R = B p (c, R) a safe region with respect to the dual problem [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]. Let à be an approximate dictionary, 1 ≤ j ≤ K, and S j and S ′ j defined by ( 17)-( 20) with ã = ãj , ǫ ≥ ãja j q , and a = a j p * .

• If µ(S j |B p (c, R)) < 1 then x ⋆ j = 0; • If µ(S ′ j |B p (c, R)) < 1 then x ⋆ j = 0. Remark 1.
In [START_REF] Herzet | Joint screening tests for lasso[END_REF], Definition 3 is used in a different context. No approximate atoms are considered and the idea is rather to simultaneously test multiple atoms which eventually lie within a same zone S, reducing the total number of tests performed.

C. Building a safe region using an approximate dictionary

As in classical screening, we have all ingredients to screen provided we can build a safe region using moderate computational overhead. The object of this section is precisely to adapt the constructions of safe regions from Section II-C, which depend on A, to define new safe regions reusing computations done during the iterations with à instead of A.

Remark 2. Let us emphasize that a safe region is a region R containing the dual solution θ ⋆ of the original dual problem [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] and not necessarily the dual solution of its approximate version with ∆ Ã instead of ∆ A . Indeed, although A is approximated to speed up computations, variable elimination needs to be guaranteed with respect to the original problem (1).

The safe regions from Section II-C are built by determining a feasible dual point θ t ∈ ∆ A . Given an arbitrary z ∈ R N , the function Θ(z|A) from ( 10) would provide such a feasible point, however it cannot be used as its computation requires multiplication by A T . A naive alternative could be to compute Θ(z| Ã), however while this always belongs to ∆ Ã (the feasible set for Ã) it does not necessarily belong to the desired feasible set ∆ A . This can be fixed using a modified dual scaling approach that we propose to call stable dual scaling.

Considering error bounds ǫ = (ǫ j ) K j=1 ∈ R K + we define

Θ ′ (z| Ã, ǫ) := z max 1, max 1≤j≤K (|ã T j z| + ǫ j z q * ) (21)
For ǫ = 0 we recover Θ(z|A) = Θ ′ (z|A, 0).

Lemma 2. Assume that a j ∈ S j = B q (ã j , ǫ j ) for 1 ≤ j ≤ K. Then, for any z ∈ R N , Θ ′ (z| Ã, ǫ) ∈ ∆ A ∩ ∆ Ã is
θ ′ t = y T ρt λ ρt 2 2 α ′ t -α ′ t ρt , with α ′ t = max 1≤j≤K |ã T j ρt | + ǫ j ρt q * -1 (22)
These principles yield new stable safe ℓ 2 -spherical regions.

1) Stable Static Safe sphere:

The static safe sphere (12) depends on λ max = max j |a T j y|. To perform only computations with Ã, we can reuse [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF] with

λ ′ max := max j |ã T j y| + ǫ j y q * ≥ λ max , (23)
leading to R = B(c ′ , R ′ ) with Stable Static Safe: c ′ = y/λ, R ′ = | 1 λ ′ max -1 λ | • y 2 . ( 24
)
2) Stable Dynamic Safe sphere: We adapt the Dynamic Safe sphere in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] to our approximate setting by replacing the dual feasible point θ t [START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF]

by θ ′ t (22), leading to R = B(c ′ , R ′ ) with Stable Dynamic Safe: c ′ = y/λ, R ′ = θ ′ t -y/λ 2 . ( 25
)
The fact that θ ′ t ∈ ∆ A directly implies that this sphere is safe by definition of the basic ℓ 2 -spherical bound [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF].

3) Stable GAP Safe sphere:

To build a safe sphere B(c ′ , R ′ ) with c ′ = θ ′ t (22) instead of c = θ t (11) we need to determine a safe radius R ′ depending only on Ã. For this reason, instead of using the (inaccessible) du-

ality gap G(x t , θ ′ t |A) = P (x t |A) -D(θ ′ t ) we will use G(x t , θ ′ t | Ã) = P (x t | Ã) -D(θ ′ t )
and add a security margin ensuring that B(θ ′ t , R ′ ) is safe. Recall the standard notation M p→q := sup u p ≤1 Mu q .

Theorem 1 (Stable GAP Safe sphere). Consider

E ≥ A -Ã r→2 where 1 ≤ r ≤ ∞. For any 1 (x, θ) ∈ R K ×∆ A , we have θ ⋆ ∈ B c ′ = θ, R ′ = 1 λ 2G(x, θ| Ã) + 2δ(x) (26) with δ(x) := y -Ãx 2 E x r + 1 2 E 2 x 2 r ( 27
)
In other words B(θ, R ′ ) is safe.

The proof was given in [START_REF]Faster and still safe: combining screening techniques and structured dictionaries to accelerate the lasso[END_REF]. By applying Theorem (1) to x t and θ ′ t in [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF] we obtain 27) with x = x t . Test expressions for µ(a|B(c, R)) -see ( 8) -and µ(S ′ |B(c ′ , R ′ )) -see [START_REF] Xiang | Learning sparse representations of high dimensional data on large scale dictionaries[END_REF] where S ′ is defined in [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF] with a = a 2 .

R = B(c ′ , R ′ ) with Stable GAP Safe: c ′ = θ ′ t , R ′ = 1 λ 2G(x t , θ ′ t | Ã) + 2δ(x t ). ( 28
)
Sphere test Center Radius Test expression -µ(a|B(c, R)) or µ(S ′ |B(c ′ , R ′ )) Static Safe y/λ |1/λmax -1/λ| • y 2 |a T y|/λ + |1/λmax -1/λ| • y 2 • a 2 Stable Static Safe y/λ |1/λ ′ max -1/λ| • y 2 |ã T y|/λ + ǫ y q * /λ + |1/λ ′ max -1/λ| • y 2 • a 2 Dynamic Safe y/λ θt -y/λ 2 |a T y|/λ + θt -y/λ 2 • a 2 Stable Dynamic Safe y/λ θ ′ t -y/λ 2 |a T y|/λ + θ ′ t -y/λ 2 • a 2 , if a T y is given 2 |ã T y|/λ + ǫ y q * /λ + θ ′ t -y/λ 2 • a 2 , otherwise GAP Safe θt 1 λ 2G(xt, θt|A) |a T θt| + 1 λ 2G(xt, θt|A) • a 2 Stable GAP Safe θ ′ t 1 λ 2G(xt, θ ′ t | Ã) + 2δ(xt) |ã T θ ′ t | + ǫ θ ′ t q * + 1 λ 2G(xt, θ ′ t | Ã) + 2δ(xt) • a 2
The first columns of Table I summarize the obtained safe regions alongside with their previously existing analogue.

Combining stable safe sphere regions with stable safe screening tests (cf Section III-B) yields screening rules which remain safe despite manipulating an approximate dictionary (last column of Table I). As before, these new rules also reuse the calculations performed by optimization algorithm with the approximate dictionary Ã, making them suitable in practice.

IV. A FAST ALGORITHM FOR ℓ 1 REGULARIZATION

We now have the tools to combine stable safe screening and the use of fast structured dictionaries. The proposed Algorithm 1 consists in incorporating a stable dynamic screening rule to iterations of a proximal algorithm. Besides the vector y and the regularization parameter λ, the input of the algorithm consists in a sequence { Ãi } I i=0 of approximate dictionaries with ÃI = A, and a corresponding sequence of error bounds {ǫ i } I i=0 , where ǫ i = (ǫ i j ) K j=1 ∈ R K + contains ℓ 2 error bounds on the atoms of the i-th dictionary Ãi , and of course ǫ I = 0.

In general terms, the proposed strategy consists in gradually switching to more accurate dictionary approximations as the optimization algorithm approaches the solution. Note that any iterative optimization technique can be used by replacing the generic update function p(x t , A, α t ) by the corresponding update step, for instance equation (3) for ISTA. A possible variation consists in performing the screening at regular intervals instead of every iteration.

To fully describe this algorithm we need to specify the stopping and switching criteria, which involve a parameter Γ. In the following subsections, we present a switching criterion which assumes the approximations { Ã} I i=0 to be arranged in increasing order of precision.

A classical stopping criterion is to stop when the duality gap falls below a threshold, recalling that strong duality holds for problem (1), i.e. G(x ⋆ , θ ⋆ |A) = 0 [START_REF] Osborne | On the lasso and its dual[END_REF]. The simplest way of computing the duality gap is to wait until the algorithm 2 Since the center c = y/λ is independent of the iteration count t, the resulting test can be simplified provided that A T y can be pre-computed.

Algorithm 1 x ⋆ = FastL1({ Ãi } I i=0 , {ǫ i } I i=0
, y, λ, Γ) 1: Initialize: t = 0, x 0 = 0, i = 0, A 0 = {1, . . . , K}, 2: α 0 according to the considered proximal algorithm 3: while stopping criterion not met do --Restrict to preserved set --

5: Ã ← Ãi [At] , x t ← (x t ) [At]
6:

--Update preserved coefficients --

7: {x t+1 , α t+1 } ← p(x t , Ã, α t ) 8:
--Dynamic Screening --

Set θ ′ t using ( 22) and c ′ , R ′ according to Table I 10:

A t+1 ← {j ∈ A t : µ(S ′ j |B(c ′ , R ′ )) ≥ 1} 11: t ← t + 1 12:
--Switching criterion --

13:
Compute γ t (cf. ( 29)) and K t (cf. (39))

14:

i = SwitchDictionary(i, I, γ t , Γ, K t )

reaches à = ÃI = A. Alternatively one could upper bound the duality gap using approximate dictionaries, this is left to future work.

A switching criterion has two main motivations:

• Convergence: to avoid diverging from the solution of the original problem (1). The higher the error bounds ǫ i , the higher the risk of moving away from the true solution, therefore the need to switch earlier in such cases. • Speed: once screened enough, the original dictionary becomes faster to apply than the approximate ones. This is crucial at high regularization regimes (λ/λ max ≈ 1) in which screening typically occurs very quickly.

A. Switching: convergence criterion

In [START_REF] Dantas | Dynamic screening with approximate dictionaries[END_REF], [START_REF]Faster and still safe: combining screening techniques and structured dictionaries to accelerate the lasso[END_REF], we proposed a empirically-tuned heuristic convergence criterion. In the following, we derive a more theoretically grounded criterion.

Given a solution estimate x t at iteration t, a natural convergence measure for the problem L(λ, A, y) is the duality gap G(x t , θ t |A) (14) with θ t calculated as in [START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF]. However, when manipulating approximate dictionaries, this quantity is not accessible without unwanted additional computations.

G(xt, θ t ′ Ã), j =10 -1 G(xt, θ t ′ Ã), j =10 -2 G(xt, θ t ′ Ã), j =10 -3 G(xt, θt′A)
As a surrogate, we propose to use a ratio

γ t = G(x t , θt | Ã) G(x t , θ ′ t | Ã) (29)
between two computable gaps arising from two different dual points θ ′ t and θt that are respectively feasible for problems L(λ, A, y) and L(λ, Ã, y):

θ ′ t ∈ ∆ A ∩ ∆ Ã -→ G(x t , θ ′ t | Ã) = P (x t | Ã) -D(θ ′ t ) (30) θt ∈ ∆ Ã -→ G(x t , θt | Ã) = P (x t | Ã) -D( θt ) (31)
where θ ′ t ∈ ∆ A ∩ ∆ Ã has been previously calculated for the stable screening with [START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF], while θt ∈ ∆ Ã is the conventional dual point [START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF] but calculated with Ã. Then again, apart from minor extra memory requirements, the computation of the quantities in [START_REF] Tibshirani | Strong rules for discarding predictors in lasso-type problems[END_REF] reuses most of the calculations3 in [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF].

If à = Ãi was kept until convergence, we would have

xt → x⋆ and θt → θ⋆ . ( 32
)
with x⋆ a solution of L(λ, Ã, y) and θ⋆ a solution of the corresponding dual problem. Then, as illustrated in Fig. 1 (left), the second quantity G(x t , θt | Ã) -which is precisely the duality gap w.r.t. the problem L(λ, Ã, y)would tend to zero, while the first quantity G(x t , θ ′ t | Ã) would typically saturate at a nonzero value since θ ′ t θ⋆ . As a result, the ratio γ t would shrink with the iterations as illustrated in Fig. 1 (right).

With that in mind, we propose to switch dictionaries when

γ t ≤ Γ (33)
where the higher Γ, the greater sensibility to the gap saturation.

B. Speed criterion

In terms of speed, the general idea is to exploit tradeoffs between approximation and computational efficiency.

Denoting C(A, A) the computational cost of multiplying matrix A [A] by a vector, it is reasonable to assume that C(A, A) ≤ C(A, A ′ ) if A ⊂ A ′ . As a result, given the nested structure of the sequence of preserved sets A t , the sequence C(A, A t ), t ≥ 0 is non-increasing.

Consider two approximate dictionaries Ãi and Ãj and assume the initial "unscreened" computational costs satisfy C( Ãi , A 0 ) < C( Ãj , A 0 ). Denote E i , E j the corresponding approximation errors (cf (4)), and assume that E i is "larger" than E j (let us denote this E i ≻ E j ) in the sense that for each column the approximation error is larger. An irrevocable switching point would be an iteration T in which the coarser approximation Ãi becomes more complex than Ãj

∀t > T, C( Ãi , A t ) ≥ C( Ãj , A t ). (34)
Note that this only happens if C( Ãi , A t ) decreases faster than C( Ãj , A t ) as screening progresses (i.e. as |A t | gets smaller).

To use the general idea just described, one needs to specify a complexity model C(A, A t ) depending on the structure of the corresponding matrix. For a generic unstructured dictionary matrix A ∈ R N ×K , a simple model is:

C(A, A t ) = N |A t | (35)
To quantify the complexity reduction on matrix-vector products entailed by an approximate dictionary we will use the concept of Relative Complexity (RC) [START_REF] Magoarou | Flexible multilayer sparse approximations of matrices and applications[END_REF], such that

C( Ã, A 0 ) = RC( Ã)×N K. (36)
In a worst-case scenario, screening does not further reduce the cost and C( Ã, A t ) = RC( Ã) × N K. In more optimistic scenarios, some fast approximate dictionary structures might still benefit from speedups upon column removal. Consider approximate dictionaries { Ã0 , . . . , ÃI }, with ÃI = A, of increasing precision and decreasing complexity E 1 ≻ . . . E I = 0 and RC( Ã0 ) < . . . < RC( ÃI ) = 1. [START_REF] Johnson | Blitz: A principled meta-algorithm for scaling sparse optimization[END_REF] The above observations suggest to switch from Ãi to

A if RC( Ãi )×N K = C( Ãi , A t ) ≥ C(A, A t ) = N |A t |, giving |A t | ≤ RC( Ãi ) K. (38)
Also, as a consequence of complexity model (36) adopted for the approximate dictionaries Ã, this criterion never causes a switching from Ãi to Ãj with i < j = I.

C. Heuristic look-ahead speed criterion

The preserved set size |A t | used in [START_REF] Massias | From safe screening rules to working sets for faster lasso-type solvers[END_REF] should ideally be the one associated to the dictionary after switching (A in this case), while in practice only the screening rate associated to the current approximation Ãi is available. This motivates the introduction of a heuristic (that preserves the safety of screening) to anticipate the potentially smaller size of the preserved set |A t | as soon as we switch to a more precise approximation (or to the original dictionary).

Indeed, stable screening generally leads to less atom eliminations due to the extra security margins (and more so for higher approximation errors). This is illustrated in Figure 2 (left) with a scatterplot in which the x-axis correspond to the "oracle" number of preserved atoms |A| obtained with a conventional Dynamic Safe screening mechanism having access to A. This is the quantity we are interested in estimating for the speed criterion [START_REF] Massias | From safe screening rules to working sets for faster lasso-type solvers[END_REF]. The y-axis corresponds to the number of preserved atoms obtained by the corresponding stable test for a given approximation Ã. Each plotted point compares the two mentioned screening ratios at a given iteration. The procedure was repeated a large number of times in order to give a full picture of how the two quantities correlate to one another. The main observation is that all points lie above the identity line which means that the stable test always overestimates (sometimes significantly) the value that |A t | will reach after switching. This phenomenon is intensified for higher approximation errors.

As a heuristic, we propose to estimate |A t | using the conventional test on the approximate atoms, i.e. µ(ã j |R).

K t = |{j ∈ {1, . . . , K} : µ(ã j |R) ≥ 1}| (39)
As shown in Fig. 2 (right graph), this correlates much better with the oracle value of |A t |, since it doesn't have the additional security margins used in the stable test. Even though the test used to compute K t is unsafe with respect to the original problem [START_REF] Dantas | Stable screening -Python code[END_REF], it has no impact on the safety of screening itself, which is still performed with the stable test. Moreover, the additional test is virtually costless since the same calculations of the first test are used.

D. Summary and example

The resulting switching strategy, shown in Algorithm 2, plays two roles:

(I) Decide when to stop using the current dictionary. (II) Choose the next dictionary. Algorithm 2 j = SwitchDictionary(i, I, γ t , Γ, K t )

1: if K t ≤ RC( Ãi ) K then ⊲ Speed criterion 2: j ← I 3: else if γ t ≤ Γ then ⊲ Convergence criterion 4: j ← i + 1 5: else 6: j ← i
The speed criterion (section IV-B) triggers a switching directly to the original dictionary ÃI = A. We suppose the Relative Complexity associated to the i-th approximation RC( Ãi ) (cf. [START_REF] Kim | Fast active-set-type algorithms for l1-regularized linear regression[END_REF]) to be known and stored in memory. The 29)) for the convergence switching criterion, with threshold Γ = 0.2. (d) Computational complexity per iteration (worst-case theoretical number of flops, cf. section IV-B). N = 2500, K = 10000, λ/λ max = 0.2. convergence criterion (section IV-A) switches to the next available approximation Ãi+1 . Note that if A is already adopted, i.e. input i = I, then A is kept.

Example: Figure 3 shows an example of the proposed algorithm, given a list of approximations { Ãi } I i=0 with I = 4, decreasing approximation error ǫ i and increasing relative complexity RC( Ãi ) = 0.15(i+1), 0 ≤ i < I. As usual, ÃI = A.

Figure 3a shows the duality gap evolution over the iterations with Γ = 0.2. As soon as the gap G(x t , θ ′ t | Ãi ) saturates, a more precise approximation is adopted (the next on the list), avoiding to excessively delay the convergence compared to the conventional case where screening is performed with A from the beginning (dotted curve). Although some delay is introduced in terms of iterations until convergence, it is largely compensated by the fact that the initial iterations with the approximate dictionaries are much faster. This is illustrated in Fig. 3b in which the duality gap is plotted as a function of the execution time.

The theoretical complexity per iteration is shown in Fig. 3d. As better approximations are adopted, the iteration complexity rises proportionally to the corresponding RC. Note that as soon as the conventional dictionary becomes faster than the current approximation, it is promptly adopted (by the speed criterion). In this example, the last approximation Ãi=3 was never used, since we switched directly from Ãi=2 to A.

E. Complexity Analysis

Existing screening tests introduce only a minor computational overhead because they primarily reuse matrix-vector multiplications either already performed in the optimization algorithm update p(x t , A, α t ) (typically the product A T ρ) or that can be precalculated once for all (A T y). The same holds for the stable tests proposed in this article. We now derive the involved number of floating point operations (flops). Since static screening rules represent only a fixed computational overhead, we concentrate on dynamic rules which could potentially (if not properly designed) lead to a significant overhead.

1) Screening cost: The most expensive computations associated to a dynamic screening rule -the ones potentially of the order of a matrix vector product O(N K) -are:

• Computation of a dual feasible point θ t in (11) (resp. θ ′ t ( 22)) requires the product A T ρ (resp. ÃT ρ) reused from the optimization algorithm.

• Sphere test µ(a|R) (resp. µ(S ′ j |R)) requires the product a T j c (resp. ãT j c) for all preserved atoms j ∈ A, which comes down to the matrix-vector product A T c (resp. ÃT c). Practical sphere regions have been designed to limit this potential overhead. While in the Dynamic Safe sphere it reduces to the precalculated product A T y, in the GAP Safe sphere it reduces to A T ρ (resp. ÃT ρ) calculated in the optimization iteration. The other required calculations are detailed in Appendix D. In short, the screening represents a rather low overhead O(N + |A|) -even its stable version-compared to the optimization update: O(N |A|) with screening or O(N K) without it, due to matrix-vector products.

2) Full iteration cost: Table IV shows the number of operations of a complete iteration in Algorithm 1 (ISTA update + screening), following [START_REF] Bonnefoy | Dynamic screening: Accelerating first-order algorithms for the lasso and grouplasso[END_REF] and adopting the complexity models in equations ( 35) and ( 36) for matrix vector multiplications. We denote flops A (t) the cost of iteration t with the conventional screening and flops Ã(t) with the stable screening. As a benchmark, we use the complexity of an ISTA iteration without screening, denoted flops N (t).

flops N (t) (K + xt 0 )N + 4K + N flops A (t) (|A| + xt 0 )N + 6|A| + 5N flops Ã(t) (RC×K + xt 0 )N + 8|A| + 7N

TABLE II: Complete iteration complexity

To obtain the total complexity of the algorithm, simply add up all iteration costs calculated with the corresponding active set size |A t | and sparsity of the solution estimate x t 0 .

V. EXPERIMENTS

In this section we demonstrate the potential of the proposed algorithm in terms of complexity reduction and time saving for ℓ 1 -minimization problems. This is done in a wide set of simulation scenarios, summarized in Table III, to evaluate the influence of the main parameters involved.

The dimension of the observed vector y is N = 2500 and that of the coefficient vector x (or, equivalently, the number of atoms in the dictionary) is K = 10000. Unitnorm input data samples y = Ax are generated from a sparse vector x with active set determined by a Bernoulli distribution with probability p = 0.02 and zero-mean standard Gaussian entries. The reported results are the average on 25 independent and identically distributed realizations of the input y. The dictionary matrix A ∈ R N ×K is generated in such a way that it is more or less efficiently approximated by a fast structured matrix according to three representative scenarios (more details are given in Section V-A).

Although other problem dimensions (N and K) and data distributions were explored, we decided to keep these parameters fixed on the reported experiments since they were observed not to decisively influence the analyzed results.

In all figures, the GAP Safe and Dynamic Safe tests are denoted respectively by the acronyms GAP and DST.

A. Specifying the Fast Structured dictionaries

The performance of the proposed algorithm is directly associated to the quality of the available approximations Ãi . A good approximation would have both a small application complexity (RC) and small approximation error (ǫ). There is a compromise, however, since each of these features usually come to the price of the other.

In the experiments, we use a particular kind of fast structured dictionaries referred to as SuKro [START_REF] Dantas | Learning dictionaries as a sum of Kronecker products[END_REF] which can be written as a sum of Kronecker products à = nkron k=1 B k ⊗C k . Its reduced multiplication cost comes from the fact that the sub-matrices

B k , C k ∈ R √ N× √
K are much smaller than A ∈ R N×K . The choice of the SuKro structure is justified by the fact that it directly provides a range of speedup-error compromises by tweaking the number of Kronecker terms (n kron ) in the sum. A higher n kron provides a more precise approximation although implying a higher RC.

We define three representative simulation scenarios -easy, moderate and hard scenarios -in which the dictionary matrix is poorly, moderately and efficiently approximated by the structured dictionary respectively. In practical terms, this translates to different approximation error decay profiles as a function of the number of Kronecker terms (n kron ) on the SuKro structure, as shown in Figure 4. The easier the scenario, the faster the approximation error decays as a function of n kron .

Although computational complexity associated to a certain SuKro operator can be calculated analytically, we measured the actual time speedup obtained in practice in order to have a more realistic value. Interestingly, the measured RC is lower than the theoretical prediction.

B. Computational complexity and Time results

Let us denote F N , F A and F Ã respectively the total computational complexity (in number of flops) of the optimization algorithm without screening, with the conventional screening and with the proposed approach, such that

F X = nit t=1 flops X (t), X ∈ {N, A, Ã} (40)
where n it is the number of iterations. When calculating F Ã , the current approximation at iteration t must be used in the expression of flops Ã(t). Likewise, we denote T N , T A and T Ã the measured running times. As a main figure of merit, we evaluate the normalized number of flops (F A /F N and F Ã/F N ) and normalized running times (T A /T N and T Ã/T N ) 5 .

In all simulated scenarios, the observed time reductions match closely the theoretical speedups predicted in terms of computational complexity, as illustrated in Figure 5 for the GAP Safe rule. It shows the correlation between theoretical complexities and measured running times in multiple independent runs. The fact that the points are well concentrated around the identity line is an empirical evidence that the predicted speedups really translate into practical accelerations. Similar results are obtained for the Dynamic Safe rule, for other convergence tolerances and switching thresholds. Given this observation, in the remainder of the paper we report only running time results, which are more relevant in practice.

C. Choosing the switching threshold

The convergence-based switching criterion proposed in Section IV-A relies on a hyperparameter: the threshold Γ. As discussed then, this parameter determines how long the approximate dictionaries are kept.

We empirically observed that the choice of Γ is mostly dependent on the quality of the approximations Ãi (represented here by the three simulation scenarios defined in Section V-A).

In Figure 6 we show the normalized times as a function of Γ ∈ [10 -2 , 0.8] for each one of the three scenarios. Each line corresponds to a different regularization level λ. For each of these lines, we are interested in the Γ value that minimizes the running times. We can see that the low regularization configurations (dark blue curves) are more sensitive to choice of Γ, especially in the Hard scenario (left plot) in which a bad choice of Γ can even lead to normalized running times greater than one -so, actually, a slowdown. This happens when a too small Γ is used, which means that the approximate dictionaries are kept longer than they should, causing an important detour in the convergence path and thus delaying convergence.

There is no reason for the optimal Γ to be the same for every regularization level. However, a common behavior is observed regardless of the regularization and, for a given simulation scenario, a single Γ value can be chosen to obtain close to optimal execution times for any λ/λ max in the tested range. Slightly different Γ values can be chosen to better adapt to each simulation scenario: Γ = 0.5, 0.25, 0.2 (indicated by the dotted vertical lines in the figure) respectively on the hard, moderate and easy scenarios. In general terms, the worse the available approximations are (i.e. the harder the scenario) the higher the Γ to be chosen. This is intuitive, since it implies being more conservative in switching earlier when the available approximations are of lower quality.

If the aimed speedup-error compromise is not fixed, a good general compromise is picking Γ = 0.5, for instance.

D. Single vs. Multiple Approximations

In our experiments, using multiple dictionary approximations proved to be always advantageous when compared to using a single approximation, as soon as the switching parameter Γ is well-tuned (see Section V-C).

Figure 7 shows the normalized running times for an entire range of regularizations and corroborates the previous statement for both GAP and Dynamic Safe tests as well as for both ISTA and FISTA algorithms. The medians among 25 runs are plotted and the shaded area contains the 25%-to-75% percentiles. Although the single-approximation version of the proposed algorithm (as introduced in [START_REF] Dantas | Dynamic screening with approximate dictionaries[END_REF], [START_REF]Faster and still safe: combining screening techniques and structured dictionaries to accelerate the lasso[END_REF]) already provides noticeable speedups to the tested gradient-based solvers, even with respect to their screening-based implementation (black curves), the use of multiple approximations (red curves) consistently leads to even better speedups. This ratifies the relevance of the generalization to multiple dictionaries introduced in this paper, by relying on more robust switching criteria. Similar results, which are omitted here, were obtained for other convergence tolerances and speedup-error compromise (these parameters are further explored in Sections V-E and V-F).

A smaller speedup is provided to the FISTA algorithm when compared to ISTA, which is expected since the former is already faster than the latter 6 . A similar argument applies to the GAP Safe rule, which has stronger screening capabilities than Dynamic Safe: our method still manages to provide some additional acceleration when combined to these already quite efficient techniques, especially in low-regularized scenarios.

For simplicity, only the results concerning multiple dictionary approximations will be reported from this point on.

E. Speedup-error compromise

It is reasonable to expect the quality of the approximations Ãi to be decisive on the performance of the proposed algorithm. Its success depends on the possibility of providing fast yet precise approximations of the dictionary A. Let us now evaluate the influence of this speedup-error compromise on the simulation time results.

Figure 8 shows the same type of normalized running times as in Figure 7, but it now compares the results of the proposed algorithm under three different scenarios in terms of speeduperror profiles (see Figure 4 for details on the approximation errors in each scenario). The impact on the algorithm performance is indeed significant. However, even on the Hard scenario, the proposed approach manages to provides nonnegligible acceleration.

F. Convergence tolerance

As previously highlighted, screening is known to work better on highly-regularized scenarios (λ/λ max ≈ 1). The speedup provided by the screening is also more pronounced when a more precise convergence is required, because the final iterations are often less costly as most atoms have been screened out. These two aspects create the triangular profiles observed in Figure 9 (right plots). This figure shows the normalized execution times in grayscale (the darker, the slower) as a function of both the regularization (x-axis) and the convergence precision in terms of the duality gap (y-axis).

Note that the proposed method (left plots) efficiently complements the screening tests in its main weaknesses, namely: low regularizations and reduced convergence requirements. In weakly regularized scenarios, while screening tests struggle to eliminate atoms especially in initial iterations, the fast approximate dictionaries come at rescue by making those iterations faster. Besides that, when higher duality gaps are targeted, less iterations using the slow dictionary A are necessary. Thus, proportionally, the accelerated part of the algorithm (with fast dictionaries) is more significant.

G. MEG source localization

In this experiment, we consider an MEG (magnetoencephalography) source localization problem, which consists in retrieving a limited set of active brain foci associated to a given MEG measured signal. A common way to achieve such a task is to solve a convex sparse inverse problem [START_REF] Matsuura | Selective minimum-norm solution of the biomagnetic inverse problem[END_REF], [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF].

We consider an MEG gain matrix (forward operator) A ∈ R 204×8193 obtained with the MNE software [START_REF] Gramfort | MNE software for processing MEG and EEG data[END_REF]. This operator can be efficiently approximated as a product of a few sparse matrices, as proposed in [START_REF] Magoarou | Flexible multilayer sparse approximations of matrices and applications[END_REF]. The speedups and corresponding approximation errors achieved by the so-called FAµST dictionaries adopted here are summarized in Table V. The obtained speedup-error compromise corresponds to a hard scenario.

The experiment consists in picking (uniformly) at random eight active sources with standard gaussian weights giving 8-sparse coefficient vectors x ∈ R 8193 to be recovered from the input signals y = Ax ∈ R 204 by solving problem (1). Figure 10 shows the required time to reach a certain duality gap precision (averaged over 50 independent runs) for a fixed regularization λ/λ max = 10 -1 , around which the correct number of sources is recovered. The proposed method is consistently advantageous, providing some speedup with respect to both the FISTA algorithm (grey dotted line) and its enhanced version with GAP Safe screening (black line). Although somewhat modest when compared to the latter, the provided gain is proportionally more pronounced at lower precision requirements in the duality gap. 7 More significant speedups would require better approximations of MEG gain matrix (in the sense of speedup-error compromise).

VI. CONCLUSION

Aiming to combine safe screening tests and fast structured operators in an effort to accelerate the solution of ℓ 1 -minimization problems, we proposed a methodology for defining safe screening tests despite an inaccurate knowledge of the explanatory variables (dictionary atoms). The proposed methodology was exemplified on some existing screening tests (stable and dynamic). The resulting stable screening tests were then employed in a fast algorithm, which exploits a series of fast approximations of the dictionary matrix. Simulation results demonstrated the two combined strategies to be quite complementary, justifying the effort to conciliate them.

In terms of longer term perspectives, the first part of the paper also provides a better understanding on which error measures are more significant to the stable screening tests, potentially leading to more suited optimization criteria in structured dictionary approximation techniques. In a broader sense, the proposed stable tests can be seen as a robust screening tool for inaccurately-defined ℓ 1 -minimization problems in which the imprecision need not be due to the use of fast approximations -for instance, due to intrinsic imprecisions on the measurement matrix of a given inverse problem. In principle, the proposed framework can also be applied to other types of safe regions and error bounds than the ones explored in this paper, such as dome tests [START_REF] Xiang | Screening tests for lasso problems[END_REF]. Finally, an interesting perspective is to extend these zone-based screening tests to off-the-grid generalizations of ℓ 1 -regularized problems [START_REF] Bredies | Inverse problems in spaces of measures[END_REF]. q ) + , which results from the Holder's inequality. Substituting [START_REF] Gramfort | MNE software for processing MEG and EEG data[END_REF] and [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF] in [START_REF] Matsuura | Selective minimum-norm solution of the biomagnetic inverse problem[END_REF] gives the bound in equation [START_REF] Dantas | Learning dictionaries as a sum of Kronecker products[END_REF].

C. Proof of Lemma 2

By definition we have θ F := Θ ′ (z| Ã, ǫ) = αz where 0 ≤ α ≤ 1 max j |ã T j z| + ǫ j z q * As a result for any 1 ≤ j ≤ K we have |a T j θ F | = |(ã j + e j ) T αz| ≤ |α| |ã T j z| + ǫ j z q * ≤ 1. This implies that the dual point θ F is feasible, i.e. θ F ∈ ∆ A . The fact that θ F ∈ ∆ Ã follows similarly (since ǫ j z q * > 0).

D. Detailed complexities

We suppose that the p * -norms of the atoms a j p * and qnorms of the atom approximation error e j q are precalculated and stored and memory. The total screening complexity is summarized in Table IV 22): N norm ρt q * (only if q * = 2). |A t | products ǫ j ρt q * on the calculation of α ′ t .

2) Safe region R: At most N operations for computing θ ty/λ 2 (if the dual objective D(θ t ) is not already calculated as a convergence criterion) for the radius of either Dynamic Safe or GAP Safe sphere.

3) Screening test: Other operations (aside from the mentioned A T c) sum up to |A t | operations for the product R a j p * ∀ j ∈ A t . The proposed stable screening tests require at most an extra N + |A t | operations: N norm c q * (if q * = 2 and c varies with t). |A t | products ǫ j c q * ∀ j ∈ A t (if c varies with t).

  [z] b a := min(max(z, a), b) denotes the projection of the scalar z onto the segment [a, b].

3 Fig. 1 :

 31 Fig. 1: Typical behaviour of duality gaps (left) and ratio γ t (right) as a function of the iteration number and error bounds.

Fig. 2 :

 2 Fig.2: Number of preserved atoms |A| using oracle conventional screening (x-axis) vs stable screening (y-axis, left plot), or vs proposed heuristic (y-axis, right plot). .

Fig. 3 :

 3 Fig. 3: Application example of the proposed algorithm. (a) and (b) Duality gap over the iterations and time respectively. (c) Gap ratio γ t (cf. (29)) for the convergence switching criterion, with threshold Γ = 0.2. (d) Computational complexity per iteration (worst-case theoretical number of flops, cf. section IV-B). N = 2500, K = 10000, λ/λ max = 0.2.

4

 4

Fig. 4 :

 4 Fig. 4: Simulated scenarios: different trade-offs between approximation error and speedup of the structured dictionaries.

  n

Fig. 5 :

 5 Fig. 5: Normalized number of flops vs. normalized running times for ISTA (tol = 10 -5 , Γ = 0.5). Each point on the graph corresponds to an independent run at a given regularization λ/λ max ∈ [10 -2 , 1]. Results for the usual GAP Safe test are shown in black and those using its stable version (as in Alg. 1) are shown in blue. The closer the dots lie to the identity line, the better the theoretical speedups are met in practice.

Fig. 6 :

 6 Fig. 6: Impact of the switching threshold Γ (x-axis) on the time results (y-axis) of Algorithm 1 using FISTA and GAP Safe. Each line corresponds to a regularization λ/λ max ∈ [10 -2 , 1].

  X = Ã, nkron = [5] X = Ã, nkron = [10] X = Ã, nkron = [15] X = Ã, nkron = [20] X = Ã, nkron =[START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Chabiron | Toward fast transform learning[END_REF][START_REF] Xiang | Learning sparse representations of high dimensional data on large scale dictionaries[END_REF]

Fig. 7 :Fig. 8 :

 78 Fig. 7: Running times normalized w.r.t. T N (ISTA or FISTA algorithms without screening). Conventional and Stable screening results respectively in black and colored lines. Left: ISTA. Right: FISTA. Top: Dynamic Safe screening. Bottom: GAP Safe. Moderate scenario, tol = 10 -5 , Γ = 0.5.

Fig. 9 :

 9 Fig. 9: Normalized execution times in grayscale (darker means slower) to reach a certain duality gap for a range of regularizations. Left: proposed technique (w/ stable screening). Right: Conventional screening.

Fig. 10 :

 10 Fig. 10: Source localization using an MEG forward operator. Duality gap as a function of time for λ/λ max = 10 -1 .

APPENDIXA. 1 p * - 1

 11 Proof of Equation[START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] Screening test for a ℓ p -ball, i.e. R = B p (c, R):sup θ∈Bp(c,R) |a T θ| = sup u∈Bp(0,1) |a T (c + Ru)| = |a T c| + R sup u∈Bp(0,1) |a T u| = |a T c| + R a p *(41)using the definition of dual norm a p * := sup u∈Bp(0,1) |a T u|.B. Proof of Equation (18)From Definition 3 and the conventional sphere test in eq. (8), the stable sphere test is given byµ(S|B p (c, R)) = sup a∈S |a T c| + R a p * ≤ sup a∈S |a T c| + R sup a∈S a p *(42)Taking S = B q (ã, ǫ), the first term gives (similarly to[START_REF] Osborne | On the lasso and its dual[END_REF])sup a∈S |a T c| = sup u∈Bq(0,1) |(ã + ǫu) T c| = |ã T c| + ǫ c q * (43)and the last term gives sup a∈Bq(ã,ǫ)a p * = sup u∈Bq(0,1) (ã + ǫu) p * = ã p * + Cǫ(44)with a constant C = sup u∈Bq(0,1) u p * = N (

  The proof is in Appendix C. Analogously to[START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF], we can now define a dual feasible point proportional to the residual ρt := y-Ãx t at iteration t, θ ′ t = Θ ′ proj ρt

	is to say more explicitly	y λ	| Ã, ǫ , that

a feasible point w.r.t. both the original dual problem

[START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]

and its modified version with à instead of A.

TABLE I :

 I Sphere center c and radius R for Static Safe, Dynamic Safe and GAP Safe (c ′ and R ′ for the stable versions). See definitions of λ max in (2), of λ ′ max in (23), of θ t in[START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF], of θ ′ t in[START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF], and of δ(x t ) in (

  θt|A) < tol ∈ [10 -6 , 10 -3 ]

	Problem parameters	Values
	Regularization λ/λmax	[10 -2 , 1] (logarithmically-spaced)
	Convergence tolerance G(xt, Data distribution y Ax with x Bernoulli-Gaussian
	Algorithm parameters	Values
	Optimization algorithm	ISTA, FISTA
	Screening Test	(stable) Dynamic, (stable) GAP Safe
	Switching threshold Γ	[0.01, 0.8] (logarithmically-spaced)
	Speedup-Error tradeoff Ãi	3 scenarios (see Section V-A)

TABLE III :

 III Pool of parameters explored in the simulations.

TABLE IV :

 IV Relative Complexities (RC)

TABLE V :

 V Approximation errors and Relative Complexities for the FAuST approximations of the MEG matrix.

  and detailed in the following.

	Screening	θt	R	Test	Total
	Dynamic Safe	3N + |A|	N	|A|	4N + 2|A|
	Stable Dynamic Safe 4N + 2|A| N	|A|	5N + 3|A|
	GAP Safe	3N + |A|	N	|A|	4N + 2|A|
	Stable GAP Safe	4N + 2|A| N	N + 2|A| 6N + 4|A|

TABLE VI :

 VI Screening: number of floating-point operations. 1) Feasible point θ θ θ t : A total of at most 3N + |A t | operations for the point in eq. (11), distributed as follows: N dot product y T ρ t . N norm ρ t 2 2 (if primal objective not calculated as a convergence criterion). N product αρ t . |A t | comparisons for the infinity norm A T ρ ∞ . At most N +|A t | extra operations for the point in eq. (

In this theorem, θ ∈ ∆ A can be arbitrary. In particular it is not necessarily expressed as in[START_REF] Fercoq | Mind the duality gap: safer rules for the lasso[END_REF], hence it does not need to belong to ∆ Ã.

Only N extra products for computing θ and |At| comparisons for ÃT ρ t ∞ are required.

Some parallelization is introduced (in the summing terms) and matrixmatrix products are faster than the naive cubic complexity.

Time results are obtained in an Intel R Core TM i7-5600U CPU @ 2.60GHz, 16GB RAM. But since mostly time ratios are reported, the results here should be relatively consistent to other machine specifications.

The results reported in the right graphs are already normalized w.r.t the FISTA time results, which are typically smaller than those of ISTA.

Neither of the two considered acceleration techniques affects the algorithms results in terms of source localization accuracy, since they solve the exact same problem at the same convergence precision. That is why taskspecific performance measures are not discussed (see[START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF] for more details on the performance of ℓ 1 -minimization techniques for brain source localization).

This work has been partially funded by the Becose ANR project (ANR-15-CE23-0021).

Code is available at: github.com/cassiofragadantas/Screening ADST [1].