Clément Dorffer

Cédric Angélique Drémeau

Efficient atom selection

come

Efficient atom selection strategy for iterative sparse approximations

Clément Dorffer1 * , Angélique Drémeau 1 and Cédric Herzet2

1 Lab-STICC UMR 6285, CNRS, ENSTA Bretagne, Brest, F-29200, France. 2 INRIA Centre Rennes-Bretagne Atlantique and Lab-STICC UMR 6285, CNRS, IMT-Atlantique, Rennes, F-35000, France.

Abstract-We propose a low-computational strategy for the efficient implementation of the "atom selection step" in sparse representation algorithms. The proposed procedure is based on simple tests enabling to identify subsets of atoms which cannot be selected. Our procedure applies on both discrete or continuous dictionaries. Experiments performed on DOA and Gaussian deconvolution problems show the computational gain induced by the proposed approach.

1 Problem statement Sparsely approximating a signal vector y ∈ R m in a dictionary A consists of finding k m coefficients x i and atoms a i ∈ A such that y ≈ k i=1 a i x i . The dictionary A can be either discrete, i.e., composed of a finite number of elements, or "continuous", i.e., having an infinite uncountable number of atoms.

Sparse approximations have proven to be relevant in many application domains and a great number of procedures to find "good" sparse approximations have been proposed in the literature: convex relaxation [START_REF] Frank | An algorithm for quadratic programming[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF][START_REF] Rao | Forward-backward greedy algorithms for atomic norm regularization[END_REF][START_REF] Ekanadham | Recovery of sparse translation-invariant signals with continuous basis pursuit[END_REF][START_REF] Tang | Compressed sensing off the grid[END_REF][START_REF] Xenaki | Grid-free compressive beamforming[END_REF][START_REF] Catala | A lowrank approach to off-the-grid sparse deconvolution[END_REF][START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF], greedy algorithms [START_REF] Stéphane | Matching pursuits with time-frequency dictionaries[END_REF][START_REF] Chandra Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF][START_REF] Tropp | Algorithms for simultaneous sparse approximation. part i: Greedy pursuit[END_REF][START_REF] Needell | Cosamp: Iterative signal recovery from incomplete and inaccurate samples[END_REF][START_REF] Karin C Knudson | Inferring sparse representations of continuous signals with continuous orthogonal matching pursuit[END_REF], Bayesian approaches [START_REF] Févotte | Blind separation of sparse sources using jeffrey's inverse prior and the em algorithm[END_REF][START_REF] Schniter | Fast bayesian matching pursuit[END_REF][START_REF] Soussen | From bernoulli-gaussian deconvolution to sparse signal restoration[END_REF][START_REF] Drémeau | Soft bayesian pursuit algorithm for sparse representations[END_REF], etc. Many popular instances of these procedures rely on the same "atom selection" step, i.e., the new atom added to the support at each iteration, says a select , verifies

a select ∈ arg max a∈A r, a , (1)
where •, • denotes the inner product 1 and r is the current "residual", i.e., the original signal from which the contributions from the previously selected atoms have been removed. The Frank-Wolfe algorithm [START_REF] Frank | An algorithm for quadratic programming[END_REF], the matching pursuit [START_REF] Stéphane | Matching pursuits with time-frequency dictionaries[END_REF] or the orthogonal matching pursuit [START_REF] Chandra Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] procedures are popular instances of algorithms using [START_REF] Frank | An algorithm for quadratic programming[END_REF]. It is worth noting that (1) constitutes the core of most sparse-approximation algorithms in the context of continuous dictionaries [START_REF] Ekanadham | Recovery of sparse translation-invariant signals with continuous basis pursuit[END_REF][START_REF] Karin C Knudson | Inferring sparse representations of continuous signals with continuous orthogonal matching pursuit[END_REF] where some standard "matrix-vector" operations available in the discrete setting are no longer possible. A brute-force evaluation 2 of (1) may become resource consuming when the number of atoms in A is large since it requires the exploration of the whole dictionary to find the atom the most correlated with r. In this work, we propose a strategy to alleviate the complexity of (1). Our procedure is inspired from work [START_REF] Herzet | Joint screening tests for lasso[END_REF]: it consists of performing simple tests allowing to identify group of atoms not attaining the maximum value of (1). Interestingly, the proposed approach provides a rigorous framework to recently-proposed procedures based on some approximations of continuous dictionaries [START_REF] Ekanadham | Recovery of sparse translation-invariant signals with continuous basis pursuit[END_REF][START_REF] Karin C Knudson | Inferring sparse representations of continuous signals with continuous orthogonal matching pursuit[END_REF]. In the rest of this abstract, we do not elaborate on the connections with [START_REF] Ekanadham | Recovery of sparse translation-invariant signals with continuous basis pursuit[END_REF][START_REF] Karin C Knudson | Inferring sparse representations of continuous signals with continuous orthogonal matching pursuit[END_REF] (details will be provided during the conference) but rather focus on the description of the proposed methodology.

Proposed strategy

Our proposed selection strategy is based on the following observations. If Ā ⊆ A,

max a∈ Ā r, a ≤ max a∈A r, a . Hence, letting τ max a∈ Ā r, a , we have ∀a ∈ A: r, a < τ ⇒ a / ∈ arg max ã∈A r, ã . (2)
In other words, if a ∈ A is an atom which satisfies the inequality in the left-hand side of (2), then this atom is surely not the one to be selected by [START_REF] Frank | An algorithm for quadratic programming[END_REF]. Elaborating on this observation, we further have:

max a∈R r, a < τ ⇒ ∀a ∈ A ∩ R : a / ∈ arg max ã∈A r, ã , (3)
where R is some arbitrary subset of R m . In the sequel we will refer to R as "region". The operational meaning of (3) is as follows: if the inequality in the left-hand side is satisfied, one is ensured that no atom in A ∩ R will attain the maximum of r, a . The entire set A ∩ R can thus be ignored, enabling us to reduce the number of candidate atoms to be tested in the selection (1). Implication (3) constitutes the basis of our complexity reduction method. More specifically, we consider (3) with some particular choices of region R. These choices are motivated by the following requirements: i) R should lead to any easy evaluation of max a∈R r, a ; ii) R should approximate as tightly as possible some part of A since larger regions typically lead to inequalities more difficult to satisfy.

In our contribution, we show that the first requirement is satisfied for some particular geometries of region R. In particular, we consider "sphere", "dome" and "slice" geometries. Sphere and dome regions can be formally expressed as

B t, {a : a -t 2 ≤ } (sphere) D t, {a : t, a ≥ , a 2 = 1} (dome),
whereas the mathematical characterization of the slice regions is more involved and not detailed in this abstract. We show that for these choices of regions, max a∈R r, a admits a simple

Algorithm 1 Efficient selection strategy

Input: residual r, subset of atoms Ā, set of regions {R l } L l=1 Init: :

A removed = ∅ Evaluate τ = max a∈ Ā r, a for all 1 ≤ l ≤ L do if max a∈R l r, a < τ then A removed = A removed ∪ (A ∩ R l ) end if end for Find a select ∈ arg max a∈A\A removed r, a
analytical expression. In particular, evaluating max a∈R r, a for sphere and dome regions basically requires the computation of one single inner product t, r . We address the second requirement by proposing a methodology to automatically adapt the size of R (via a tuning of ) in order to satisfy the inequality in the left-hand side of (3). We do not detail this procedure here but mention that the evaluation of the "optimal" value of has a negligible complexity.

We thus propose the following strategy (summarized in Algorithm 1) to speed up the computation of (1). We select a set Ā ⊂ A and L regions {R l } L l=1 , and apply test (3) for each region. Each test allows to identify a set 3 of atoms which do not attain the maximum of r, a . We evaluate (1) by working on a reduced dictionary:

a select ∈ arg max a∈A\A removed r, a , (4)
where A removed denotes the set of atoms which have been removed by the tests (3). For discrete dictionaries, the computational complexity of the proposed method (for sphere and dome regions) evolves as

O( card( Ā)m evaluation of τ + Lm evaluation of (3)
+ card(A\A removed )m evaluation of ( 4)

).

This has to be compared to the complexity required by a brute evaluation of (1), i.e., O(card(A)m). In the next section, we propose to compare the efficiency and the computational gain allowed by the proposed methodology.

Experiments

We propose to challenge the proposed selection procedure and the classical exhaustive search on two different problems: a direction-of-arrival (DOA) estimation problem using a discrete dictionary and a Gaussian deconvolution problem with a continuous dictionary. Within the DOA estimation framework, we examine the number of scalar products required to achieve the selection step (1), with and without the proposed method. This provides a quantitative assessment of the computational gain induced by our method. The Gaussian deconvolution problem acts as a proof of concept, illustrating the interest of our method in continuous dictionaries.

In the DOA estimation problem, we consider a dictionary composed of n = 1000 normalized steering vectors of size m = 100, each corresponding to a different angle of arrival in [-π/2, π/2]. In the Gaussian deconvolution problem, 3 Which may be empty in some cases. In Fig. 1(a), we considered the DOA estimation problem. We applied 5 iterations of orthogonal matching pursuit on y (which requires to solve (1) at each iteration). Each column of the figure represents the (cumulated) number of inner products which have been computed with an exhaustive search of a select in the whole dictionary (blue), and with the proposed method described in Algorithm 1 (red). We set L = 100 and define the centers of the (dome) regions {t l } L l=1 by a regular subsampling of the dictionary. We let Ā = {t l } L l=1 . We see that the proposed method allows for a gain of complexity of one order of magnitude with respect to a brute-force approach.

In Fig. 1(b) we consider the Gaussian deconvolution problem and show the set of removed atoms, A ∩ R l , for each (sphere) region {R l } L=100 l=1 . We choose the centers of the regions {t l } L l=1 by a regular subsampling of µ ∈ [0, 100]. The size of each region (value of ) is automatically tuned to verify (if possible) the inequality in the left-hand side (3) as discussed in Section 2. We set Ā = {t l } L=100 l=1 . The figure shows the interval of value of µ which can be removed from the dictionary for each test region. We see in Fig. 1(b) that by the end of the proposed elimination procedure, the search space is reduced to a small interval, around µ ≈ 90.

Figure 1 :

 1 Figure 1: (a) : "discrete" DOA estimation problem -number of scalar products evaluated for the atom selection. (b) : "continuous" Gaussian deconvolution problem -interval of atoms which can be ignored in problem (1).

In case of possibly negative coefficients in x, it is common to use the absolute value of the inner product. However, one can also deal with negative coefficients using eq. (1) with a doubled size dictionary containing atoms a i and their negatives -a i .

In the discrete setting, the standard approach consists of evaluating r, a for all a ∈ A, leading to a complexity scaling as O(card(A)m). In the continuous setting, a classic approach consists of running gradient ascent algorithms initialized on a fine discretization of the dictionary.

The authors thank the DGA/MRIS, the ONR (N62909-17-1-2007) and the ANR (ANR-15-CE23-0021) for their financial support.