Stable safe screening and structured dictionaries for faster L1 regularization

TL;DR

This paper introduces stable safe screening tests combined with structured dictionary approximations to accelerate L1-regularized least squares, achieving faster computation while maintaining solution accuracy.

Contribution

It proposes a new family of stable screening tests that handle approximation errors, enabling efficient dictionary approximations during optimization.

Findings

Significant reduction in computational complexity.

Faster execution times demonstrated across various scenarios.

Effective combination of safe screening with structured dictionary approximations.

Abstract

In this paper, we propose a way to combine two acceleration techniques for the $\ell\_{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.

Tables

Table 1. TABLE I : Sphere center 𝐜 𝐜 \bm{\mathrm{c}} and radius R 𝑅 R for Static Safe, Dynamic Safe and GAP Safe ( 𝐜 ′ superscript 𝐜 ′ \bm{\mathrm{c}}^{\prime} and R ′ superscript 𝑅 ′ R^{\prime} for the stable versions). See definitions of λ max subscript 𝜆 \lambda_{\max} in ( 2 ), of λ max ′ subscript superscript 𝜆 ′ \lambda^{\prime}_{\max} in ( 23 ), of 𝜽 t subscript 𝜽 𝑡 \bm{\mathrm{\theta}}_{t} in ( 11 ), of 𝜽 t ′ subscript superscript 𝜽 ′ 𝑡 \bm{\mathrm{\theta}}^{\prime}_{t} in ( III-C ), and of δ ​ ( 𝐱 t ) 𝛿 subscript 𝐱 𝑡 \delta(\bm{\mathrm{x}}_{t}) in ( 27 ) with 𝐱 = 𝐱 t 𝐱 subscript 𝐱 𝑡 \bm{\mathrm{x}}=\bm{\mathrm{x}}_{t} . Test expressions for μ ​ ( 𝐚 | B ​ ( 𝐜 , R ) ) 𝜇 conditional 𝐚 𝐵 𝐜 𝑅 \mu(\bm{\mathrm{a}}|B(\bm{\mathrm{c}},R)) – see ( 8 ) – and μ ​ ( 𝒮 ′ | B ​ ( 𝐜 ′ , R ′ ) ) 𝜇 conditional superscript 𝒮 ′ 𝐵 superscript 𝐜 ′ superscript 𝑅 ′ \mu(\mathcal{S}^{\prime}|B(\bm{\mathrm{c}}^{\prime},R^{\prime})) – see ( 20 ) where 𝒮 ′ superscript 𝒮 ′ \mathcal{S}^{\prime} is defined in ( 19 ) with a = ‖ 𝐚 ‖ 2 𝑎 subscript norm 𝐚 2 a=\|\bm{\mathrm{a}}\|_{2} .

Sphere test	Center	Radius	Test expression - $\mu \left(𝐚|B(𝐜,R)\right)$ or $\mu \left({𝒮}^{\prime }|B({𝐜}^{\prime },R^{\prime })\right)$
Static Safe	$𝐲/\lambda$	$\left|1/{\lambda }_{\max }-1/\lambda \right|\cdot {\Vert 𝐲\Vert }_2$	$\left|{𝐚}^T𝐲\right|/\lambda +\left|1/{\lambda }_{\max }-1/\lambda \right|\cdot {\Vert 𝐲\Vert }_2\cdot {\Vert 𝐚\Vert }_2$
Stable Static Safe	$𝐲/\lambda$	$\left|1/{\lambda }_{\max }^{\prime }-1/\lambda \right|\cdot {\Vert 𝐲\Vert }_2$	$\left|{\widetilde{𝐚}}^T𝐲\right|/\lambda +ϵ{\Vert 𝐲\Vert }_{q^{\ast }}/\lambda +\left|1/{\lambda }_{\max }^{\prime }-1/\lambda \right|\cdot {\Vert 𝐲\Vert }_2\cdot {\Vert 𝐚\Vert }_2$
Dynamic Safe	$𝐲/\lambda$	${\Vert {𝜽}_t-𝐲/\lambda \Vert }_2$	$\left|{𝐚}^T𝐲\right|/\lambda +{\Vert {𝜽}_t-𝐲/\lambda \Vert }_2\cdot {\Vert 𝐚\Vert }_2$
Stable Dynamic Safe	$𝐲/\lambda$	${\Vert {𝜽}_t^{\prime }-𝐲/\lambda \Vert }_2$	$\left|{𝐚}^T𝐲\right|/\lambda +{\Vert {𝜽}_t^{\prime }-𝐲/\lambda \Vert }_2\cdot {\Vert 𝐚\Vert }_2$, if ${𝐚}^T𝐲$ is given222Since the center $𝐜=𝐲/\lambda$ is independent of the iteration count $t$, the resulting test can be simplified provided that ${𝐀}^T𝐲$ can be pre-computed.
			$\left|{\widetilde{𝐚}}^T𝐲\right|/\lambda +ϵ{\Vert 𝐲\Vert }_{q^{\ast }}/\lambda +{\Vert {𝜽}_t^{\prime }-𝐲/\lambda \Vert }_2\cdot {\Vert 𝐚\Vert }_2$,   otherwise
GAP Safe	${𝜽}_t$	$\frac{1}{\lambda }\sqrt{2G({𝐱}_t,{𝜽}_t|𝐀)}$	$\left|{𝐚}^T{𝜽}_t\right|+\frac{1}{\lambda }\sqrt{2G({𝐱}_t,{𝜽}_t|𝐀)}\cdot {\Vert 𝐚\Vert }_2$
Stable GAP Safe	${𝜽}_t^{\prime }$	$\frac{1}{\lambda }\sqrt{2G({𝐱}_t,{𝜽}_t^{\prime }|\widetilde{𝐀})+2\delta \left({𝐱}_t\right)}$	$\left|{\widetilde{𝐚}}^T{𝜽}_t^{\prime }\right|+ϵ{\Vert {𝜽}_t^{\prime }\Vert }_{q^{\ast }}+\frac{1}{\lambda }\sqrt{2G({𝐱}_t,{𝜽}_t^{\prime }|\widetilde{𝐀})+2\delta \left({𝐱}_t\right)}\cdot {\Vert 𝐚\Vert }_2$

Table 2. TABLE II : Complete iteration complexity

${\text{flops}}_N\left(t\right)$	$\left(K+{\Vert {𝐱}_t\Vert }_0\right)N+4K+N$
${\text{flops}}_{𝐀}\left(t\right)$	$\left(\left|𝒜\right|+{\Vert {𝐱}_t\Vert }_0\right)N+6\left|𝒜\right|+5N$
${\text{flops}}_{\widetilde{𝐀}}\left(t\right)$	$\left(\text{RC}\times K+{\Vert {𝐱}_t\Vert }_0\right)N+8\left|𝒜\right|+7N$

Table 3. TABLE III : Pool of parameters explored in the simulations.

Problem parameters	Values
Regularization $\lambda /{\lambda }_{max}$	$[{10}^{-2},1]$ (logarithmically-spaced)
Convergence tolerance
Data distribution $𝐲$	$𝐀𝐱$ with $𝐱$ Bernoulli-Gaussian
Algorithm parameters	Values
Optimization algorithm	ISTA, FISTA
Screening Test	(stable) Dynamic, (stable) GAP Safe
Switching threshold $\Gamma$	$[0.01,0.8]$ (logarithmically-spaced)
Speedup-Error tradeoff ${\widetilde{𝐀}}^i$	3 scenarios (see Section V-A)

Table 4. TABLE IV : Relative Complexities (RC)

$n_{\text{kron}}$	5	10	15	20
Measured RC	0.08	0.15	0.21	0.28
Theoretical RC	0.15	0.30	0.45	0.60

Table 5. TABLE V : Approximation errors and Relative Complexities for the FAuST approximations of the MEG matrix.

Mean ${ϵ}_j$	$2.1\cdot {10}^{-1}$	$1.2\cdot {10}^{-1}$	$6.6\cdot {10}^{-2}$	$5.9\cdot {10}^{-2}$
Measured RC	0.22	0.33	0.45	0.63

Table 6. TABLE VI : Screening: number of floating-point operations.

Screening	${𝜽}_t$	$ℛ$	Test	Total
Dynamic Safe	$3N+\left|𝒜\right|$	$N$	$\left|𝒜\right|$	$4N+2\left|𝒜\right|$
Stable Dynamic Safe	$4N+2\left|𝒜\right|$	$N$	$\left|𝒜\right|$	$5N+3\left|𝒜\right|$
GAP Safe	$3N+\left|𝒜\right|$	$N$	$\left|𝒜\right|$	$4N+2\left|𝒜\right|$
Stable GAP Safe	$4N+2\left|𝒜\right|$	$N$	$N+2\left|𝒜\right|$	$6N+4\left|𝒜\right|$

Equations

$$\mathcal{L}(\lambda ,\mathbf{A},\mathrm{y}):$$

$$\lambda \le {\lambda }_{\max }:=\Vert {\mathbf{A}}^T\mathrm{y}{\Vert }_{\infty }$$

$${\mathrm{x}}_{t+1}\leftarrow p({\mathrm{x}}_t,\mathbf{A},L_t)={\mathrm{ST}}_{\frac{\lambda }{L_t}}\left({\mathrm{x}}_t+\frac{1}{L_t}{\mathbf{A}}^T(\mathrm{y}-\mathbf{A}{\mathrm{x}}_t)\right)$$

$$\mathbf{A}=\tilde{\mathbf{A}}+\mathbf{E}$$

$${\theta }^{\star }=$$

$${\Delta }_{\mathbf{A}}$$

$$\bm{\mathrm{a}}^{T}_{j}\bm{\mathrm{\theta}}^{\star}=\left\{\begin{array}[]{lcc}\operatorname{sign}(\bm{\mathrm{x}}^{\star}(j))&\text{if}&\bm{\mathrm{x}}^{\star}(j)\neq 0\\ \gamma\in[-1,1]&\text{if}&\bm{\mathrm{x}}^{\star}(j)=0\end{array}\right.,\quad\forall j\in\{1,\dots,K\}.$$

$$\mu (\mathrm{a}|\mathcal{R})=\underset{\theta \in \mathcal{R}}{\sup }|{\mathrm{a}}^T\theta |$$

$$\mu ({\mathrm{a}}_j|\mathcal{R})<1⟹{\mathrm{x}}^{\star }(j)=0.$$

$$\mu (\mathrm{a}|B_p(\mathrm{c},R))=|{\mathrm{a}}^T\mathrm{c}|+R\Vert \mathrm{a}{\Vert }_{p^{\ast }}$$

$$\mathcal{R}=B(\mathrm{c}=\mathrm{y}/\lambda ,R=\Vert {\theta }_F-\mathrm{y}/\lambda {\Vert }_2),$$

$$\Theta (\mathrm{z}|\mathbf{A}):=\frac{\mathrm{z}}{\max \left(1,\Vert {\mathbf{A}}^T\mathrm{z}{\Vert }_{\infty }\right)}.$$

$$\Theta \left({\mathrm{proj}}_{\mathrm{z}}\frac{\mathrm{y}}{\lambda }|\mathbf{A}\right)=\underset{{\mathrm{z}}^{\prime }=\alpha \mathrm{z}}{\mathrm{argmin}}\Vert {\mathrm{z}}^{\prime }-\mathrm{y}/\lambda {\Vert }_2\text{s.t. }\Vert {\mathbf{A}}^T{\mathrm{z}}^{\prime }{\Vert }_{\infty }\le 1$$

$${\theta }_t=$$

$$\text{Static Safe: }\mathrm{c}=\mathrm{y}/\lambda ,R=|\frac{1}{{\lambda }_{\max }}-\frac{1}{\lambda }|\cdot \Vert \mathrm{y}{\Vert }_2$$

$$\text{Dynamic Safe: }\mathrm{c}=\mathrm{y}/\lambda ,R=\Vert {\theta }_t-\mathrm{y}/\lambda {\Vert }_2.$$

$$G({\mathrm{x}}_t,{\theta }_t|\mathbf{A}):=P({\mathrm{x}}_t|\mathbf{A})-D({\theta }_t)$$

$$\text{GAP Safe: }\mathrm{c}={\theta }_t,R=\frac{1}{\lambda }\sqrt{2G({\mathrm{x}}_t,{\theta }_t|\mathbf{A})}.$$

$$\mu (\mathcal{S}|\mathcal{R}):=\underset{\mathrm{a}\in \mathcal{S}}{\sup }\mu (\mathrm{a}|\mathcal{R})=\underset{\mathrm{a}\in \mathcal{S}}{\sup }\underset{\theta \in \mathcal{R}}{\sup }|{\mathrm{a}}^T\theta |$$

$$\mu (\mathcal{S}|\mathcal{R})<1⟹{\mathrm{x}}^{\star }(j)=0.$$

$$\mathcal{S}:=B_q(\tilde{\mathrm{a}},ϵ)$$

$$\mu (\mathcal{S}|B_p(\mathrm{c},R))\le |{\tilde{\mathrm{a}}}^T\mathrm{c}|+ϵ\Vert \mathrm{c}{\Vert }_{q^{\ast }}+R\Vert \tilde{\mathrm{a}}{\Vert }_{p^{\ast }}+CRϵ.$$

$${\mathcal{S}}^{\prime }:=\mathcal{S}\cap \{\mathrm{a}:\Vert \mathrm{a}{\Vert }_{p^{\ast }}=a\}.$$

$$\mu ({\mathcal{S}}^{\prime }|B_p(\mathrm{c},R))=|{\tilde{\mathrm{a}}}^T\mathrm{c}|+ϵ\Vert \mathrm{c}{\Vert }_{q^{\ast }}+Ra.$$

$${\Theta }^{\prime }(\mathrm{z}|\tilde{\mathbf{A}},ϵ):=\frac{\mathrm{z}}{\max \left(1,{\max }_{1\le j\le K}(|{\tilde{\mathrm{a}}}_j^T\mathrm{z}|+{ϵ}_j\Vert \mathrm{z}{\Vert }_{q^{\ast }})\right)}$$

$${\theta }_t^{\prime }={\left[\frac{{\mathrm{y}}^T{\tilde{\rho }}_t}{\lambda \Vert {\tilde{\rho }}_t{\Vert }_2^2}\right]}_{-{\alpha }_t^{\prime }}^{{\alpha }_t^{\prime }}{\tilde{\rho }}_t,$$

$$\text{with }{\alpha }_t^{\prime }={\left(\underset{1\le j\le K}{\max }\left(|{\tilde{\mathrm{a}}}_j^T{\tilde{\rho }}_t|+{ϵ}_j\Vert {\tilde{\rho }}_t{\Vert }_{q^{\ast }}\right)\right)}^{-1}$$

$${\lambda }_{\max }^{\prime }:=\underset{j}{\max }\left(|{\tilde{\mathrm{a}}}_j^T\mathrm{y}|+{ϵ}_j\Vert \mathrm{y}{\Vert }_{q^{\ast }}\right)\ge {\lambda }_{\max },$$

$$\text{Stable Static Safe: }{\mathrm{c}}^{\prime }=\mathrm{y}/\lambda ,R^{\prime }=|\frac{1}{{\lambda }_{\max }^{\prime }}-\frac{1}{\lambda }|\cdot \Vert \mathrm{y}{\Vert }_2.$$

$$\text{Stable Dynamic Safe:}{\mathrm{c}}^{\prime }=\mathrm{y}/\lambda ,R^{\prime }=\Vert {\theta }_t^{\prime }-\mathrm{y}/\lambda {\Vert }_2.$$

$$\mathcal{E}\ge \Vert \mathbf{A}-\tilde{\mathbf{A}}{\Vert }_{r\to 2}$$

Full text

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 [2], source localization [3] or spike deconvolution [4], 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 [5] [6], and convex optimization approaches relying on $\ell_{1}$-norm minimization, where the regression problem is expressed as an $\ell_{1}$ regularized optimization problem –known as Basis Pursuit [7] or LASSO [8]– which solution is computed using iterative convex optimization algorithms such as ISTA or its variants [9, 10, 11, 12, 13].

For very high-dimensional problems, however, iterative algorithms to solve $\ell_{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:

Fast structured operators [14, 15, 16, 17, 18] which provide faster matrix-vector products (see Section III-C); 2. 2.

Safe screening tests [19, 20, 21, 22, 23, 24], which safely eliminate unused explanatory variables (see Section II).

This paper extends the results in [25] [26]. In [25] we introduced a safe screening sphere test that manipulates an approximate dictionary and in [26] the GAP Safe sphere test [22] was extended to this new setting. Here, we revisit the previous results under a broader formalism and propose a fast algorithm for $\ell_{1}$-minimization which combines safe screening and (potentially multiple) fast approximate dictionaries.

I-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 $\bm{\mathrm{y}}\in\mathbb{R}^{N}$ and modeled as $\bm{\mathrm{y}}\approx\mathbf{A}\bm{\mathrm{x}}_{0}$ where $\bm{\mathrm{x}}_{0}\in\mathbb{R}^{K}$ is sparse and ${\mathbf{A}=[\bm{\mathrm{a}}_{1},\dots,\bm{\mathrm{a}}_{K}]}\in\mathbb{R}^{N\times K}$. In the context of linear inverse problems, $\mathbf{A}$ is the so-called measurement matrix [27], while for sparse signal representations $\mathbf{A}$ would be the dictionary matrix [7], and in statistics the design matrix [8]. In this paper, we adopt the terminology of sparse signal representations hence $\mathbf{A}$ is called a dictionary and its columns, denoted $\bm{\mathrm{a}}_{j}\in\mathbb{R}^{N}$, are called atoms.

The $\ell_{1}$-regularized least squares, referred to as Lasso or Basis Pursuit, consists in finding a sparse coefficient vector $\bm{\mathrm{x}}\in\mathbb{R}^{K}$, solution of the following optimization problem:

[TABLE]

where $P(\bm{\mathrm{x}}|\mathbf{A})$ is called the primal objective and the parameter $\lambda\!>\!0$ controls the trade-off between data fidelity and sparsity of the solution. We suppose that

[TABLE]

since otherwise $\bm{\mathrm{x}}^{\star}=\bm{\mathrm{0}}\in\mathbb{R}^{K}$ is the unique solution.

I-B Iterative algorithms

First-order iterative algorithms, especially proximal-gradient methods (ISTA [9], FISTA [10], TwIST [11], SpaRSA [12], Chambolle-Pock [13]), are popular approaches for solving (1).

One can use an abstract notations to represent the update step for a generic iterative algorithm for the problem $\mathcal{L}(\lambda,\mathbf{A},\bm{\mathrm{y}})$: $\{\bm{\mathrm{x}}_{t+1},\bm{\mathrm{\alpha}}_{t+1}\}\leftarrow p(\mathbf{x}_{t},\mathbf{A},\bm{\mathrm{\alpha}}_{t})$ where $\mathbf{x}_{t}$ is the current estimate of the primal variable and $\bm{\mathrm{\alpha}}_{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

[TABLE]

where $\operatorname{ST}_{u}(z)\!=\!{\operatorname{sign}(z)(|z|-u)_{{}+{}}}$ denotes the soft-thresholding operation, which is the proximal operator associated to the $\ell_{1}$-norm.

I-C Computational bottlenecks

The main bottleneck of existing iterative algorithms in terms of computational complexity is the cost of the required matrix-vector 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 $\mathbf{A}^{T}\mathbf{A}$ is pre-computed and stored).

I-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 $\bm{\mathrm{x}}^{\star}$ before even computing it.

I-D1 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 [14], 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 [17] 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 [16]. 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 sub-dictionaries of smaller size. It was proposed in [28] 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 [18] and [29] (respectively for 2-dimensional and higher-order tensorial data), increasing its flexibility in exchange of some additional complexity.

In practical applications, the dictionary matrix $\mathbf{A}$ is not necessarily structured. A possible strategy is to replace certain iterations $\bm{\mathrm{x}}_{t+1}=p(\mathbf{x}_{t},\mathbf{A},\bm{\mathrm{\alpha}}_{t})$ with $\bm{\mathrm{x}}_{t+1}=p(\bm{\mathrm{x}}_{t},\tilde{\mathbf{A}},\bm{\mathrm{\alpha}}_{t})$ where $\tilde{\mathbf{A}}$ is a structured approximation of $\mathbf{A}$, i.e.,

[TABLE]

with approximation error $\mathbf{E}\!=\![\bm{\mathrm{e}}_{1},\dots,\bm{\mathrm{e}}_{K}]\!\in\!\mathbb{R}^{N\times K}$.

I-D2 Acceleration with screening rules

For a given instance $\mathcal{L}(\lambda,\mathbf{A},\bm{\mathrm{y}})$ of (1), characterized by a regularization parameter $\lambda$, a dictionary $\mathbf{A}$, and an input vector $\bm{\mathrm{y}}$, the safe screening approach, introduced by [19], consists in identifying and removing from the dictionary a subset of atoms which are guaranteed to have zero weight in a solution $\bm{\mathrm{x}}^{\star}$, 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 [30, 31] based on correlation measures between the atoms $\bm{\mathrm{a}}_{j}$ and the input signal $\bm{\mathrm{y}}$, safe screening has zero risk of false rejections. Moreover, screening techniques are transparent to the underlying $\ell_{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 [19] [20] [32], 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 [22], are repeatedly applied during the iterative optimization algorithm leveraging its current solution estimate and gradually increasing the set of rejected atoms.

Sequential rules [33, 34, 20, 35] exploit the fact that (1) is commonly solved over a grid of regularization parameters $\lambda_{i}$ and reuse the results from a previous configuration $\lambda_{i-1}$ to improve the screening performance.

I-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 $\ell_{p}$-balls, which allows us to exemplify the stable tests on some existing $\ell_{2}$-sphere tests –one of which being the state-of-the-art GAP Safe [22]– for a general $\ell_{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 $\ell_{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.

I-F Related Work

Apart from the aforementioned structured dictionaries and safe screening tests –as well as other preceding correlation-based feature selection heuristics [30, 31]– 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 [36, 37] start with small restricted problems and progressively include more promising features. In [38], the authors combine a working set strategy with safe screening and in [39] 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 [40] 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.

I-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 $\mathcal{A}$ of non-repeating integers, ${\mathbf{A}_{[\mathcal{A}]}:=[\bm{\mathrm{a}}_{i}]_{i\in\mathcal{A}}}$ denotes a sub-matrix of $\mathbf{A}$ composed of the columns indexed by the elements in $\mathcal{A}$. The notation extends to vectors: $\bm{\mathrm{x}}_{[\mathcal{A}]}:=[\bm{\mathrm{x}}(i)]_{i\in\mathcal{A}}$, where $\bm{\mathrm{x}}(i)$ denotes the $i$-th component of $\bm{\mathrm{x}}$. Screening rules allow to identify such sets of reduced cardinality $|\mathcal{A}_{t}|<K$, so that the complexity of applying $\mathbf{A}_{[\mathcal{A}]}$ or its transpose is reduced compared to that of applying $\mathbf{A}$, while ensuring that $\mathcal{A}_{t}$ contains the support of the solution $\{i:\bm{\mathrm{x}}^{\star}(i)\neq 0\}$. This is achieved using the dual of (1) [19]:

[TABLE]

where

[TABLE]

is the dual feasible set and $D(\bm{\mathrm{\theta}})$ is the dual objective. The dual and primal solutions ($\bm{\mathrm{\theta}}^{\star}$ and $\bm{\mathrm{x}}^{\star}$) are linked through the relation $\bm{\mathrm{y}}=\mathbf{A}\bm{\mathrm{x}}^{\star}+\lambda\bm{\mathrm{\theta}}^{\star}$. Optimality conditions (KKT) at the dual solution $\bm{\mathrm{\theta}}^{\star}$ read (see for instance [24] for more details)

[TABLE]

Hence, every dictionary atom for which $|\bm{\mathrm{a}}_{j}^{T}\bm{\mathrm{\theta}}^{\star}|\!<\!1$ is inactive.

II-A Notion of safe region

Since the optimal solution $\bm{\mathrm{\theta}}^{\star}$ of the dual problem (5) is not known, the inner products $\bm{\mathrm{a}}_{j}^{T}\bm{\mathrm{\theta}}^{\star}$ cannot be evaluated. Fortunately, given only $\mathbf{A}$ and $\bm{\mathrm{y}}$, it is possible to identify at a moderate computational cost a region $\mathcal{R}\subset\mathbb{R}^{N}$, called safe region, which is guaranteed to contain the optimal $\bm{\mathrm{\theta}}^{\star}$.

Definition 1** (Safe region).**

A region $\mathcal{R}\subset\mathbb{R}^{N}$ is safe (with respect to the dual problem (5)) if and only if it contains the dual solution $\bm{\mathrm{\theta}}^{\star}$, i.e. $\bm{\mathrm{\theta}}^{\star}\in\mathcal{R}$.

II-B Screening an atom given a safe region

Consider an atom $\bm{\mathrm{a}}_{j}$. If the inequality $|\bm{\mathrm{a}}_{j}^{T}\bm{\mathrm{\theta}}|<1$ holds for all $\bm{\mathrm{\theta}}\in\mathcal{R}$ where $\mathcal{R}$ is a safe region, then the above analysis ensures that $\bm{\mathrm{a}}_{j}$ is inactive. This gives rise to a screening test.

Definition 2** (Screening test for an atom).**

Given a region $\mathcal{R}$, a screening test for the atom $\bm{\mathrm{a}}\in\mathbb{R}^{N}$ is given by:

[TABLE]

A sufficient condition for an atom $\bm{\mathrm{a}}_{j}$ to be inactive can be expressed as follows: if $\mathcal{R}$ is safe then

[TABLE]

In practice, for each atom $\bm{\mathrm{a}}_{j}$, computing the test $\mu(\bm{\mathrm{a}}_{j}|\mathcal{R})$ allows to eliminate or not the atom.

Formally, given a safe region, the atoms can be partitioned into a preserved set $\mathcal{A}=\{j\in\{1,\dots,K\}:\mu(\bm{\mathrm{a}}_{j}|\mathcal{R})\geq 1\}$ and its complement, the rejection set $\mathcal{A^{\mathsf{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 $\mathbf{A}$, $\bm{\mathrm{y}}$, and possibly the current iterate $\bm{\mathrm{x}}_{t}$ of an iterative algorithm addressing (1).

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 $\mathcal{R}$ which is a closed $\ell_{p}$-ball with center $\bm{\mathrm{c}}$ and radius $R$, denoted $B_{p}(\bm{\mathrm{c}},R)=\{\bm{\mathrm{z}}:\|\bm{\mathrm{z}}-\bm{\mathrm{c}}\|_{p}\leq R\}$. The screening test for this region has a closed form (a proof is given in Appendix -A)

[TABLE]

where $\|\cdot\|_{p^{*}}$ denotes the dual norm associated to the $p$-norm, with $\frac{1}{p}+\frac{1}{p^{*}}=1$. For simplicity, we omit subscripts for the $\ell_{2}$-ball ($p=p^{*}=2$) denoted $B(\bm{\mathrm{c}},R)$.

II-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 $\|\bm{\mathrm{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 $\bm{\mathrm{a}}^{T}\bm{\mathrm{c}}$ for all atoms that have not been screened out yet. This calls for techniques where the choice of $R$ and $\bm{\mathrm{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 [19], follows from the simple observation that the solution $\bm{\mathrm{\theta}}^{\star}$ of the dual problem (5) is the Euclidean projection of $\bm{\mathrm{y}}/\lambda$ on the feasible set $\Delta_{\mathbf{A}}$. As a result, if some feasible point $\bm{\mathrm{\theta}}_{F}\in\Delta_{\mathbf{A}}$ is known, then $\bm{\mathrm{\theta}}^{\star}$ can’t be further away from $\bm{\mathrm{y}}/\lambda$ than $\bm{\mathrm{\theta}}_{F}$ in the $\ell_{2}$ sense, i.e. $\|\bm{\mathrm{\theta}}^{\star}-\bm{\mathrm{y}}/\lambda\|_{2}\leq\|\bm{\mathrm{\theta}}_{F}-\bm{\mathrm{y}}/\lambda\|_{2}~{}\forall\bm{\mathrm{\theta}}_{F}\in\Delta_{\mathbf{A}}$.

This leads to an $\ell_{2}$-spherical safe region ($p\!=\!p^{*}\!=\!2$),

[TABLE]

which full specification requires a feasible point $\bm{\mathrm{\theta}}_{F}\in\Delta_{\mathbf{A}}$.

To generate a feasible point $\bm{\mathrm{\theta}}_{F}\in\Delta_{\mathbf{A}}$, one could compute the Euclidean projection of any given point $\bm{\mathrm{z}}\in\mathbb{R}^{N}$ over the closed convex set $\Delta_{\mathbf{A}}$. As this is too computationally demanding, simpler scaling approaches are preferred: given $\bm{\mathrm{z}}\in\mathbb{R^{N}}$ we determine a scaling factor $\alpha$ so that $\Theta(\bm{\mathrm{z}}|\mathbf{A}):=\alpha\bm{\mathrm{z}}\in\Delta_{\mathbf{A}}$.

By definition of $\Delta_{\mathbf{A}}$, if $\|\mathbf{A}^{T}\bm{\mathrm{z}}\|_{\infty}\leq 1$ no scaling is needed, otherwise one only needs to divide $\bm{\mathrm{z}}$ by $\|\mathbf{A}^{T}\bm{\mathrm{z}}\|_{\infty}$. This yields

[TABLE]

While the scaling in 10 selects a feasible point to build a safe sphere (9) given an arbitrary $\bm{\mathrm{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 [19] formalized as follows (a proof is given in [21, Lemma 8])

Proposition 1**.**

Denoting $\mathrm{proj}_{\bm{\mathrm{z}}}\dfrac{\bm{\mathrm{y}}}{\lambda}=\dfrac{\bm{\mathrm{y}}^{T}\bm{\mathrm{z}}}{\lambda\|\bm{\mathrm{z}}\|_{2}^{2}}\bm{\mathrm{z}}$ the Euclidean projection of $\bm{\mathrm{y}}/\lambda$ onto the direction of $\bm{\mathrm{z}}$, we have

[TABLE]

i.e. among all points which are both dual feasible and proportional to $\bm{\mathrm{z}}$, $\Theta\left(\mathrm{proj}_{\bm{\mathrm{z}}}\dfrac{\bm{\mathrm{y}}}{\lambda}~{}|~{}\mathbf{A}\right)$ is the one closest to $\bm{\mathrm{y}}/\lambda$.

To obtain dynamic spherical regions, we need dual feasible points that change over iterations ($\bm{\mathrm{\theta}}_{F}=\bm{\mathrm{\theta}}_{t}$ at iteration $t$) and leverage the current state of the optimization algorithm. Following [21], this can be achieved by defining $\bm{\mathrm{\theta}}_{t}$ proportional to the current residual $\bm{\mathrm{\rho}}_{t}:=\bm{\mathrm{y}}-\mathbf{A}\bm{\mathrm{x}}_{t}$, since the dual solution $\bm{\mathrm{\theta}}^{\star}$ is proportional to $\bm{\mathrm{y}}-\mathbf{A}\bm{\mathrm{x}}^{\star}$. Taking $\bm{\mathrm{z}}=\bm{\mathrm{\rho}}_{t}$ and $\bm{\mathrm{\theta}}_{t}=\Theta\left(\mathrm{proj}_{\bm{\mathrm{\rho}}_{t}}\dfrac{\bm{\mathrm{y}}}{\lambda}~{}|~{}\mathbf{A}\right)$ yields

[TABLE]

where $[z]_{a}^{b}:=\min(\max(z,a),b)$ denotes the projection of the scalar $z$ onto the segment $[a,b]$.

The above principles yield existing safe spherical regions.

II-C1 Static Safe sphere [19]

With $\bm{\mathrm{z}}=\bm{\mathrm{y}}/\lambda$ in (10) one gets $\theta_{F}=\Theta(\bm{\mathrm{y}}/\lambda|\mathbf{A})=\bm{\mathrm{y}}/\lambda_{\max}$, and (9) reads $\mathcal{R}=B(\bm{\mathrm{c}},R)$ with

[TABLE]

where we recall that $\lambda_{\max}$ is defined in (2).

II-C2 Dynamic Safe sphere [21]

With the dual feasible point defined in (11), (9) yields $\mathcal{R}=B(\bm{\mathrm{c}},R)$ with

[TABLE]

II-C3 GAP Safe sphere [22]

Given any primal-dual feasible pair $(\bm{\mathrm{x}}_{t},\bm{\mathrm{\theta}}_{t})\in\mathbb{R}^{K}\!\times\!\Delta_{\mathbf{A}}$ (say the solution estimations at iteration $t$) and denoting

[TABLE]

the corresponding duality gap, the GAP Safe sphere is $\mathcal{R}=B(\bm{\mathrm{c}},R)$ with

[TABLE]

This region is provably safe for any dual feasible point $\bm{\mathrm{\theta}}_{t}\!\in\!\Delta_{\mathbf{A}}$ [22, Theorem 2]. The authors, following [21], propose to use $\bm{\mathrm{\theta}}_{t}$ defined by (11).

Note that the cited safe regions reuse quantities computed along the iterations of classical proximal algorithms such as the residual $\bm{\mathrm{\rho}}_{t}$, the inner products $\mathbf{A}^{T}\bm{\mathrm{\rho}}_{t}$, or the duality gap (which typically serves as a stopping criterion).

III Stable Safe Screening

When iterating with $\tilde{\mathbf{A}}$ instead of $\mathbf{A}$, the computational overhead of accessing inner products $\bm{\mathrm{a}}^{T}\bm{\mathrm{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 $\mathbf{A}$ (instead of $\tilde{\mathbf{A}}$) to implement the above screening tests, but the whole purpose of re-using quantities already computed with a structured (computationally efficient) matrix $\tilde{\mathbf{A}}$ 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 $\tilde{\mathbf{A}}$ and $\mathbf{A}$.

III-A Screening a zone given a safe region

In conventional screening tests (Definition 2) the dual solution $\bm{\mathrm{\theta}}^{\star}$ is assumed to belong to a known safe region $\mathcal{R}$. This allows to screen (or not) an atom $\bm{\mathrm{a}}$ that is perfectly known, (re)using computations of inner products with this atom.

When iterating with $\tilde{\mathbf{A}}$, efficient screening tests are only entitled to reuse inner products with approximated atoms, $\tilde{\mathbf{A}}^{T}\bm{\mathrm{c}}$. Yet, they should be able to screen certain atoms $\bm{\mathrm{a}}$, using the only knowledge that these are “close” to the corresponding approximated atoms $\tilde{\bm{\mathrm{a}}}$.

One way to capture this knowledge is to consider $\mathcal{S}\subset\mathbb{R}^{N}$ some neighborhood of $\bm{\mathrm{a}}$, assumed to be known. We will call $\mathcal{S}$ a zone. This gives rise to screening tests for zones.

Definition 3** (Screening test for a zone).**

Given a safe region $\mathcal{R}$, a screening test for the zone $\mathcal{S}\subset\mathbb{R}^{N}$ is given by:

[TABLE]

A sufficient condition for an atom $\bm{\mathrm{a}}_{j}$ to be inactive can be expressed as follows: if $\mathcal{R}$ is safe and if $\bm{\mathrm{a}}_{j}\in\mathcal{S}$, then

[TABLE]

III-B Stable* screening tests

Given an approximate dictionary $\tilde{\mathbf{A}}$, error bounds $\epsilon_{j}$ such that $\|\tilde{\bm{\mathrm{a}}}_{j}-\bm{\mathrm{a}}_{j}\|_{q}\leq\epsilon_{j}$, $1\leq j\leq K$, can be pre-computed for some choice of $q\in[1,\infty]$, yielding spherical zones $\mathcal{S}_{j}:=B_{q}(\tilde{\bm{\mathrm{a}}}_{j},\epsilon_{j})$ known to contain the atoms: $\bm{\mathrm{a}}_{j}\in\mathcal{S}_{j}$. Using (8), a spherical zone

[TABLE]

gives rise to the following estimate:

[TABLE]

where $C=C_{p,q}:=N^{\left(1-\frac{1}{p}-\frac{1}{q}\right)_{+}}$ and $q^{*}\in[1,\infty]$ is such that $\frac{1}{q}+\frac{1}{q^{*}}=1$. See proof in Appendix -B.

In practice, the norms $a_{j}\!=\!\|\bm{\mathrm{a}}_{j}\|_{p^{*}}$ can be pre-computed and stored with negligible overhead, yielding restricted spherical zones

[TABLE]

For such zones we get

[TABLE]

Lemma 1**.**

Consider a dictionary $\mathbf{A}$ and $\mathcal{R}=B_{p}(\mathbf{c},R)$ a safe region with respect to the dual problem (5). Let $\tilde{\mathbf{A}}$ be an approximate dictionary, $1\leq j\leq K$, and $\mathcal{S}_{j}$ and $\mathcal{S}^{\prime}_{j}$ defined by (17)-(20) with $\tilde{\bm{\mathrm{a}}}=\tilde{\bm{\mathrm{a}}}_{j}$, $\epsilon\geq\|\tilde{\bm{\mathrm{a}}}_{j}-\bm{\mathrm{a}}_{j}\|_{q}$, and $a=\|\bm{\mathrm{a}}_{j}\|_{p^{*}}$.

• •

If $\mu(\mathcal{S}_{j}|B_{p}(\mathbf{c},R))<1$ then $\bm{\mathrm{x}}_{j}^{\star}=0$;

• •

If $\mu(\mathcal{S}^{\prime}_{j}|B_{p}(\mathbf{c},R))<1$ then $\bm{\mathrm{x}}_{j}^{\star}=0$.

Remark 1**.**

In [40], 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 $\mathcal{S}$, reducing the total number of tests performed.

III-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 $\mathbf{A}$, to define new safe regions reusing computations done during the iterations with $\tilde{\mathbf{A}}$ instead of $\mathbf{A}$.

Remark 2**.**

Let us emphasize that a safe region is a region $\mathcal{R}$ containing the dual solution $\bm{\mathrm{\theta}}^{\star}$ of the original dual problem (5) and not necessarily the dual solution of its approximate version with $\Delta_{\tilde{\mathbf{A}}}$ instead of $\Delta_{\mathbf{A}}$. Indeed, although $\mathbf{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 $\bm{\mathrm{\theta}}_{t}\!\in\!\Delta_{\mathbf{A}}$. Given an arbitrary $\bm{\mathrm{z}}\in\mathbb{R}^{N}$, the function $\Theta(\bm{\mathrm{z}}|\mathbf{A})$ from (10) would provide such a feasible point, however it cannot be used as its computation requires multiplication by $\mathbf{A}^{T}$. A naive alternative could be to compute $\Theta(\bm{\mathrm{z}}|\tilde{\mathbf{A}})$, however while this always belongs to $\Delta_{\tilde{\mathbf{A}}}$ (the feasible set for $\tilde{\mathbf{A}}$) it does not necessarily belong to the desired feasible set $\Delta_{\mathbf{A}}$. This can be fixed using a modified dual scaling approach that we propose to call stable dual scaling.

Considering error bounds $\bm{\mathrm{\epsilon}}=(\epsilon_{j})_{j=1}^{K}\in\mathbb{R}_{+}^{K}$ we define

[TABLE]

For $\bm{\mathrm{\epsilon}}=\bm{\mathrm{0}}$ we recover $\Theta(\bm{\mathrm{z}}|\mathbf{A})=\Theta^{\prime}(\bm{\mathrm{z}}|\mathbf{A},\bm{\mathrm{0}})$.

Lemma 2**.**

Assume that $\bm{\mathrm{a}}_{j}\in\mathcal{S}_{j}=B_{q}(\tilde{\bm{\mathrm{a}}}_{j},\epsilon_{j})$ for $1\leq j\leq K$. Then, for any $\bm{\mathrm{z}}\in\mathbb{R}^{N}$, $\Theta^{\prime}(\bm{\mathrm{z}}|\tilde{\mathbf{A}},\bm{\mathrm{\epsilon}})\in\Delta_{\mathbf{A}}\cap\Delta_{\tilde{\mathbf{A}}}$ is a feasible point w.r.t. both the original dual problem (5) and its modified version with $\tilde{\mathbf{A}}$ instead of $\mathbf{A}$.

The proof is in Appendix -C. Analogously to (11), we can now define a dual feasible point proportional to the residual $\tilde{\bm{\mathrm{\rho}}}_{t}\!:=\!\bm{\mathrm{y}}\!-\!\tilde{\mathbf{A}}\bm{\mathrm{x}}_{t}$ at iteration $t$, $\bm{\mathrm{\theta}}^{\prime}_{t}=\Theta^{\prime}\left(\mathrm{proj}_{\tilde{\bm{\mathrm{\rho}}}_{t}}\dfrac{\bm{\mathrm{y}}}{\lambda}~{}|~{}\tilde{\mathbf{A}},\bm{\mathrm{\epsilon}}\right)$, that is to say more explicitly

[TABLE]

These principles yield new stable safe $\ell_{2}$-spherical regions.

III-C1 Stable Static Safe sphere

The static safe sphere (12) depends on $\lambda_{\max}=\max_{j}|\bm{\mathrm{a}}_{j}^{T}\bm{\mathrm{y}}|$. To perform only computations with $\tilde{\mathbf{A}}$, we can reuse (12) with

[TABLE]

leading to $\mathcal{R}=B(\bm{\mathrm{c}}^{\prime},R^{\prime})$ with

[TABLE]

III-C2 Stable Dynamic Safe sphere

We adapt the Dynamic Safe sphere in (13) to our approximate setting by replacing the dual feasible point $\bm{\mathrm{\theta}}_{t}$ (11) by $\bm{\mathrm{\theta}}^{\prime}_{t}$ (III-C), leading to $\mathcal{R}=B(\bm{\mathrm{c}}^{\prime},R^{\prime})$ with

[TABLE]

The fact that $\bm{\mathrm{\theta}}^{\prime}_{t}\in\Delta_{\mathbf{A}}$ directly implies that this sphere is safe by definition of the basic $\ell_{2}$-spherical bound (9).

III-C3 Stable GAP Safe sphere

To build a safe sphere $B(\bm{\mathrm{c}}^{\prime},R^{\prime})$ with $\bm{\mathrm{c}}^{\prime}=\bm{\mathrm{\theta}}^{\prime}_{t}$ (III-C) instead of $\bm{\mathrm{c}}=\bm{\mathrm{\theta}}_{t}$ (11) we need to determine a safe radius $R^{\prime}$ depending only on $\tilde{\mathbf{A}}$. For this reason, instead of using the (inaccessible) duality gap $G(\bm{\mathrm{x}}_{t},\bm{\mathrm{\theta}}^{\prime}_{t}|\mathbf{A})=P(\bm{\mathrm{x}}_{t}|\mathbf{A})-D(\bm{\mathrm{\theta}}^{\prime}_{t})$ we will use $G(\bm{\mathrm{x}}_{t},\bm{\mathrm{\theta}}^{\prime}_{t}|\tilde{\mathbf{A}})=P(\bm{\mathrm{x}}_{t}|\tilde{\mathbf{A}})-D(\bm{\mathrm{\theta}}^{\prime}_{t})$ and add a security margin ensuring that $B(\bm{\mathrm{\theta}}^{\prime}_{t},R^{\prime})$ is safe. Recall the standard notation $\|\mathbf{M}\|_{p\to q}:=\sup_{\|u\|_{p}\leq 1}\|\mathbf{M}\bm{\mathrm{u}}\|_{q}$.

Theorem 1** (Stable GAP Safe sphere).**

Consider

[TABLE]

where $1\leq r\leq\infty$. For any111In this theorem, $\theta\in\Delta_{\mathbf{A}}$ can be arbitrary. In particular it is not necessarily expressed as in (III-C), hence it does not need to belong to $\Delta_{\tilde{\mathbf{A}}}$. $(\bm{\mathrm{x}},\bm{\mathrm{\theta}})\in\mathbb{R}^{K}\!\times\!\Delta_{\mathbf{A}}$, we have

[TABLE]

In other words $B(\bm{\mathrm{\theta}},R^{\prime})$ is safe.

The proof was given in [26]. By applying Theorem (1) to $\bm{\mathrm{x}}_{t}$ and $\bm{\mathrm{\theta}}^{\prime}_{t}$ in (III-C) we obtain $\mathcal{R}=B(\bm{\mathrm{c}}^{\prime},R^{\prime})$ with

[TABLE]

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 $\tilde{\mathbf{A}}$, making them suitable in practice.

IV A Fast Algorithm for $\ell_{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 $\bm{\mathrm{y}}$ and the regularization parameter $\lambda$, the input of the algorithm consists in a sequence $\{\tilde{\mathbf{A}}^{i}\}_{i=0}^{I}$ of approximate dictionaries with $\tilde{\mathbf{A}}^{I}=\mathbf{A}$, and a corresponding sequence of error bounds $\{\bm{\mathrm{\epsilon}}^{i}\}_{i=0}^{I}$, where $\bm{\mathrm{\epsilon}}^{i}=(\epsilon^{i}_{j})_{j=1}^{K}\in\mathbb{R}_{+}^{K}$ contains $\ell_{2}$ error bounds on the atoms of the $i$-th dictionary $\tilde{\mathbf{A}}^{i}$, and of course $\bm{\mathrm{\epsilon}}^{I}=\bm{\mathrm{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(\bm{\mathrm{x}}_{t},\mathbf{A},\bm{\mathrm{\alpha}}_{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 $\Gamma$. In the following subsections, we present a switching criterion which assumes the approximations $\{\tilde{\mathbf{A}}\}_{i=0}^{I}$ 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(\bm{\mathrm{x}}^{\star},\bm{\mathrm{\theta}}^{\star}|\mathbf{A})\!=\!0$ [41]. The simplest way of computing the duality gap is to wait until the algorithm reaches $\tilde{\mathbf{A}}=\tilde{\mathbf{A}}^{I}=\mathbf{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 $\bm{\mathrm{\epsilon}}^{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 ($\lambda/\lambda_{\max}\!\approx\!1$) in which screening typically occurs very quickly.

IV-A Switching: convergence criterion

In [25, 26], we proposed a empirically-tuned heuristic convergence criterion. In the following, we derive a more theoretically grounded criterion.

Given a solution estimate $\bm{\mathrm{x}}_{t}$ at iteration $t$, a natural convergence measure for the problem $\mathcal{L}(\lambda,\mathbf{A},\bm{\mathrm{y}})$ is the duality gap $G(\bm{\mathrm{x}}_{t},\bm{\mathrm{\theta}}_{t}|\mathbf{A})$ (14) with $\bm{\mathrm{\theta}}_{t}$ calculated as in (11). However, when manipulating approximate dictionaries, this quantity is not accessible without unwanted additional computations.

As a surrogate, we propose to use a ratio

[TABLE]

between two computable gaps arising from two different dual points $\bm{\mathrm{\theta}}^{\prime}_{t}$ and $\tilde{\bm{\mathrm{\theta}}}_{t}$ that are respectively feasible for problems $\mathcal{L}(\lambda,\mathbf{A},\bm{\mathrm{y}})$ and $\mathcal{L}(\lambda,\tilde{\mathbf{A}},\bm{\mathrm{y}})$:

[TABLE]

where $\bm{\mathrm{\theta}}^{\prime}_{t}\in\Delta_{\mathbf{A}}\cap\Delta_{\tilde{\mathbf{A}}}$ has been previously calculated for the stable screening with (III-C), while $\tilde{\bm{\mathrm{\theta}}}_{t}\in\Delta_{\tilde{\mathbf{A}}}$ is the conventional dual point (11) but calculated with $\tilde{\mathbf{A}}$. Then again, apart from minor extra memory requirements, the computation of the quantities in (31) reuses most of the calculations333Only $N$ extra products for computing $\|\tilde{\bm{\mathrm{\theta}}}\|$ and $|\mathcal{A}_{t}|$ comparisons for $\|\tilde{\mathbf{A}}^{T}\bm{\mathrm{\rho}}_{t}\|_{\infty}$ are required. in (30).

If $\tilde{\mathbf{A}}=\tilde{\mathbf{A}}^{i}$ was kept until convergence, we would have

[TABLE]

with $\tilde{\bm{\mathrm{x}}}^{\star}$ a solution of $\mathcal{L}(\lambda,\tilde{\mathbf{A}},\bm{\mathrm{y}})$ and $\tilde{\bm{\mathrm{\theta}}}^{\star}$ a solution of the corresponding dual problem. Then, as illustrated in Fig. 1 (left), the second quantity $G(\bm{\mathrm{x}}_{t},\tilde{\bm{\mathrm{\theta}}}_{t}|\tilde{\mathbf{A}})$ –which is precisely the duality gap w.r.t. the problem $\mathcal{L}(\lambda,\tilde{\mathbf{A}},\bm{\mathrm{y}})$– would tend to zero, while the first quantity $G(\bm{\mathrm{x}}_{t},\bm{\mathrm{\theta}}^{\prime}_{t}|\tilde{\mathbf{A}})$ would typically saturate at a nonzero value since $\bm{\mathrm{\theta}}^{\prime}_{t}\nrightarrow\tilde{\bm{\mathrm{\theta}}}^{\star}$.

As a result, the ratio $\gamma_{t}$ would shrink with the iterations as illustrated in Fig. 1 (right).

With that in mind, we propose to switch dictionaries when

[TABLE]

where the higher $\Gamma$, the greater sensibility to the gap saturation.

IV-B Speed criterion

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

Denoting $\mathcal{C}(\mathbf{A},\mathcal{A})$ the computational cost of multiplying matrix $\mathbf{A}_{[\mathcal{A}]}$ by a vector, it is reasonable to assume that $\mathcal{C}(\mathbf{A},\mathcal{A})\leq\mathcal{C}(\mathbf{A},\mathcal{A}^{\prime})$ if $\mathcal{A}\subset\mathcal{A}^{\prime}$. As a result, given the nested structure of the sequence of preserved sets $\mathcal{A}_{t}$, the sequence $\mathcal{C}(\mathbf{A},\mathcal{A}_{t})$, $t\geq 0$ is non-increasing.

Consider two approximate dictionaries $\tilde{\mathbf{A}}^{i}$ and $\tilde{\mathbf{A}}^{j}$ and assume the initial “unscreened” computational costs satisfy $\mathcal{C}(\tilde{\mathbf{A}}^{i},\mathcal{A}_{0})<\mathcal{C}(\tilde{\mathbf{A}}^{j},\mathcal{A}_{0})$. Denote $\mathbf{E}^{i}$, $\mathbf{E}^{j}$ the corresponding approximation errors (cf (4)), and assume that $\mathbf{E}^{i}$ is “larger” than $\mathbf{E}^{j}$ (let us denote this $\mathbf{E}^{i}\succ\mathbf{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 $\tilde{\mathbf{A}}^{i}$ becomes more complex than $\tilde{\mathbf{A}}^{j}$

[TABLE]

Note that this only happens if $\mathcal{C}(\tilde{\mathbf{A}}^{i},\mathcal{A}_{t})$ decreases faster than $\mathcal{C}(\tilde{\mathbf{A}}^{j},\mathcal{A}_{t})$ as screening progresses (i.e. as $|\mathcal{A}_{t}|$ gets smaller).

To use the general idea just described, one needs to specify a complexity model $\mathcal{C}(\mathbf{A},\mathcal{A}_{t})$ depending on the structure of the corresponding matrix. For a generic unstructured dictionary matrix $\mathbf{A}\in\mathbb{R}^{N\times K}$, a simple model is:

[TABLE]

To quantify the complexity reduction on matrix-vector products entailed by an approximate dictionary we will use the concept of Relative Complexity (RC) [16], such that

[TABLE]

In a worst-case scenario, screening does not further reduce the cost and $\mathcal{C}(\tilde{\mathbf{A}},\mathcal{A}_{t})=\text{RC}(\tilde{\mathbf{A}})\!\times\!NK$. In more optimistic scenarios, some fast approximate dictionary structures might still benefit from speedups upon column removal.

Consider approximate dictionaries $\{\tilde{\mathbf{A}}^{0},\dots,\tilde{\mathbf{A}}^{I}\}$, with $\tilde{\mathbf{A}}^{I}=\mathbf{A}$, of increasing precision and decreasing complexity

[TABLE]

The above observations suggest to switch from $\tilde{\mathbf{A}}^{i}$ to $\mathbf{A}$ if $\text{RC}(\tilde{\mathbf{A}}^{i})\!\times\!NK=\mathcal{C}(\tilde{\mathbf{A}}^{i},\mathcal{A}_{t})\geq\mathcal{C}(\mathbf{A},\mathcal{A}_{t})=N|\mathcal{A}_{t}|$, giving

[TABLE]

Also, as a consequence of complexity model (36) adopted for the approximate dictionaries $\tilde{\mathbf{A}}$, this criterion never causes a switching from $\tilde{\mathbf{A}}^{i}$ to $\tilde{\mathbf{A}}^{j}$ with $i<j\neq I$.

IV-C Heuristic look-ahead speed criterion

The preserved set size $|\mathcal{A}_{t}|$ used in (38) should ideally be the one associated to the dictionary after switching ($\mathbf{A}$ in this case), while in practice only the screening rate associated to the current approximation $\tilde{\mathbf{A}}^{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 $|\mathcal{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 $|\mathcal{A}|$ obtained with a conventional Dynamic Safe screening mechanism having access to $\mathbf{A}$. This is the quantity we are interested in estimating for the speed criterion (38). The y-axis corresponds to the number of preserved atoms obtained by the corresponding stable test for a given approximation $\tilde{\mathbf{A}}$. 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 $|\mathcal{A}_{t}|$ will reach after switching. This phenomenon is intensified for higher approximation errors.

As a heuristic, we propose to estimate $|\mathcal{A}_{t}|$ using the conventional test on the approximate atoms, i.e. $\mu(\tilde{\bm{\mathrm{a}}}_{j}|\mathcal{R})$.

[TABLE]

As shown in Fig. 2 (right graph), this correlates much better with the oracle value of $|\mathcal{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 (1), 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.

IV-D Summary and example

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

1. (I)

Decide when to stop using the current dictionary. 2. (II)

Choose the next dictionary.

The speed criterion (section IV-B) triggers a switching directly to the original dictionary $\tilde{\mathbf{A}}^{I}=\mathbf{A}$. We suppose the Relative Complexity associated to the i-th approximation $\text{RC}(\tilde{\mathbf{A}}^{i})$ (cf. (36)) to be known and stored in memory. The convergence criterion (section IV-A) switches to the next available approximation $\tilde{\mathbf{A}}^{i+1}$. Note that if $\mathbf{A}$ is already adopted, i.e. input $i=I$, then $\mathbf{A}$ is kept.

Example

Figure 3 shows an example of the proposed algorithm, given a list of approximations $\{\tilde{\mathbf{A}}^{i}\}_{i=0}^{I}$ with ${I=4}$, decreasing approximation error $\bm{\mathrm{\epsilon}}^{i}$ and increasing relative complexity $\text{RC}(\tilde{\mathbf{A}}^{i})\!=\!0.15({i\!+\!1})$, $0\!\leq\!i\!<\!I$. As usual, $\tilde{\mathbf{A}}^{I}\!=\!\mathbf{A}$.

Figure 3(a) shows the duality gap evolution over the iterations with $\Gamma=0.2$. As soon as the gap $G(\bm{\mathrm{x}}_{t},\bm{\mathrm{\theta}}^{\prime}_{t}|\tilde{\mathbf{A}^{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 $\mathbf{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. 3(b) in which the duality gap is plotted as a function of the execution time.

The theoretical complexity per iteration is shown in Fig. 3(d). 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 $\tilde{\mathbf{A}}^{i=3}$ was never used, since we switched directly from $\tilde{\mathbf{A}}^{i=2}$ to $\mathbf{A}$.

IV-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(\mathbf{\bm{\mathrm{x}}_{t},\mathbf{A},\bm{\mathrm{\alpha}}_{t}})$ (typically the product $\mathbf{A}^{T}\bm{\mathrm{\rho}}$) or that can be precalculated once for all ($\mathbf{A}^{T}\bm{\mathrm{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.

IV-E1 Screening cost

The most expensive computations associated to a dynamic screening rule – the ones potentially of the order of a matrix vector product $\mathcal{O}(NK)$ – are:

• •

Computation of a dual feasible point $\bm{\mathrm{\theta}}_{t}$ in (11) (resp. $\bm{\mathrm{\theta}}^{\prime}_{t}$(III-C)) requires the product $\mathbf{A}^{T}\bm{\mathrm{\rho}}$ (resp. $\tilde{\mathbf{A}}^{T}\tilde{\bm{\mathrm{\rho}}}$) reused from the optimization algorithm.

• •

Sphere test $\mu(\bm{\mathrm{a}}|\mathcal{R})$ (resp. $\mu(\mathcal{S}^{\prime}_{j}|\mathcal{R})$) requires the product $\bm{\mathrm{a}}^{T}_{j}\bm{\mathrm{c}}$ (resp. $\tilde{\bm{\mathrm{a}}}^{T}_{j}\bm{\mathrm{c}}$) for all preserved atoms $j\in\mathcal{A}$, which comes down to the matrix-vector product $\bm{\mathrm{A}}^{T}\bm{\mathrm{c}}$ (resp. $\tilde{\bm{\mathrm{A}}}^{T}\bm{\mathrm{c}}$). Practical sphere regions have been designed to limit this potential overhead. While in the Dynamic Safe sphere it reduces to the precalculated product $\mathbf{A}^{T}\bm{\mathrm{y}}$, in the GAP Safe sphere it reduces to $\mathbf{A}^{T}\bm{\mathrm{\rho}}$ (resp. $\tilde{\mathbf{A}}^{T}\tilde{\bm{\mathrm{\rho}}}$) calculated in the optimization iteration.

The other required calculations are detailed in Appendix -D. In short, the screening represents a rather low overhead $\mathcal{O}(N+|\mathcal{A}|)$ –even its stable version– compared to the optimization update: $\mathcal{O}(N|\mathcal{A}|)$ with screening or $\mathcal{O}(NK)$ without it, due to matrix-vector products.

IV-E2 Full iteration cost

Table IV shows the number of operations of a complete iteration in Algorithm 1 (ISTA update + screening), following [21] and adopting the complexity models in equations (35) and (36) for matrix vector multiplications. We denote $\text{flops}_{\mathbf{A}}(t)$ the cost of iteration $t$ with the conventional screening and $\text{flops}_{\tilde{\mathbf{A}}}(t)$ with the stable screening. As a benchmark, we use the complexity of an ISTA iteration without screening, denoted $\text{flops}_{N}(t)$.

To obtain the total complexity of the algorithm, simply add up all iteration costs calculated with the corresponding active set size $|\mathcal{A}_{t}|$ and sparsity of the solution estimate $\|\bm{\mathrm{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 $\ell_{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 $\bm{\mathrm{y}}$ is $N=2500$ and that of the coefficient vector $\bm{\mathrm{x}}$ (or, equivalently, the number of atoms in the dictionary) is $K=10000$. Unit-norm input data samples $\bm{\mathrm{y}}=\mathbf{A}\bm{\mathrm{x}}$ are generated from a sparse vector $\bm{\mathrm{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 $\bm{\mathrm{y}}$. The dictionary matrix $\mathbf{A}\in\mathbb{R}^{N\times 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.

V-A Specifying the Fast Structured dictionaries

The performance of the proposed algorithm is directly associated to the quality of the available approximations $\tilde{\mathbf{A}}^{i}$. A good approximation would have both a small application complexity (RC) and small approximation error ($\bm{\mathrm{\epsilon}}$). 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 [18] which can be written as a sum of Kronecker products $\tilde{\mathbf{A}}\!=\!\sum_{k=1}^{n_{\text{kron}}}\mathbf{B}_{k}\!\otimes\!\mathbf{C}_{k}$. Its reduced multiplication cost comes from the fact that the sub-matrices $\mathbf{B}_{k}$, $\mathbf{C}_{k}\in\mathbb{R}^{\sqrt{N}\!\times\!\sqrt{K}}$ are much smaller than $\mathbf{A}\in\mathbb{R}^{N\!\times\!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_{\text{kron}}$) in the sum. A higher $n_{\text{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_{\text{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_{\text{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. 444Some parallelization is introduced (in the summing terms) and matrix-matrix products are faster than the naive cubic complexity.

V-B Computational complexity and Time results

Let us denote $F_{N}$, $F_{\mathbf{A}}$ and $F_{\tilde{\mathbf{A}}}$ 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

[TABLE]

where $n_{\text{it}}$ is the number of iterations. When calculating $F_{\tilde{\mathbf{A}}}$ , the current approximation at iteration $t$ must be used in the expression of $\text{flops}_{\tilde{\mathbf{A}}}(t)$. Likewise, we denote $T_{N}$, $T_{\mathbf{A}}$ and $T_{\tilde{\mathbf{A}}}$ the measured running times.

As a main figure of merit, we evaluate the normalized number of flops ($F_{\mathbf{A}}/F_{N}$ and $F_{\tilde{\mathbf{A}}}/F_{N}$) and normalized running times ($T_{\mathbf{A}}/T_{N}$ and $T_{\tilde{\mathbf{A}}}/T_{N}$) 555Time results are obtained in an Intel® Core™ 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. .

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.

V-C Choosing the switching threshold

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

We empirically observed that the choice of $\Gamma$ is mostly dependent on the quality of the approximations $\tilde{\mathbf{A}}^{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 $\Gamma\in[10^{-2},0.8]$ for each one of the three scenarios. Each line corresponds to a different regularization level $\lambda$. For each of these lines, we are interested in the $\Gamma$ value that minimizes the running times. We can see that the low regularization configurations (dark blue curves) are more sensitive to choice of $\Gamma$, especially in the Hard scenario (left plot) in which a bad choice of $\Gamma$ can even lead to normalized running times greater than one –so, actually, a slowdown. This happens when a too small $\Gamma$ 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 $\Gamma$ 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 $\Gamma$ value can be chosen to obtain close to optimal execution times for any $\lambda/\lambda_{\max}$ in the tested range. Slightly different $\Gamma$ values can be chosen to better adapt to each simulation scenario: $\Gamma=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 $\Gamma$ 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 $\Gamma=0.5$, for instance.

V-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 $\Gamma$ 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 [25, 26]) 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 latter666The results reported in the right graphs are already normalized w.r.t the FISTA time results, which are typically smaller than those of ISTA.. 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.

V-E Speedup-error compromise

It is reasonable to expect the quality of the approximations $\tilde{\mathbf{A}}^{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 $\mathbf{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 speedup-error 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 non-negligible acceleration.

V-F Convergence tolerance

As previously highlighted, screening is known to work better on highly-regularized scenarios ($\lambda/\lambda_{\max}\approx 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 $\mathbf{A}$ are necessary. Thus, proportionally, the accelerated part of the algorithm (with fast dictionaries) is more significant.

V-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 [42, 43].

We consider an MEG gain matrix (forward operator) $\mathbf{A}\!\in\!\mathbb{R}^{204\!\times\!8193}$ obtained with the MNE software [44]. This operator can be efficiently approximated as a product of a few sparse matrices, as proposed in [16]. The speedups and corresponding approximation errors achieved by the so-called FA$\mu$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 $\bm{\mathrm{x}}\in\mathbb{R}^{8193}$ to be recovered from the input signals $\bm{\mathrm{y}}\!=\!\mathbf{A}\bm{\mathrm{x}}\in\mathbb{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 $\lambda/\lambda_{\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.777 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 task-specific performance measures are not discussed (see [43] for more details on the performance of $\ell_{1}$-minimization techniques for brain source localization).

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 $\ell_{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 $\ell_{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 [24]. Finally, an interesting perspective is to extend these zone-based screening tests to off-the-grid generalizations of $\ell_{1}$-regularized problems [45].

-A Proof of Equation (8)

Screening test for a $\ell_{p}$-ball, i.e. $\mathcal{R}=B_{p}(\bm{\mathrm{c}},R)$:

[TABLE]

using the definition of dual norm $\|\bm{\mathrm{a}}\|_{p^{*}}\!:=\!\sup_{\bm{\mathrm{u}}\in B_{p}(\bm{\mathrm{0}},1)}|\bm{\mathrm{a}}^{T}\bm{\mathrm{u}}|$.

-B Proof of Equation (18)

From Definition 3 and the conventional sphere test in eq. (8), the stable sphere test is given by

[TABLE]

Taking $\mathcal{S}=B_{q}(\tilde{\bm{\mathrm{a}}},\epsilon)$, the first term gives (similarly to (-A))

[TABLE]

and the last term gives

[TABLE]

with a constant $C=\sup_{\bm{\mathrm{u}}\in B_{q}(\bm{\mathrm{0}},1)}\|\bm{\mathrm{u}}\|_{p^{*}}=N^{\left(\frac{1}{p^{*}}-\frac{1}{q}\right)_{+}}$, which results from the Holder’s inequality. Substituting (-B) and (43) in (-B) gives the bound in equation (18).

-C Proof of Lemma 2

By definition we have $\bm{\mathrm{\theta}}_{F}:=\Theta^{\prime}(\bm{\mathrm{z}}|\tilde{\mathbf{A}},\bm{\mathrm{\epsilon}})=\alpha\bm{\mathrm{z}}$ where

[TABLE]

As a result for any $1\leq j\leq K$ we have

[TABLE]

This implies that the dual point $\bm{\mathrm{\theta}}_{F}$ is feasible, i.e. $\bm{\mathrm{\theta}}_{F}\in\Delta_{\mathbf{A}}$. The fact that $\bm{\mathrm{\theta}}_{F}\in\Delta_{\tilde{\mathbf{A}}}$ follows similarly (since $\epsilon_{j}\|\bm{\mathrm{z}}\|_{q^{*}}>0$).

-D Detailed complexities

We suppose that the $p^{*}$-norms of the atoms $\|\bm{\mathrm{a}}_{j}\|_{p^{*}}$ and $q$-norms of the atom approximation error $\|\bm{\mathrm{e}}_{j}\|_{q}$ are precalculated and stored and memory. The total screening complexity is summarized in Table IV and detailed in the following.

-D1 Feasible point $\boldsymbol{\theta}_{t}$

A total of at most $3N+|\mathcal{A}_{t}|$ operations for the point in eq. (11), distributed as follows:

$N$

dot product $\bm{\mathrm{y}}^{T}\bm{\mathrm{\rho}}_{t}$.

$N$

norm $\|\bm{\mathrm{\rho}}_{t}\|^{2}_{2}$ (if primal objective not calculated as a convergence criterion).

$N$

product $\alpha\bm{\mathrm{\rho}}_{t}$.

$|\mathcal{A}_{t}|$

comparisons for the infinity norm $\|\mathbf{A}^{T}\bm{\mathrm{\rho}}\|_{\infty}$.

At most $N\!+\!|\mathcal{A}_{t}|$ extra operations for the point in eq. (III-C):

$N$

norm $\|\tilde{\bm{\mathrm{\rho}}}_{t}\|_{q^{*}}$ (only if $q^{*}\neq 2$).

$|\mathcal{A}_{t}|$

products $\epsilon_{j}\|\tilde{\bm{\mathrm{\rho}}}_{t}\|_{q^{*}}$ on the calculation of $\alpha^{\prime}_{t}$.

-D2 Safe region $\mathcal{R}$

At most $N$ operations for computing $\|\bm{\mathrm{\theta}}_{t}-\bm{\mathrm{y}}/\lambda\|_{2}$ (if the dual objective $D(\bm{\mathrm{\theta}}_{t})$ is not already calculated as a convergence criterion) for the radius of either Dynamic Safe or GAP Safe sphere.

-D3 Screening test

Other operations (aside from the mentioned $\mathbf{A}^{T}\bm{\mathrm{c}}$) sum up to $|\mathcal{A}_{t}|$ operations for the product $R\|\bm{\mathrm{a}}_{j}\|_{p^{*}}~{}\forall~{}j\in\mathcal{A}_{t}$. The proposed stable screening tests require at most an extra $N+|\mathcal{A}_{t}|$ operations:

$N$

norm $\|\bm{\mathrm{c}}\|_{q^{*}}$ (if $q^{*}\neq 2$ and $\bm{\mathrm{c}}$ varies with $t$).

$|\mathcal{A}_{t}|$

products $\epsilon_{j}\|\bm{\mathrm{c}}\|_{q^{*}}~{}\forall~{}j\in\mathcal{A}_{t}$ (if $\bm{\mathrm{c}}$ varies with $t$).