Charles Soussen
email: charles.soussen@univ-lorraine.fr

Jérome Idier
email: jerome.idier@irccyn.ec-nantes.fr.

Junbo Duan
email: junbo.duan@mail.xjtu.edu.cn

David Brie
email: david.brie@univ-lorraine.fr.

Jérôme Idier

J Idier

Homotopy based algorithms for ℓ 0 -regularized least-squares

HAL is

Index Terms

Sparse signal estimation; ℓ 0 -regularized least-squares; ℓ 0 -homotopy; ℓ 1 -homotopy; stepwise algorithms; orthogonal least squares; model order selection.

I. INTRODUCTION

Sparse approximation from noisy data is traditionally addressed as the constrained least-square problems

min x y -Ax 2 2 subject to x 0 ≤ k (1)
or

min x x 0 subject to y -Ax 2 2 ≤ ε (2)
where x 0 is the ℓ 0 -"norm" counting the number of nonzero entries in x, and the quadratic fidelity-todata term y -Ax 2 2 measures the quality of approximation. Formulation [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF] is well adapted when one has a knowledge of the maximum number k of atoms to be selected in the dictionary A. On the contrary, the choice of ( 2) is more appropriate when k is unknown but one has a knowledge of the variance of the observation noise. The value of ε may then be chosen relative to the noise variance. Since both [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF] and ( 2) are subset selection problems, they are discrete optimization problems. They are known to be NP-hard except for specific cases [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF].

When no knowledge is available on either k or ε, the unconstrained formulation

min x {J (x; λ) = y -Ax 2 2 + λ x 0 } (3)
is worth being considered, where λ expresses the trade-off between the quality of approximation and the sparsity level [START_REF] Nikolova | Description of the minimizers of least squares regularized with ℓ0 norm. Uniqueness of the global minimizer[END_REF]. In a Bayesian viewpoint, (3) can be seen as a (limit) maximum a posteriori formulation where y-Ax 2 2 and the penalty x 0 are respectively related to a Gaussian noise distribution and a prior distribution for sparse signals (a limit Bernoulli-Gaussian distribution with infinite Gaussian variance) [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF].

A. Classification of methods 1) ℓ 0 -constrained least-squares: The discrete algorithms dedicated to problems (1)-( 2) can be categorized into two classes. First, the forward greedy algorithms explore subsets of increasing cardinalities starting from the empty set. At each iteration, a new atom is appended to the current subset, therefore gradually improving the quality of approximation [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]. Greedy algorithms include, by increasing order of complexity: Matching Pursuit (MP) [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF], Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF], and Orthogonal Least [START_REF] Friedman | Multivariate adaptive regression splines[END_REF] DRAFT Squares (OLS) [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF], also referred to as forward selection in statistical regression [START_REF] Miller | Subset selection in regression[END_REF] and known as Order

Recursive Matching Pursuit (ORMP) [START_REF] Cotter | Forward sequential algorithms for best basis selection[END_REF] and Optimized Orthogonal Matching Pursuit (OOMP) [START_REF] Rebollo-Neira | Optimized orthogonal matching pursuit approach[END_REF]. The second category are thresholding algorithms, where each iteration delivers a subset of same cardinality k. Popular thresholding algorithms include Iterative Hard Thresholding [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF], Subspace Pursuit [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF] and CoSaMP [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF].

Among these two categories, greedy algorithms are well-adapted to the resolution of ( 1) and (2) for variable sparsity levels. Indeed, they yield a series of subsets for consecutive k (i.e., for decreasing approximation errors ε) since at each iteration, the current subset is increased by one element.

2) ℓ 0 -penalized least-squares: In [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF], we evidenced that the minimization of J (x; λ) using a descent algorithm leads to bidirectional extensions of forward (orthogonal) greedy algorithms. To be more specific, consider a candidate subset S corresponding to the support of x. Including a new element into S yields a decrease of the square error, defined as the minimum of y -Ax 2 2 for x supported by S. On the other hand, the penalty term λ x 0 is increased by λ. Overall, the cost function J (x; λ) decreases as soon as the square error variation exceeds λ. Similarly, a decrease of J (x; λ) occurs when an element is removed from S provided that the squared error increment is lower than λ. Because both inclusion and removal operations can induce a decrease of J , the formulation (3) allows one to design descent schemes allowing a "forward-backward" search strategy, where each iteration either selects a new atom (forward selection) or de-selects an atom that was previously selected (backward elimination). The Bayesian OMP [START_REF] Herzet | Bayesian pursuit algorithms[END_REF] and Single Best Replacement (SBR) [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] algorithms have been proposed in this spirit. They are extensions of OMP and OLS, respectively. Their advantage over forward greedy algorithms is that an early wrong atom selection may be later cancelled. Forward-backward algorithms include the so-called stepwise regression algorithms which are OLS extensions [START_REF] Miller | Subset selection in regression[END_REF], [START_REF] Efroymson | Multiple regression analysis[END_REF], [START_REF] Berk | Forward and backward stepping in variable selection[END_REF], and OMP based algorithms of lower complexity [START_REF] Herzet | Bayesian pursuit algorithms[END_REF], [START_REF] Zhang | Adaptive forward-backward greedy algorithm for learning sparse representations[END_REF].

3) Connection with the continuous relaxation of the ℓ 0 norm: The algorithms described so far are discrete search strategies dedicated to ℓ 0 -regularized least-squares. A classical alternative consists in relaxing the ℓ 0 -norm by a continuous function that is nondifferentiable at 0, and optimizing the resulting cost function. See, e.g., [START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF], [START_REF] Zibulevsky | ℓ1ℓ2 optimization in signal and image processing[END_REF] and [START_REF] Candès | Enhancing sparsity by reweighted ℓ1 minimization[END_REF]- [START_REF] Selesnick | Sparse signal estimation by maximally sparse convex optimization[END_REF] for convex (ℓ 1 ) and nonconvex relaxation, respectively.

The convex problem min x y -Ax 2 2 s.t. x 1 ≤ t is referred to as both Basis Pursuit Denoising (BPDN) and the LASSO. It is noticeable that BPDN leads to stepwise algorithms [START_REF] Figueiredo | Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems[END_REF], [START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF] including the popular ℓ 1 -homotopy [START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF]- [START_REF] Efron | Least angle regression[END_REF], a forward-backward greedy search whose complexity is close to that of OMP. ℓ 1 -homotopy is closely connected to the Least Angle Regression (LARS), a simpler forward strategy allowing only atom selections. It is referred to as "LARS with the LASSO modification" in [START_REF] Efron | Least angle regression[END_REF]. [START_REF] Friedman | Multivariate adaptive regression splines[END_REF] DRAFT

λ S ⋆ 2 E(S ⋆ 0 ) E(S ⋆ 1 ) E(S) + λ|S| λ ⋆ 2 λ ⋆ 1 E(S ⋆ 2 ) S ⋆ 1 S ⋆ 0 = ∅ λ ⋆ 3 = 0
Fig. 1. Representation of lines λ → E (S) + λ|S| for various subsets S. The ℓ0-curve, in plain line, is the minimal curve λ → minS{E (S) + λ|S|}. It is continuous, concave, and piecewise affine with a finite number of pieces. The ℓ0-penalized regularization path is composed of the supports (here, S ⋆ 0 , S ⋆ 1 , S ⋆ 2 ) that are optimal for some λ-values. For instance, S ⋆ 1 is optimal for λ ∈ [λ ⋆ 2 , λ ⋆ 1 ]. These supports S ⋆ induce global minimizers of J (x; λ), defined as the least-square solutions xS⋆. For instance, x S ⋆ 1 is a global minimizer of J (x; λ) with respect to

x whenever λ ∈ [λ ⋆ 2 , λ ⋆ 1 ].
Importantly, ℓ 1 -homotopy solves the BPDN for a continuum of values of t.

B. Main idea

Our approach is dedicated to ℓ 0 -penalized least-squares. It is based on the following geometrical interpretation.

First, for any subset S, we can define a linear function λ → E(S) + λ|S|, where

E(S) = y -Ax 2 2
is the corresponding least-square error and |S| stands for the cardinality of S. For each subset S, this function yields a line in the 2D domain (λ, J ), as shown on Fig. 1.

Second, the set of solutions to (3) is piecewise constant with respect to λ (see Appendix A for a proof).

Geometrically, this result can be easily understood by noticing that the minimum of J (x; λ) with respect to x is obtained for all λ-values by considering the concave envelope of the set of lines λ → E(S) + λ|S| for all subsets S. The resulting piecewise affine curve is referred to as the ℓ 0 -curve (see Fig. 1). Its edges are related to the supports of the sparse solutions for all λ, and its vertices yield the breakpoints λ ⋆ i around which the set of optimal solutions arg min x J (x; λ) is changing.

We take advantage of this interpretation to propose two suboptimal greedy algorithms that address [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] for a continuum of λ-values. Continuation Single Best Replacement (CSBR) repeatedly minimizes J (x; λ)

with respect to x for decreasing λ's. ℓ 0 Regularization Path Descent (ℓ 0 -PD) is a more complex algorithm maintaining a list of subsets so as to improve (decrease) the current approximation of the ℓ 0 curve.

C. Related works 1) Bi-objective optimization:

The formulations (1), ( 2) and ( 3) can be interpreted as the same biobjective problem because they all intend to minimize both the approximation error y -Ax 2 2 and the sparsity measure x 0 . Although x is continuous, the bi-objective optimization problem should rather be considered as a discrete one where both objectives reread E(S) and |S|. Indeed, the continuous solutions deduce from the discrete solutions, x reading as a least-square minimizer among all vectors supported by S. Fig. 2 is a classical bi-objective representation where each axis is related to a single objective [START_REF] Das | A closer look at drawbacks of minimizing weighted sums of objectives for Pareto set generation in multicriteria optimization problems[END_REF], namely |S| and E(S). In bi-objective optimization, a point S is called Pareto optimal when no other point S ′ can decrease both objectives [START_REF] Marler | Survey of multi-objective optimization methods for engineering[END_REF]. In the present context, |S| takes integer values, thus the Pareto solutions are the minimizers of E(S) subject to |S| ≤ k for consecutive values of k. Equivalently, they minimize |S| subject to E(S) ≤ ε for some ε. They are usually classified as supported or non-supported.

The former lay on the convex envelope of the Pareto frontier (the bullet points in Fig. 2) whereas the latter lay in the nonconvex areas (the square point). It is well known that a supported solution can be reached when minimizing the weighted sum of both objectives, i.e., when minimizing E(S) + λ|S| with respect to S for some weight λ. On the contrary, the non-supported solutions cannot [START_REF] Marler | Survey of multi-objective optimization methods for engineering[END_REF]. Choosing between the weighting sum method and a more complex method is a nontrivial question. The answer depends on the problem at-hand and specifically, on the size of the nonconvex areas in the Pareto frontier. 2) ℓ 1 and ℓ 0 -homotopy seen as a weighted sum method: It is important to notice that for convex objectives, the Pareto solutions are all supported. Consider the BPDN; because y -Ax 2 2 and x 1 are convex functions of x, the set of minimizers of y -Ax 2 2 + λ x 1 for all λ coincides with the set of minimizers of y -Ax 2 2 s.t. x 1 ≤ t for all t [START_REF] Van Den | Probing the Pareto frontier for basis pursuit solutions[END_REF]. Both sets are referred to as the (unique) "ℓ 1 -regularization path". The situation is different with ℓ 0 -regularization. Now, the weighted sum formulation (3) may not yield the same solutions as the constrained formulations (1) and ( 2) because the ℓ 0 -norm is nonconvex [START_REF] Nikolova | Description of the minimizers of least squares regularized with ℓ0 norm. Uniqueness of the global minimizer[END_REF]. This will lead us to define two ℓ 0 -regularization paths, namely the "ℓ 0 -penalized path" and the "ℓ 0 -constrained path" (Section II).

On the algorithmic side, the ℓ 0 problems are acknowledged to be difficult. Many authors actually discourage the direct optimization of J because there are a very large number of local minimizers [START_REF] Candès | Enhancing sparsity by reweighted ℓ1 minimization[END_REF], [START_REF] Wipf | Iterative reweighted ℓ1 and ℓ2 methods for finding sparse solutions[END_REF]. In [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF], however, we showed that forward-backward extensions of OLS are able to escape from some local minimizers of J (x; λ) for a given λ. This motivates us to propose efficient OLS-based strategies for minimizing J for variable λ-values.

3) Positioning with respect to other stepwise algorithms: In statistical regression, the word "stepwise" originally refers to Efroymson's algorithm [START_REF] Efroymson | Multiple regression analysis[END_REF], proposed in 1960 as an empirical extension of forward selection (i.e., OLS). Other stepwise algorithms were proposed in the 1980's [8, Chapter 3] among which Berk's and Broersen's algorithms [START_REF] Berk | Forward and backward stepping in variable selection[END_REF], [START_REF] Broersen | Subset regression with stepwise directed search[END_REF]. All these algorithms perform a single replacement per iteration, i.e., a forward selection or a backward elimination. They were originally applied to overdetermined problems in which the number of columns of A is lower than the number of rows. Recent stepwise algorithms were designed as either OMP [START_REF] Herzet | Bayesian pursuit algorithms[END_REF], [START_REF] Zhang | Adaptive forward-backward greedy algorithm for learning sparse representations[END_REF] or OLS extensions [START_REF] Haugland | A bidirectional greedy heuristic for the subspace selection problem[END_REF], [START_REF] Chatterjee | Projection-based and look-ahead strategies for atom selection[END_REF]. They all aim to find subsets of cardinality k yielding a low approximation error E(S) for all k. Although our algorithms share the same objective, they are inspired by (i) the ℓ 1 -homotopy algorithm; and (ii) the structural properties of the ℓ 0 -regularization paths. To the best of our knowledge, the idea of reconstructing an ℓ 0 -regularization path using ℓ 0 -homotopy procedures is novel.

CSBR and ℓ 0 -PD both read as descent algorithms in different senses: CSBR, first sketched in [START_REF] Duan | A continuation approach to estimate a solution path of mixed L2-L0 minimization problems[END_REF], repeatedly minimizes J (x; λ) for decreasing λ's. On the contrary, ℓ 0 -PD minimizes J (x; λ) for any λvalue simultaneously by maintaining a list of candidate subsets. The idea of maintaining a list of support candidates was recently developed within the framework of forward selection [START_REF] Kwon | Multipath matching pursuit[END_REF], [START_REF] Maymon | The Viterbi algorithm for subset selection[END_REF]. Our approach is different, because a family of optimization problems are being addressed together. In contrast, the supports in the list are all candidate solutions to solve the same problem in [START_REF] Kwon | Multipath matching pursuit[END_REF], [START_REF] Maymon | The Viterbi algorithm for subset selection[END_REF].

4) Positioning with respect to continuation algorithms:

The principle of continuation is to handle a difficult problem by solving a sequence of simpler problems with warm start initialization, and gradually [START_REF] Friedman | Multivariate adaptive regression splines[END_REF] DRAFT tuning some continuous hyperparameter [START_REF] Wasserstrom | Numerical solutions by the continuation method[END_REF]. In sparse approximation, the word continuation is used in two opposite contexts.

First, the BDPN problem involving the ℓ 1 -norm. BPDN is solved for decreasing hyperparameter values using the solution for each value as a warm starting point for the next value [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]. ℓ 1 -homotopy [START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF], [START_REF] Efron | Least angle regression[END_REF], [START_REF] Malioutov | Homotopy continuation for sparse signal representation[END_REF] exploits that the ℓ 1 regularization path is piecewise affine and tracks the breakpoints between consecutive affine pieces. CSBR is designed in a similar spirit and can be interpreted as an "ℓ 0 -homotopy" procedure (although the ℓ 0 minimization steps are solved in a sub-optimal way) working for decreasing λ-values.

Second, the continuous approximation of the (discrete) ℓ 0 pseudo-norm [START_REF] Trzasko | Highly undersampled magnetic resonance image reconstruction via homotopic ℓ0minimization[END_REF] using a Graduated Non Convexity (GNC) approach [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF]: a series of continuous concave metrics is considered leading to the resolution of continuous optimization problems with warm start initialization. Although the full reconstruction of the ℓ 0 -regularization path has been rarely addressed, it is noticeable that a GNC-like approach, called SparseNet, aims to gradually update some estimation of the regularization path induced by increasingly non-convex sparsity measures [START_REF] Mazumder | SparseNet: Coordinate descent with nonconvex penalties[END_REF]. This strategy relies on the choice of a grid of λ-values. Because the influence of the grid is critical [START_REF] Van Den | Probing the Pareto frontier for basis pursuit solutions[END_REF], it is useful to adapt the grid while the nonconvex measure is modified [START_REF] Mazumder | SparseNet: Coordinate descent with nonconvex penalties[END_REF]. On the contrary, our approach does not rely on a grid definition. The λ-values are rather adaptively computed similar to the ℓ 1 -homotopy principle [START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF], [START_REF] Efron | Least angle regression[END_REF].

The paper is organized as follows. In Section II, we define the ℓ 0 -regularization paths and establish their main properties. The CSBR and ℓ 0 -PD algorithms are respectively proposed in Sections III and IV.

In Section V, both algorithms are analyzed and compared with the state-of-art algorithms based on nonconvex penalties for difficult inverse problems. Additionally, we investigate the automatic choice of the cardinality k using classical order selection rules.

II. ℓ 0 -REGULARIZATION PATHS

A. Definitions, terminology and working assumptions

Let m×n denote the size of the dictionary A (usually, m ≤ n in sparse approximation). The observation signal y and the weight vector x are of size m×1 and n×1, respectively. We assume that any min(m, n) columns of A are linearly independent so that for any subset S ⊂ {1, . . . , n}, the submatrix of A gathering the columns indexed by S is full rank, and the least-square error E(S) can be numerically computed.

This assumption is however not necessary for the theoretical results provided hereafter.

We denote by |S| the cardinality of a subset S. We use the alternative notations "S +{i}" and "S -{i}"

for the forward selection S ∪ {i} and backward elimination S \ {i}. We can then introduce the generic March 18, 2015 DRAFT notation S ± {i} for single replacements: S ± {i} stands for S + {i} if i / ∈ S, and S -{i} if i ∈ S. We will frequently resort to the geometrical interpretation of Fig. 1. With a slight abuse of terminology, the line λ → E(S) + λ|S| will be simply referred to as "the line S".

Hereafter, we start by defining the ℓ 0 -regularized paths as the set of supports of the solutions to problems (1), ( 2) and (3) for varying hyperparameters. As seen in Section I, the solutions may differ whether the ℓ 0 -regularization takes the form of a bound constraint or a penalty. This will lead us to distinguish the "ℓ 0 -constrained path" and the "ℓ 0 -penalized path". We will keep the generic terminology "ℓ 0 -regularization paths" for statements that apply to both. The solutions delivered by our greedy algorithms will be referred to as the "approximate ℓ 0 -penalized path" since they are suboptimal algorithms.

B. Definition and properties of the ℓ 0 -regularized paths

The continuous problems (1), ( 2) and ( 3) can be converted as the discrete problems:

min S E(S) subject to |S| ≤ k, (4) min
S |S| subject to E(S) ≤ ε, (5) min
S Ĵ (S; λ) E(S) + λ|S| , (6)
where S stands for the support of x. The optimal solutions x to problems (1), ( 2) and (3) can indeed be simply deduced from those of (4), ( 5) and [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF], respectively, x reading as the least-square minimizers among all vectors supported by S. In the following, the formulation (5) will be omitted because it leads to the same ℓ 0 -regularization path as formulation (4) [START_REF] Nikolova | Description of the minimizers of least squares regularized with ℓ0 norm. Uniqueness of the global minimizer[END_REF].

Let us first define the set of solutions to (4) and ( 6) and the ℓ 0 -curve, related to the minimum value in [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] for all λ > 0.

Definition 1 For k ≤ min(m, n), let S ⋆ C (k) be the set of minimizers of the constrained problem (4). For λ > 0, let S ⋆ P (λ) be the set of minimizers of the penalized problem [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF]. Additionally, we define the ℓ 0 -curve as the function λ → min S { Ĵ (S; λ)}. It is the concave envelope of a finite number of linear functions. Thus, it is concave and piecewise affine. Let λ

⋆ I+1 0 < λ ⋆ I < . . . < λ ⋆ 1 < λ ⋆ 0 +∞ delimit
the affine intervals (I + 1 contiguous intervals; see Fig. 1 in the case where I = 2).

Each set S ⋆ C (k) or S ⋆ P (λ) can be thought of as a single support (e.g., S ⋆ C (k) is reduced to the support S ⋆a in the example of Fig. 2). They are defined as sets of supports because the minimizers of ( 4) and [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] might not be always unique. Let us now provide a key property of the set S ⋆ P (λ).

March 18, 2015 DRAFT

Theorem 1 S ⋆ P (λ) is a piecewise constant function of λ, being constant on each interval λ ∈ (λ ⋆ i+1 , λ ⋆ i ). Proof: See Appendix A.
This property allows us to define the ℓ 0 -regularization paths in a simple way.

Definition 2 The ℓ 0 -constrained path is the set (of sets) S ⋆ C = {S ⋆ C (k), k = 0, . . . , min(m, n)}. The ℓ 0 -penalized path is defined as S ⋆ P = {S ⋆ P (λ), λ > 0}. According to Theorem 1, S ⋆ P is composed of (I + 1) distinct sets S ⋆ P (λ), one for each interval λ ∈ (λ ⋆ i+1 , λ ⋆ i ).
S ⋆ C gathers the solutions to (4) for all k. As illustrated on Fig. 2, the elements of S ⋆ C are the Pareto solutions whereas the elements of S ⋆ P correspond to the convex envelope of the Pareto frontier. Therefore, both ℓ 0 -regularization paths may not coincide [START_REF] Nikolova | Description of the minimizers of least squares regularized with ℓ0 norm. Uniqueness of the global minimizer[END_REF], [START_REF] Das | A closer look at drawbacks of minimizing weighted sums of objectives for Pareto set generation in multicriteria optimization problems[END_REF]. As stated in Theorem 2, S ⋆ P ⊂ S ⋆ C , but the reverse inclusion is not guaranteed.

Theorem 2 S ⋆ P ⊂ S ⋆ C . Moreover, for any λ / ∈ {λ ⋆ I , . . . , λ ⋆ 0 }, there exists k such that S ⋆ P (λ) = S ⋆ C (k).
Proof: See Appendix A.

C. Approximate ℓ 0 -penalized regularization path

Let us introduce notations for the approximate ℓ 0 -penalized path delivered by our heuristic search algorithms. Throughout the paper, the ⋆ notation is reserved for optimal solutions (e.g., S ⋆ P ). It is removed when dealing with numerical solutions. The outputs of our algorithms will be composed of a list λ = {λ 1 , . . . , λ J+1 } of decreasing λ-values, and a list S = {S 0 , . . . , S J } of candidate supports, with S 0 = ∅. S j is a suboptimal solution to (6) for λ ∈ (λ j+1 , λ j ). In the first interval λ > λ 1 , the solution is S 0 = ∅. The reader shall keep in mind that each output S j induces a suboptimal solution x j to (3) for λ ∈ (λ j+1 , λ j ). This vector is the least-square solution supported by S j . It can be computed using the pseudo-inverse of the subdictionary indexed by the set of atoms in S j .

Geometrically, each support S j yields a line segment. Appending these segments yields an approximate ℓ 0 -curve covering the domain (λ J+1 , +∞), as illustrated on Fig. 3.

III. GREEDY CONTINUATION ALGORITHM (CSBR)

Our starting point is the Single Best Replacement algorithm [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] dedicated to the minimization of J (x; λ) with respect to x, or equivalently to Ĵ (S; λ) = E(S) + λ|S| with respect to S. We first describe SBR for a given λ. Then, the CSBR extension is presented for decreasing and adaptive λ's.

3) Update S ← S best .

SBR terminates when Ĵ (S best ; λ) ≥ Ĵ (S; λ), i.e., when no single replacement can decrease the cost function. This occurs after a finite number of iterations because SBR is a descent algorithm and there are a finite number of possible subsets S ⊂ {1, . . . , n}. In the limit case λ = 0, we have Ĵ (S; 0) = E(S).

Only insertions can be performed since any removal increases the squared error E(S). SBR coincides with the well-known OLS algorithm [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF]. Generally, the n replacement trials necessitate to compute E(S + {i})

for all insertion trials and E(S -{i}) for all removals. In [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF], we proposed a fast and stable recursive implementation based on the Cholesky factorization of the Gram matrix A T S A S when S is modified by one element (where A S stands for the submatrix of A gathering the active columns). SBR is summarized in Tab. I. The optional output parameters ℓ add and δE add are unnecessary in the standard version. Their knowledge will be useful to implement the extended CSBR algorithm.

Let us illustrate the behavior of SBR on a simple example using the geometrical interpretation of Fig. 4, where a single replacement is represented by a vertical displacement (from top to bottom) between the two lines S and S ± {ℓ}. S init = ∅ yields an horizontal line since Ĵ (∅; λ) = y 2 2 does not depend on

B. Principle of the continuation search

Our continuation strategy is inspired by ℓ 1 -homotopy which recursively computes the minimizers of y -Ax 2 2 + λ x 1 when λ is continuously decreasing [START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF]- [START_REF] Efron | Least angle regression[END_REF]. An iteration of ℓ 1 -homotopy consists in two steps:

• Find the next value λ new < λ cur for which the ℓ 1 optimality conditions are violated with the current active set S (λ cur denotes the current value); CSBR follows the same principle. The first step is now related to some local ℓ 0 -optimality conditions, and the second step consists in calling SBR at λ new with the current active set as initial solution; see Fig. 5 for a sketch. A main difference with ℓ 1 -homotopy is that the ℓ 0 solutions are suboptimal, i.e., they are local minimizers of J (x; λ) with respect to x.

1) Local optimality conditions:

Let us first reformulate the stopping conditions of SBR at a given λ.

SBR terminates when a local minimum of Ĵ (S; λ) has been found:

∀i ∈ {1, . . . , n}, Ĵ (S ± {i}; λ) ≥ Ĵ (S; λ). ( 8
)
This condition is illustrated on Fig. 6(a): all lines related to single replacements S ± {i} lay above the black point representing the value of Ĵ (S; λ) for the current λ. By separating the conditions related to insertions S + {i} and removals S -{i}, (8) rereads as the interval condition:

λ ∈ [δE add (S), δE rmv (S)], (9)
where

δE add (S) max i / ∈S E(S) -E(S + {i}) (10a) δE rmv (S) min i∈S E(S -{i}) -E(S) (10b)
refer to the maximum variation of the squared error when an atom is added in the support S (respectively, removed from S). with the initial subset S + {ℓ add }. The line S + {ℓ add } lays below all other lines S + {i} (dotted lines). Here, the λ-axis has been stretched by an arbitrary factor for improved readability. The horizontal length λnew does not match the vertical length δE add (S), as it should without any stretching. The same stretching process will be done in Fig. 7. In practice, λ new can be set to the limit value

2) Violation of the local optimality conditions:

λ new = δE add (S) (12)
provided that S is replaced with S + {ℓ add }.

As illustrated on Fig. 6 Because S and S +{ℓ add } both lead to the same value of Ĵ ( . ; λ new ), the de-selection of ℓ add is forbidden in the first iteration of SBR.

C. CSBR algorithm

CSBR is summarized in Tab. II. The repeated calls to SBR deliver subsets S j for decreasing λ j . As shown on Fig. 5, the solution S j covers the interval (λ j+1 , λ j ]. At the very first iteration, we have S 0 = ∅, and ( 11)-( 12) reread:

ℓ add ∈ arg max i∈{1,...,n} | y, a i | a i 2 and λ 1 = y, a ℓadd 2 a ℓadd 2 2 . ( 13
)
According to Tab. II, CSBR stops when λ j = 0, i.e., the whole domain λ ∈ Ê + has been scanned.

However, this choice may not be appropriate when dealing with noisy data and overcomplete dictionaries.

In such cases, ad hoc early stopping rules can be considered [START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF], [START_REF] Donoho | About SparseLab[END_REF]. A natural rule takes the form λ j ≤ λ stop with λ stop > 0. Alternative rules involve a maximum cardinality (|S j | ≥ k stop ) and/or a minimum squared error (E(S j ) ≤ ε stop ).

Fig. 5 shows a step-by-step illustration with the early stop λ j ≤ λ stop . The initial support S init = {ℓ add } and λ 1 are precomputed in [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF]. In the first call S 1 = SBR(S init ; λ 1 ), a number of single replacements updates S ← S ± {ℓ} are carried out leading to S 1 = S. This process is represented by the plain vertical arrow at λ 1 linking both lines S 0 and S 1 (the line S init is not shown for readability reasons). Once S 1 is obtained, the next value λ 2 is computed. This process is represented by an oblique, dashed arrow joining λ 1 and λ 2 . These two processes are being repeated alternatively at the second and third iterations of CSBR. Finally, CSBR terminates after λ 4 has been computed because λ 4 ≤ λ stop .

IV. ℓ 0 -REGULARIZATION PATH DESCENT (ℓ 0 -PD)

On the theoretical side, the ℓ 0 -penalized regularization path is piecewise constant (Theorem 1). It yields the ℓ 0 curve which is piecewise affine, continuous and concave (Fig. 1). The curve related to the CSBR outputs does not fulfill this property since: (i) there might be jumps in this curve; and (ii) the slope of the line S j is not necessarily increasing with j (see Fig. 5). This motivates us to propose another algorithm whose outputs are consistent with the structural properties of the ℓ 0 -curve.

We propose to gradually update a list S of candidate subsets S j while imposing that the related curve is a concave polygon, obtained as the concave envelope of the set of lines S j (see Fig. 7(a)). The subsets in S are updated so as to decrease at most the concave polygonal curve. In particular, we impose that the least value is λ J+1 = 0, so that the concave envelope is computed over the whole domain λ ∈ Ê + .

A. Descent of the concave polygon

The principle of ℓ 0 -PD is to perform a series of descent steps, where a new candidate subset S new is considered and included in the list S only if the resulting concave polygon can be decreased. This descent test is illustrated on Fig. 7 for two examples (top and bottom subfigures). For each example, the initial polygon is represented in (a). It is updated when its intersection with the line S new is non-empty

0 λ λ inf = δErmv λsup Snew S j Ĵ (S; λ) 0 λ λ j+1 λ j S j Ĵ (S; λ) (a) (b) S new = S j -{ℓ rmv } (c)

B. Selection of the new candidate support

We first need to assign a Boolean label S j .expl to each subset S j . It equals 1 if S j has already been "explored" and 0 otherwise. The following exploration process is being carried out given a subset S = S j : all the possible single replacements S ± {i} are tested. The best insertion ℓ add and removal ℓ rmv are both kept in memory, with ℓ add defined in [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF] and similarly,

ℓ rmv ∈ arg min i∈S {E(S -{i}) -E(S)}. ( 14
)
At any ℓ 0 -PD iteration, the unexplored subset S j of lowest cardinality (i.e., of lowest index j) is selected. ℓ 0 -PD attempts to include S add = S j + {ℓ add } and S rmv = S j -{ℓ rmv } into S, so that the concave polygon can be decreased at most. The CCV Descent procedure (Tab. III) is first called with S new ← S add leading to possible updates of S and λ. It is called again with S new ← S rmv . Fig. 7 illustrates each of these calls: the slope of S new is Sort the subsets in S by increasing cardinality;

|S j | + 1 and |S j | -
Sort λ in the decreasing order; end S j has been explored, the new supports that have been included in S (if any) are tagged as unexplored.

C. ℓ 0 -PD algorithm ℓ 0 -PD is stated in Tab. IV. Initially, S is formed of the empty support S 0 = ∅. The resulting concave polygon is reduced to a single horizontal edge. The corresponding endpoints are λ 1 = 0 and (by extension) λ 0 +∞. In the first iteration, S 0 is explored: the best insertion S add = {ℓ add } is computed in [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF], and included in S during the call to CCV Descent. The updated set S is now composed of S 0 = ∅ (explored) and S 1 = S add (unexplored). The new concave polygon has two edges delimited by λ 2 = 0, λ 1 and λ 0 = +∞, with λ 1 given in [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF]. Generally, either 0, 1, or 2 new unexplored supports S add and S rmv may be included in S at a given iteration while a variable number of supports may be removed from S.

ℓ 0 -PD terminates when all supports in S have been explored. When this occurs, the concave polygon cannot decrease anymore with any single replacement S j ±{i}, with S j ∈ S. Practically, the early stopping rule λ j ≤ λ stop can be adopted, where j denotes the unexplored subset having the least cardinality. This rule ensures that all candidate subsets S j corresponding to the interval (λ stop , +∞) have been explored.

Similar to CSBR, alternative stopping conditions of the form |S j | ≥ k stop or E(S j ) ≤ ε stop can be adopted.

March 18, 2015 DRAFT

D. Fast implementation

The CCV Descent procedure calls the function intersect to compute the intersection between a concave polygon S and a line S new . Lemma 1 states that this intersection is empty in two simple situations. Hence, the call to intersect is not needed in these situations. This implementation detail is omitted in Tab. III for brevity reasons.

Lemma 1 Let S = {S j , j = 0, . . . , J} be a list of supports associated to a continuous, concave polygon λ → min j Ĵ (S j ; λ) with J + 1 edges, delimited by λ = {λ 0 , . . . , λ J+1 }. The following properties hold for all j:

• If δE add (S j ) < λ j+1 , then the line S add = S j + {ℓ add } lays above the current concave polygon.

• If δE rmv (S j ) > λ j , then the line S rmv = S j -{ℓ rmv } lays above the current concave polygon.

March 18, 2015 DRAFT

Proof: We give a sketch of proof using geometrical arguments. Firstly, δE add (S j ) is the λ-value of the intersection point between lines S j and S new = S j + {ℓ add }; see Fig. This proves the first result. A similar sketch applies to the second result.

E. Main differences between CSBR and ℓ 0 -PD

First, we stress that contrary to CSBR, the index j in λ j does not identify with the iteration number anymore for ℓ 0 -PD. Actually, the current iteration of ℓ 0 -PD is related to an edge of the concave polygon, i.e., a whole interval (λ j+1 , λ j ), whereas the current iteration of CSBR is dedicated to a single value λ j which is decreasing when the iteration number j increases.

Second, the computation of the next value λ j+1 ≤ λ j in CSBR is only based on the violation of the lower bound of ( 9), corresponding to atom selections. In ℓ 0 -PD, the upper bound is considered as well. This is the reason why the λ-values are not scanned in a decreasing order anymore. This may improve the very sparse solutions found in the early iterations within an increased computation time, as we will see hereafter.

V. NUMERICAL RESULTS

The algorithms are evaluated on two kinds of problems involving ill-conditioned dictionaries. The behavior of CSBR and ℓ 0 -PD is first analyzed for simple examples. Then, we provide a detailed comparison with other nonconvex algorithms for many scenarii.

A. Two generic problems

The sparse deconvolution problem takes the form y = h * x ⋆ + n where the impulse response h is a Gaussian filter of standard deviation σ, and the noise n is assumed i. )) where σ 2 n is the variance of the noise process n. The jump detection problem is illustrated on Fig. 8(c,d). Here, A is the squared dictionary (m = n) defined by A i,j = 1 if i ≥ j, and 0 otherwise. The atom a j codes for a jump at location j, and x ⋆ j matches the height of the jump. When x ⋆ is k-sparse, Ax ⋆ yields a piecewise constant signal with k pieces, x ⋆ being the first-order derivative of the signal Ax ⋆ .

Both generic problems involve either square or slightly undercomplete dictionaries. The case of overcomplete dictionaries will be discussed as well, e.g., by considering the deconvolution problem with undersampled observations y. The generic problems are already difficult because neighboring columns of A are highly correlated, and a number of fast algorithms that are efficient for well-conditioned dictionaries may fail to recover the support of x ⋆ . The degree of difficulty of the deconvolution problem is controlled by the width σ of the Gaussian impulse response and the sparsity k: for large values of k and/or σ, the Gaussian features resulting from the convolution h * x ⋆ strongly overlap. For the jump detection problem, 2) Model order selection: It may often be useful to select a single solution x j . The proposed algorithms are compatible with most classical methods of model order selection [START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF], [START_REF] Wang | Model selection[END_REF] because they are greedy algorithms. Assuming that the variance of the observation noise is unknown, we distinguish two categories of cost functions for estimation of the order x j 0 = |S j |. The first take the form min j {m log E(S j ) + α|S j |} where α equals 2, log m, and 2 log log m for the Akaike, Minimum Description Length (MDL)

and Hannan and Quinn criteria, respectively [START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF]. The second are cross-validation criteria [START_REF] Wahba | Practical approximate solutions to linear operator equations when the data are noisy[END_REF], [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF]. The sparse approximation framework allows one to derive simplified expressions of the latter up to the storage of intermediate solutions of greedy algorithms for consecutive cardinalities [START_REF] Miller | Subset selection in regression[END_REF], [START_REF] Wang | Model selection[END_REF], [START_REF] Austin | On the relation between sparse reconstruction and parameter estimation with model order selection[END_REF].

For the sparse deconvolution and jump detection problems, we found that the Akaike and cross validation criteria severely over-estimate the expected number of spikes. On the contrary, the MDL criterion yields quite accurate results. We found that the modified MDLc version dedicated to short data records (i.e., when the number of observations is moderately larger than the model order) [START_REF] De Ridder | Modified AIC and MDL model selection criteria for short data records[END_REF] yields the best results for all the scenarii we have tested. It reads: for noisy data, the spikes of smallest amplitudes are drowned in the noise. One cannot expect to detect them.

min j log E(S j ) + log(m)(|S j | + 1) m -|S j | -2 . ( 15
3) Further empirical observations: Fig. 11 is a typical display of the approximate ℓ 0 -curves yielded by CSBR and ℓ 0 -PD. The ℓ 0 -PD curve is structurally continuous and concave whereas for the CSBR curve, there are two kinds of breakpoints depicted with black and white circles. The former are "continuous" breakpoints. This occurs when no single replacement is done during the call to SBR (SBR(S init ; λ j ) returns S j = S init ; see Tab. II). Otherwise, a discontinuity breakpoint (white circle) appears. In Fig. 11, the CSBR and ℓ 0 -PD curves almost coincide for large λ's, where only continuous breakpoints can be observed. For low λ's, the ℓ 0 -PD curve lays below the CSBR curve, and discontinuity breakpoints appear in the latter curve. (c,d)). From iterations 25 to 40, the very sparse solutions previously found (k ≤ 20) are improved as a series of atom de-selections is performed. They are being improved again around iteration 80. On the contrary, the sparsest solutions are never improved with CSBR, which works for decreasing λ's (Figs. 12(a,b)). For ℓ 0 -PD, the early stopping parameter λ stop may have a strong influence on the improvement of the sparsest solutions and the overall computation time. This point will be further discussed below.

C. Extensive comparisons

The proposed algorithms are compared with popular nonconvex algorithms for both problems introduced in subsection V-A with various parameter settings: problem dimension (m, n), ratio m/n, signal-tonoise ratio, cardinality of x ⋆ , and width σ of the Gaussian impulse response for the deconvolution problem.

The settings are listed on Table V for 10 scenarii. Because the proposed algorithms are orthogonal greedy algorithms, they are better suited to problems in which the level of sparsity is moderate to high. We therefore restrict ourselves to the case where k = x ⋆ 0 ≤ 30.

1) Competing algorithms:

We focus on the comparison with algorithms based on nonconvex penalties.

It is indeed increasingly acknowledged that the BPDN estimates are less accurate than sparse approximation estimates based on nonconvex penalties. We do not consider forward greedy algorithms either; we already showed that SBR is (unsurprisingly) more efficient than the simpler OMP and OLS algorithms [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF].

Among the popular nonconvex algorithms, we consider:

1) Iterative Reweighted Least Squares (IRLS) for ℓ q minimization, q < 1 [52];

2) Iterative Reweighted ℓ 1 (IRℓ 1 ) coupled with the penalty log(|x i | + ε) [START_REF] Candès | Enhancing sparsity by reweighted ℓ1 minimization[END_REF], [START_REF] Wipf | Iterative reweighted ℓ1 and ℓ2 methods for finding sparse solutions[END_REF], [START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF]; 3) ℓ 0 penalized least squares for cyclic descent (L0LS-CD) [START_REF] Seneviratne | Sparse coloured system identification with guaranteed stability[END_REF]; 4) Smoothed ℓ 0 (SL0) [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF], [START_REF] Eftekhari | Robust-sl0 for stable sparse representation in noisy settings[END_REF].

We resort to a penalized least-square implementation for all algorithms, the only algorithm directly working with the ℓ 0 penalty being L0LS-CD. We do not consider simpler thresholding algorithms (Iterative Hard Thresholding, CoSaMP, Subspace Pursuit) proposed in the context of compressive sensing since we found that SBR behaves much better than these algorithms for ill-conditioned dictionaries [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF]. We found that L0LS-CD is more efficient than thresholding algorithms. Moreover, the cyclic descent approach is becoming very popular in the recent sparse approximation literature [START_REF] Mazumder | SparseNet: Coordinate descent with nonconvex penalties[END_REF], [START_REF] Marjanovic | lq sparsity penalized linear regression with cyclic descent[END_REF] although its speed of convergence is sensitive to the quality of the initial solution. Here, we use the BPDN initial solution

arg min x { y -Ax 2 2 + µ x 1 }
where µ is set to half of the maximum tested λ-value (more details will be given hereafter). This simple ad hoc setting allows us to get a rough initial solution that is nonzero and very sparse within a fast computation time.

The three other considered algorithms work with sparsity measures depending on an arbitrary parameter.

Regarding IRLS, we set q = 0.5 or 0.1 as suggested in [START_REF] Lai | Improved iteratively reweighted least squares for unconstrained smoothed ℓq minimization[END_REF]. We chose to run IRLS twice, with q = 0.5 and then q = 0.1 (with the previous output at q = 0.5 as initial solution) so that IRLS is less sensitive to local solutions at q = 0.1. SL0 is a GNC-like algorithm working for increasingly non-convex penalties [START_REF] Friedman | Multivariate adaptive regression splines[END_REF] DRAFT (i.e., Gaussian functions of decreasing widths). For simplicity reasons, we set the lowest width relative to the knowledge of the smallest nonzero amplitude of the ground truth solution x ⋆ . The basic SL0 implementation is dedicated to noise-free problems [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF]. There exist several adaptations in the noisy setting [START_REF] Eftekhari | Robust-sl0 for stable sparse representation in noisy settings[END_REF], [START_REF] Ye | Sparse channel estimation of mimo-ofdm systems with unconstrained smoothed l0-norm-regularized least squares compressed sensing[END_REF] including the precursory work [START_REF] Saito | Superresolution of noisy band-limited data by data adaptive regularization and its application to seismic trace inversion[END_REF]. We chose the efficient implementation of [START_REF] Ye | Sparse channel estimation of mimo-ofdm systems with unconstrained smoothed l0-norm-regularized least squares compressed sensing[END_REF] in which the original pseudo-inverse calculations are replaced by a quasi-Newton strategy using limited memory BFGS updates. Finally, the IRℓ 1 implementation depends on both the choice of parameter ε (which controls the degree of nonconvexity) and the ℓ 1 solver. We have tested two ℓ 1 solvers: the incrowd algorithm [START_REF] Gill | The in-crowd algorithm for fast basis pursuit denoising[END_REF] together with an empirical setting of ε > 0, and ℓ 1 homotopy in the limit case ε → 0, following [START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF]. We found that ℓ 1 homotopy is faster than in-crowd, mainly because the Matlab implementation of in-crowd (provided by the authors) makes calls to the quadprog built-in function, which is computationally expensive for large dimension problems.

2) Numerical protocol: Because the competing algorithms work for a single λ we need to define a grid, denoted by {λ G i , i = 1, . . . , N λ }, for comparison purposes. Such grid is defined in logscale for each of the 10 scenarii (k, A, SNR) defined in Table V. The number of grid points is N λ = 11. For a given scenario, T = 30 trials are being performed in which k-sparse vectors x ⋆ (t) and noise vector n(t) are randomly drawn. This leads us to simulate T observation vectors y(t) = Ax ⋆ (t) + n(t) with t ∈ {1, . . . , T }. Specifically, the location of the nonzero amplitudes in x ⋆ (t) are uniformly distributed and the amplitude values are drawn according to an i.i.d. Gaussian distribution. For each trial t, all competing algorithms need to be run N λ times with y(t) and λ G i as inputs whereas CSBR and ℓ 0 -PD are run only once since they deliver estimates for a continuum of values of λ. Their solution for each λ G i directly deduces from their set of output supports and the knowledge of both breakpoints surrounding λ G i . The algorithms are first evaluated in the optimization viewpoint: the related criteria are their capacity to reach a low value of J (x; λ) and the corresponding CPU time. In this viewpoint, the proposed methods might be somehow favored since they are more directly designed with the criterion J (x; λ) in mind.

On the other hand, J (x; λ) appears to be a natural indicator because solving either ℓ 0 -minimization problem (1), ( 2) or ( 3) is the ultimate goal of any sparse approximation method. As detailed below, some post-processing will be applied to the outputs of algorithms that do not rely on the ℓ 0 -norm so that they are not strongly disadvantaged. Practically, we store the value of J (x; λ G i ) found for each trial and each λ G i . Averaging this value over the trials t yields a table TabJ(a, λ G i ) where a denotes a candidate algorithm. Similarly, the CPU time is averaged over the trials t, leading to another table TabCPU(a, λ G i ). Each table is represented separately as a 2D plot with a specific color for each algorithm: see, e.g., continuum of λ's, and the CPU time is computed only once.

The algorithms are also evaluated in terms of support recovery accuracy. For this purpose, let us first define the "support error" as the minimum over i of the distance

|S ⋆ (t)\S(t, a, λ G i )| + |S(t, a, λ G i )\S ⋆ (t)| (16)
between the support S ⋆ (t) of the unknown sparse vector x ⋆ (t) and the support S(t, a, λ G i ) of the sparse reconstruction at λ G i with algorithm a. (16) takes into account both numbers of false negatives |S ⋆ (t)\S(t, a, λ G i )| and of false positives |S(t, a, λ G i )\S ⋆ (t)|. Denoting by S(t, a, λ G opt ) ← S(t, a, λ G i ) the solution support that is the closest to S ⋆ (t) according to [START_REF] Berk | Forward and backward stepping in variable selection[END_REF], we further consider the number of true positives in S(t, a, λ G opt ), defined as |S ⋆ (t) ∩ S(t, a, λ G opt )|. We will thus report:

• the support error;

• the corresponding number of true positives; The underlying idea in this analysis is that when SE is small (respectively, TP is high), the algorithms are likely to perform well provided that λ is appropriately chosen. However, in practical applications, only one estimate is selected using a suitable model selection criterion. We therefore provide additional evaluations of the MDLc estimate accuracy. For CSBR and ℓ 0 -PD, all output supports are considered to compute the MDLc estimate as described in subsection V-B. For other algorithms, it is equal to one of the sparse reconstructions obtained at λ G i for i ∈ {1, . . . , N λ }. The same three measures as above are computed for the MDLc estimate and averaged over T trials. They are denoted by MDLc-SE(a), MDLc-TP(a) and MDLc-Order(a).

3) Technical adaptations for comparison purposes: Because IRLS and SL0 do not deliver sparse vectors in the strict sense, it is necessary to sparsify their outputs before computing their SE(a) score. This is done by running one iteration of cyclic descent (L0LS-CD): most small nonzero amplitudes are then thresholded to 0. Regarding the values of J (x; λ), a post-processing is performed for algorithms that do not rely on the ℓ 0 -norm. This post-processing can be interpreted as a local descent of J (x; λ).

It consists in:

(i) running one iteration of cyclic descent (L0LS-CD); (ii) computing the squared error related to the output support. L0LS-CD is indeed a local descent algorithm dedicated to J (x; λ) but the convergence towards a least-square minimizer is not reached in one iteration. 4) Analysis in the optimization viewpoint: CSBR and ℓ 0 -PD are always among the most accurate to minimize the cost function, as illustrated on Figs. 13, 14 and 15. We can clearly distinguish two groups of algorithms on these figures: IRLS, L0LS-CD and SL0 one the one hand, and the OLS-based algorithms (SBR, CSBR, ℓ 0 -PD) and IRℓ 1 on the other hand, which are the most accurate. We cannot clearly discriminate the accuracy of SBR and CSBR: one may behave slightly better than the other depending on the scenarii. On the contrary, SBR and CSBR are often outperformed by ℓ 0 -PD. The obvious advantage of CSBR and ℓ 0 -PD over SBR and IRℓ 1 is that they are ℓ 0 -homotopy algorithms, i.e., a set of solutions are delivered for many sparsity levels, and the corresponding λ-values are adaptively found. On the contrary, the SBR output is related to a single λ whose tuning may be tricky. Another advantage over IRℓ 1 is that the structure of forward-backward algorithms is simpler, as no call to any ℓ 1 solver is required. Moreover, the number of parameters to tune is lower: there is a single (early) stopping parameter λ stop . [START_REF] Friedman | Multivariate adaptive regression splines[END_REF] DRAFT min x J (x; λ) CPU Time (seconds) min x J (x; λ) CPU Time (seconds) The price to pay for a better performance is an increase of the computation burden. On Figs. [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF][START_REF] Herzet | Bayesian pursuit algorithms[END_REF] and 15, two lines are drawn for CSBR (respectively, for ℓ 0 -PD). They are horizontal because the algorithm is run only once per trial, so there is a single computation time measurement. The first line corresponds to the overall computation time, i.e., from the start to the termination of CSBR / ℓ 0 -PD. This time is often more expensive than for other algorithms. However, the latter times refer to a single execution for some λ G i value. If one wants to recover sparse solutions for many λ G i 's, they must be cumulated. This is the reason why we have drawn a second line for CSBR and ℓ 0 -PD corresponding to a normalization (by N λ = 11) of the overall computation time. In this viewpoint, the CPU time of CSBR and ℓ 0 -PD are very reasonable.

(E) (E) (F) (F) (G) (G)
The computation time depends on many factors among which the implementation of algorithms (including the memory storage) and the chosen stopping rules. We have followed an homogeneous implementation of algorithms to make the CPU time comparisons meaningful. We have defined two sets of stopping rules depending on the problem dimension. The default parameters apply to medium size problems (m = 300). They are relaxed for problems of larger dimension (m > 500) to avoid huge computational costs. The stopping rule of CSBR and ℓ 0 -PD is always λ ≤ λ stop = αλ G 1 with α = 1 for CSBR and 0.5 (medium size) or 0.8 (large size) for ℓ 0 -PD. For L0LS-CD, the maximum number of cyclic descents (update of every amplitude x i ) is set to 60 or 10 depending on the dimension. For SL0, we have followed the default setting of [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF] for the rate of deformation of the nonconvex penalty.

The number of BFGS iterations done in the local minimization steps for each penalty is set to L = 40 or 5. It is set to 5L for the last penalty which is the most nonconvex. Regarding IRLS and IRLℓ 1 , we keep the same settings whatever the dimension since the computation times remain reasonable for large dimensions. Finally, SBR does not require any arbitrary stopping rule. The problems of large dimensions correspond to scenarii C and D. We observe on Fig. 13 several of the considered algorithms. Here, CSBR and ℓ 0 -PD are among the best algorithms in terms of support error.

The results of Tab. VII and Fig. 15 correspond to the deconvolution problem in noise-free case. The data y are undersampled so that the dictionary A is overcomplete. The undersampling rate ∆ ≈ m/n is set to 2 in scenarii H and I and 4 in scenario J. Again, CSBR and ℓ 0 -PD are among the best (SE, TP, MDLc-order) especially for the most difficult problem J.

6) Overcomplete dictionaries with noise:

We now provide arguments indicating that the proposed algorithms are competitive as well for noisy problems with overcomplete dictionaries. The detailed experiments commented below are not reported for space reasons.

We have first considered the noisy deconvolution problem with ∆ = 2 or 4 leading to overcomplete dictionaries, the other parameters being set as in scenarii A to D. Although the data approximation is March 18, 2015 DRAFT over SBR is that sparse solutions are provided for a continuum of λ-values, enabling the utilization of any classical order selection method. We found that the MDL criterion yields very accurate estimates of the cardinality x 0 .

Our perspectives include the proposal of forward-backward search algorithms that will be faster than SBR and potentially more efficient. In the standard version of SBR, CSBR and ℓ 0 -PD, a single replacement refers to the insertion or removal of a dictionary element. The cost of an iteration is essentially related to the n linear system resolutions done to test single replacements for all dictionary atoms. The proposed algorithms obviously remain valid when working with a larger neighborhood, e.g., when testing the replacement of two atoms simultaneously, but their complexity becomes huge. To avoid such numerical explosion, one may rather choose not to carry out all the replacement tests, but only some tests that are likely to be effective. Extensions of OMP and OLS were recently proposed in this spirit [START_REF] Chatterjee | Projection-based and look-ahead strategies for atom selection[END_REF] and deserve consideration for proposing efficient forward-backward algorithms.

APPENDIX A PROPERTIES OF THE ℓ 0 REGULARIZATION PATHS In this appendix, we prove that the ℓ 0 -penalized path S ⋆ P (see Definition 2) is piecewise constant (Theorem 1) and is a subset of the ℓ 0 -constrained regularization path S ⋆ C (Theorem 2). We will denote the ℓ 0 -curve by λ → J ⋆ (λ) = min S { Ĵ (S; λ)}. Let us recall that this function is concave and affine on each interval (λ ⋆ i+1 , λ ⋆ i ), with i ∈ {0, . . . , I} (Definition 1). Moreover, λ ⋆ I+1 = 0 and λ ⋆ 0 = +∞.

A. Proof of Theorem 1

We prove Theorem 1 together with the following lemma, which is informative about the content of S ⋆ P (λ) for the breakpoints λ = λ ⋆ i .

Lemma 2 Let i ∈ {1, . . . , I -1}. Then, for all λ ∈ (λ ⋆ i+1 , λ ⋆ i ), S ⋆ P (λ) ⊂ S ⋆ P (λ ⋆ i+1 ) ∩ S ⋆ P (λ ⋆ i ). For the first and last intervals, we have:

• For all λ ∈ (0, λ ⋆ I ), S ⋆ P (λ) ⊂ S ⋆ P (λ ⋆ I ).

• For all λ ∈ (λ ⋆ 1 , +∞), S ⋆ P (λ) = {∅} ⊂ S ⋆ P (λ ⋆ 1 ).

Proof of Theorem 1: By definition, the ℓ 0 -curve is the concave envelope of the (finite) set of lines S for all possible subsets S. Because it is affine on the i-th interval (λ ⋆ i+1 , λ ⋆ i ), J ⋆ (λ) coincides with Ĵ (S i ; λ) = E(S i ) + λ|S i |, where S i is some optimal subset for all λ ∈ (λ ⋆ i+1 , λ ⋆ i ).

Let λ ∈ (λ ⋆ i+1 , λ ⋆ i ) and S ∈ S ⋆ P (λ). Then, Ĵ (S; λ) = Ĵ (S i ; λ). It follows that both lines S and S i necessarily coincide; otherwise, they would intersect at λ, and line S would lay below S i on either interval (λ ⋆ i+1 , λ) or (λ, λ ⋆ i ), which contradicts the definition of S i . We conclude that S ∈ S ⋆ P (λ ′ ) for all λ ′ ∈ (λ ⋆ i+1 , λ ⋆ i ). We have shown that the content of S ⋆ P (λ) does not depend on λ when λ ∈ (λ ⋆ i+1 , λ ⋆ i ). Proof of Lemma 2: The first result S ⋆ P (λ) ⊂ S ⋆ P (λ ⋆ i+1 ) ∩ S ⋆ P (λ ⋆ i ) is obtained by slightly adapting the proof of Theorem 1: replace (λ ⋆ i+1 , λ ⋆ i ) by the closed interval [λ ⋆ i+1 , λ ⋆ i ], and set λ ′ to both endpoints of this interval.

The second and third results are obtained similarly, by considering the intervals (0, λ ⋆ I ] and [λ ⋆ 1 , +∞), and setting λ ′ ← λ ⋆ I and λ ′ ← λ ⋆ 1 , respectively. It is obvious that S ⋆ P (λ) reduces to the empty support for λ > λ ⋆ 1 since the ℓ 0 -curve is constant for λ > λ ⋆ 1 .

B. Proof of Theorem 2

The first result is straightforward: for any λ and for S ∈ S ⋆ P (λ), we have S ∈ S ⋆ C (|S|). Otherwise, there would exist S ′ with |S ′ | ≤ |S| and E(S ′ ) < E(S). Then, Ĵ (S ′ ; λ) < Ĵ (S; λ) would contradict S ∈ S ⋆ P (λ). To prove the second result, let us first show that for any i, ∃k i : ∀λ ∈ (λ ⋆ i+1 , λ ⋆ i ), S ⋆ P (λ) ⊂ S ⋆ C (k i ). Let S ∈ S ⋆ P (λ) for some λ ∈ (λ ⋆ i+1 , λ ⋆ i ). Theorem 1 implies that S ∈ S ⋆ P (λ) for any λ ∈ (λ ⋆ i+1 , λ ⋆ i ). Therefore, J ⋆ (λ) = Ĵ (S; λ) for λ ∈ (λ ⋆ i+1 , λ ⋆ i ), and the slope of line S, i.e., |S|, is constant whatever S ∈ S ⋆ P (λ) and λ ∈ (λ ⋆ i+1 , λ ⋆ i ). Let us denote this constant by k i = |S|. According to the first paragraph of the proof, S ∈ S ⋆ P (λ) implies that S ∈ S ⋆ C (k i ). Let us prove the reverse inclusion S ⋆ C (k i ) ⊂ S ⋆ P (λ). Let λ ∈ (λ ⋆ i+1 , λ ⋆ i ) and S ∈ S ⋆ C (k i ). First, we have |S| ≤ k i . Second, for any S ′ ∈ S ⋆ P (λ), we have |S ′ | = k i by definition of k i . We also have that E(S ′ ) = E(S) because S ⋆ P (λ) ⊂ S ⋆ C (k i ). Finally, Ĵ (S ′ ; λ) ≥ Ĵ (S; λ). S ′ ∈ S ⋆ P (λ) implies that S ∈ S ⋆ P (λ). This concludes the proof of the second result.

Fig. 3 .

 3 Fig. 3. Notations relative to our heuristic search algorithms. Their outputs are: (i) a sequence of values λj sorted in the decreasing order; (ii) as many supports Sj, Sj being the solution associated to all λ ∈ (λj+1, λj). By extension, S0 = ∅ for λ > λ1.

λ.

  At the first SBR iteration, a new dictionary atom ℓ = a is selected. The line related to the updated support S ← {a} is of slope |S| = 1. Similarly, some new dictionary atoms b and c are being selected in the next iterations, yielding the supports S ← {a, b} and S ← {a, b, c}. On Fig.4, the dotted lines related to the latter supports have slopes equal to 2 and 3. At iteration 4, the single best replacement is the removal ℓ = a. The resulting support S ← {b, c} is of cardinality 2, and the related line is parallel to the line {a, b} found at iteration 2. During the fifth iteration, none of the n single replacements decreases Ĵ ({b, c}; λ). SBR stops with output S = {b, c}.

• 2 2Fig. 4 .

 24 Fig. 4. Step-by-step illustration of the call S = SBR(∅; λ). Each single replacement is represented by a vertical displacement (from top to bottom) from lines S to S ± {ℓ}. The symbols '+' and '-' respectively refer to the selection and de-selection of atoms a, b and c. Four SBR iterations are done from the initial support Sinit = ∅: the selection of a, b and c, and the de-selection of a. The final output S ← {b, c} is of cardinality 2.

Fig. 5 .

 5 Fig. 5. Step-by-step illustration of CSBR with the early stopping condition λj ≤ λstop. The initial support is S0 = ∅. SBR is called for three decreasing values (plain vertical arrows), with output Sj at λj. The search for the next value λj+1 is represented by an oblique displacement along the line Sj .

Fig. 6 .

 6 Fig. 6. Termination of SBR and next call to SBR. (a) When SBR terminates, no single replacement S ± {i} can decrease Ĵ (S; λ). The dotted lines S + {i} (of slope |S| + 1) lay above the black point (λ, Ĵ (S; λ)). Similarly, all lines S -{i}, of slope |S| -1, lay above this point. (b) Here, S is the SBR output at λcur. The next call to SBR is done at λnew = δE add (S)

  Consider the current output S = SBR(S init ; λ cur ). The local optimality condition (9) is then met for λ = λ cur , but also for any λ ∈ [δE add (S), λ cur ]. The new value for which (9) is violated is λ new = δE add (S)c where c > 0 is arbitrarily small. The violation occurs for i = ℓ add , with ℓ add ∈ arg max i / ∈S {E(S) -E(S + {i})}.

  (b), the line S + {ℓ add } lays below all other parallel lines S + {i}. It intersects line S at λ new . The vertical arrow represents the new call to SBR with inputs S + {ℓ add } and λ new .

  (b). The new concave polygon (c) is obtained as the concave envelope of the former polygon and the line S new . All subsets in S whose edges lay above the line S new are removed from S.This procedure is formally presented in Tab. III. Let us now specify how the new candidate subsets S new are built. S new = S j + {ℓ add }

Fig. 7 .

 7 Fig. 7. ℓ0-PD algorithm: descent of the concave polygon when a new support Snew = Sj + {ℓ add } (top) or Snew = Sj -{ℓrmv} (bottom) is included. (a) Initial configuration. (b) The intersection with line Snew is computed. This yields an interval [λ inf , λsup] for which Snew lays below the concave polygon. (c) When this interval is non-empty, the supports Sj whose related edges lay above the line Snew are removed while Snew is included in S. The list of values λj (corresponding to the vertices of the new concave polygon) is being updated. Their number may either increase or decrease.

  7(b). Secondly, we notice that|S j | ≤ |S add | ≤ |S j+1 | because the concave polygon is concave and |S add | = |S j | + 1.It follows from these two facts that if δE add (S j ) < λ j+1 , the line S add lays above S j+1 for λ ≤ λ j+1 , and above S j for λ ≥ λ j+1 .

Fig. 8 .

 8 Fig. 8. Generic deconvolution (a,b) and jump detection (c,d) problems. The data vectors y and the k nonzero amplitudes of x ⋆ are represented in plain lines and with small circles, respectively. (a) Sparse deconvolution problem with k = 10 spikes, SNR = 25 dB, σ = 24 (Gaussian impulse response), and m = 900, n = 756 (size of dictionary A). (b) Sparse deconvolution problem with k = 30, SNR = 10 dB, σ = 3, m = 300, n = 252. (c) Jump detection problem with k = 30, SNR = 25 dB, m = n = 300. (d) Jump detection problem with k = 10, SNR = 10 dB, m = n = 300.

4 Fig. 9 .

 49 Fig. 9. Jump detection example: processing of the data of Fig. 8(c) using CSBR. Three sparse solutions xj are shown, each being related to some CSBR output Sj, with xj 0 = |Sj |. The original vector y is represented in dashed lines and the approximation Axj is in solid line.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Model order selection using MDLc: display of the selected sparse solution xj and the related data approximation signal. The data of Fig. 8(c,b) (k = 30 true spikes) are processed using CSBR and ℓ0-PD, respectively. (a) corresponds to the simulation shown on Fig. 9. The MDLc solution is the CSBR output support S27 of cardinality 27. (b) is related to the ℓ0-PD output S15, with |S15| = 16.

)Fig. 10 (

 10 Fig. 10(a) illustrates that the number of spikes found using MDLc is very accurate for high SNRs (27 spikes are found, the unknown order being 30). It is underestimated for low SNRs: 16 spikes are found (instead of 30) for the simulation of Fig. 10(b) where SNR = 10 dB. This behavior is relevant because

Fig. 12

 12 Fig. 12 provides some insight on the CSBR and ℓ 0 -PD iterations for a sparse deconvolution problem with x ⋆ 0 = 17 and SNR = 20 dB. In the CSBR subfigures, the horizontal axis represents the number of single replacements: 60 replacements are being performed from the initial empty support during the successive calls to SBR. For ℓ 0 -PD, the horizontal axis shows the iteration number. At most two new supports are being included in the list of candidate subsets at each iteration. The number of effective single replacements is therefore increased by 0, 1 or 2. During the first 25 iterations, ℓ 0 -PD mainly operates atom selections similar to CSBR. The explored subsets are thus of increasing cardinality and λ is decreasing (Figs. 12(c,d)). From iterations 25 to 40, the very sparse solutions previously found

Fig. 12 .

 12 Fig. 12. Series of single replacements performed by CSBR and ℓ0-PD. (a) CSBR: cardinality of the current support found after each single replacement. (b) Breakpoints λj found by CSBR, represented in log-scale. SBR is executed for each λj, and the number of single replacements for fixed λj matches the length of the horizontal steps in the figure. (c) ℓ0-PD: cardinality of the supports appended to the regularization path during the iterations. At each iteration, 0, 1 or 2 supports are included. Vertical steps appear whenever two supports are simultaneously included. (d) ℓ0-PD: representation in log-scale of the current interval (λj+1, λj) (grey color). When the grey bars reach the bottom of the image, the lower bound equals λj+1 = 0.

Fig. 13 .Fig. 13 .

 1313 Fig.[START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF]. CSBR and ℓ 0 -PD are represented with continuous curves because J (x; λ) is computed for a

•

  the corresponding model order |S(t, a, λ G opt )|. Averaging these measures over T trials yields the support error score SE(a), the true positive score TP(a) and the model order, denoted by Order(a). The numbers of false positives (FP) and of true/false negatives can be directly deduced, e.g., FP(a) = Order(a) -TP(a).

Fig. 14 .

 14 Fig. 14. Comparison of algorithms for the jump detection problem for the scenarii E, F and G of Table V.

Fig. 15 .

 15 Fig. 15. Comparison of algorithms for the noise-free deconvolution problem, i.e., for the scenarii H, I and J of Table V. Some markers do not appear for low λ's (L0LS-CD, SL0) in the left figures because they do not lay in the zoom-in window (their performance is poor).

  Fig. 2. Sparse approximation seen as a bi-objective optimization problem. The Pareto frontier gathers the non-dominated points: no other point can strictly decrease both |S| and E (S). Bullets and squares are all Pareto solutions. A supported solution is a

	E (S)
				S ′
	E (S ⋆a ) E (S ⋆b )	S ⋆a	S ⋆b
	E (S ⋆c )			S ⋆c
		k	k + 1	k + 2	|S|
	March 18, 2015				DRAFT

minimizer of E (S) + λ|S| with respect to S for some λ. S ⋆a and S ⋆c are supported, contrary to S ⋆b .

TABLE I SBR

 I ALGORITHM FOR MINIMIZATION OF Ĵ (S; λ) FOR FIXED λ [3]. BY DEFAULT, Sinit = ∅. THE OUTPUTS δE add AND ℓ add ARE OPTIONAL. THE SINGLE REPLACEMENT TESTS APPEAR IN THE FOR LOOP.

inputs : A, y, λ, Sinit outputs: S, δE add , ℓ add S best ← Sinit; repeat S ← S best ; for i = 1 to n do Compute Ĵ (S ± {i}; λ); end S best ← S ± {ℓ} with ℓ computed from (7); until Ĵ (S best ; λ) ≥ Ĵ (S; λ); Compute ℓ add according to (11); Set δE add = E (S) -E (S + {ℓ add });

TABLE II CSBR

 II ALGORITHM: SBR IS CALLED REPEATEDLY FOR DECREASING λj 'S. AT ITERATION j, BOTH THE NEXT VALUE λj+1 AND THE NEXT INITIAL SUBSET Sj + {ℓ add } ARE PROVIDED AS SBR OUTPUTS.

	Compute λ1 according to (13);
	j ← 1;
	while λj > 0 do
	Call [Sj , δE

inputs : A, y outputs: S: list of supports Sj; λ: list of λj S0 ← ∅; Sinit ← {ℓ add } with ℓ add computed from (13); add , ℓ add ] = SBR(Sinit; λj); Sinit ← Sj + {ℓ add }; λj+1 ← δE add ; j ← j + 1; end

TABLE III CONCAVE

 III POLYGON DESCENT PROCEDURE. WHEN A NEW SUBSET IS INCLUDED, BOTH LISTS S AND λ ARE UPDATED. THE FUNCTION intersect COMPUTES THE INTERSECTION BETWEEN A LINE AND A CONCAVE POLYGON. THIS YIELDS AN INTERVAL [λ inf , λsup]. BY CONVENTION, λ inf > λsup WHEN THE INTERSECTION IS EMPTY. Remove any subset Sj from S such that [λj+1, λj] ⊂ [λ inf , λsup];

Procedure: CCV Descent(S, Snew, λ)

Call [λ inf , λsup] = intersect(S, Snew); if λ inf < λsup then

Include Snew as an unexplored support in S;

TABLE IV ℓ0

 IV -PD ALGORITHM. THE ALGORITHM MAINTAINS A LIST S OF SUPPORTS Sj WHOSE CARDINALITY IS INCREASING WITH j. THE UNEXPLORED SUPPORT HAVING THE LOWEST CARDINALITY IS EXPLORED AT EACH ITERATION. THE LISTS S AND λ ARE UPDATED DURING THE CALLS TO CCV Descent; λ IS SORTED IN THE DECREASING ORDER, WITH λJ+1 = 0. DURING THE FIRST ITERATION, j = 0 LEADS TO Srmv ← ∅.

	inputs : A, y
	outputs: S, λ
	λ ← {λ1} with λ1 ← 0;
	S0 ← ∅, S0.expl ← 0;
	S ← {S0};

while {∃j : Sj.expl = 0} do Set j as the lowest index such that Sj.expl = 0; Sj.expl ← 1; Compute S add ← Sj + {ℓ add } from (11); if j = 0 then Srmv ← ∅; else Compute Srmv ← Sj -{ℓrmv} from (14); end Call CCV Descent(S, S add , λ); Call CCV Descent(S, Srmv, λ); end

TABLE V SETTINGS

 V RELATED TO EACH SCENARIO: k IS THE SPARSITY. f CONTROLS THE DICTIONARY SIZE: m = f mDEF, n = f nDEF WITH nDEF ≈ mDEF = 300. BY DEFAULT, f = 1. THE UNDERSAMPLING PARAMETER ∆ EQUALS 1 BY DEFAULT (mDEF ≥ nDEF). IT IS INCREASED TO GENERATE PROBLEMS WITH OVERCOMPLETE DICTIONARIES (m ≈ n/∆). THE GAUSSIAN IMPULSE RESPONSE WIDTH IS SET TO σ = f σDEF WITH σDEF = 3 OR 8.

	Scenario	Type	SNR	k f ∆	m	n	σ
	A	Deconv.	25 30 1	1	300	282	3
	B	Deconv.	10 10 1	1	300	252	8
	C	Deconv.	25 10 3	1	900	756 24
	D	Deconv.	25 30 6	1 1800 1692 18
	E	Jumps	25 10 1	1	300	300	∅
	F	Jumps	25 30 1	1	300	300	∅
	G	Jumps	10 10 1	1	300	300	∅
	H	Deconv. +∞ 10 3	2	450	756 24
	I	Deconv. +∞ 30 3	2	450	756 24
	J	Deconv. +∞ 10 1	4	75	252	8

TABLE VI JUMP

 VI DETECTION PROBLEM IN THE NOISY SETTING. THE ALGORITHMS ARE EVALUATED IN TERMS OF SUPPORT ERROR (SE) AND NUMBER OF TRUE POSITIVES (TP). THE NUMBER OF JUMPS THAT ARE FOUND IS REPORTED (ORDER) TOGETHER WITH THE "TRUE ORDER" CORRESPONDING TO THE GROUND TRUTH k. THE SCORES RELATED TO THE MDLC ESTIMATE ARE INDICATED SIMILARLY.

	Scenario E	ℓ0-PD CSBR SBR ℓ0LS-CD Sℓ0 IRℓ1 IRLS
	SE	1.6	1.6 1.6	5.3 4.0 1.5	1.8
	TP	8.6	8.7 8.6	5.2 7.8 8.7	8.7
	Order (true: 10)	8.8	9.0 8.9	5.7 9.6 8.9	9.1
	MDLc-SE	4.7	4.3 4.1	22.7 5.6 4.1	3.5
	MDLc-TP	8.7	8.8 8.8	6.9 8.6 8.8	8.8
	MDLc-Order	12.2	11.9 11.6	26.6 12.7 11.7 11.0
	Scenario F	ℓ0-PD CSBR SBR ℓ0LS-CD Sℓ0 IRℓ1 IRLS
	SE	11.1	11.9 11.8	22.5 11.6 10.9 11.6
	TP	21.2	20.6 20.7	9.2 20.3 20.8 20.7
	Order (true: 30)	23.6	23.2 23.2	10.9 22.2 22.5 23.1
	MDLc-SE	13.7	13.4 13.4	39.2 14.0 13.1 13.3
	MDLc-TP	21.8	21.8 21.4	12.9 21.6 22.1 21.5
	MDLc-Order	27.3	27.0 26.3	35.0 27.2 27.2 26.4
	Scenario G	ℓ0-PD CSBR SBR ℓ0LS-CD Sℓ0 IRℓ1 IRLS
	SE	7.3	7.5 7.5	8.9 10.3 7.2	7.5
	TP	4.0	3.6 3.6	3.1 2.9 3.9	4.0
	Order (true: 10)	5.23	4.73 4.73	5.17 6.07 4.97 5.57
	MDLc-SE	11.4	10.7 10.9	11.7 15.1 11.2 10.7
	MDLc-TP	4.2	4.2 4.2	3.0 3.9 4.2	4.4
	MDLc-Order	9.8	9.1 9.3	7.6 12.8 9.6	9.5

TABLE VII SPARSE

 VII DECONVOLUTION PROBLEM IN THE NOISE-FREE SETTING: EXACT SUPPORT RECOVERY.

	Scenario H	ℓ0-PD CSBR SBR ℓ0LS-CD Sℓ0 IRℓ1 IRLS
	SE	2.5	3.6 4.8	11.4 13.0 0.8	6.1
	TP	8.3	8.2 6.8	0.4	0.1 9.5	9.4
	Order (true: 10)	9.1	10.0 8.3	2.2	3.2 9.8 14.9
	MDLc-SE	3.6	3.8 5.8	168.5 343.8 1.1	9.0
	MDLc-TP	8.6	8.6 7.9	3.3	6.6 9.5	9.6
	MDLc-Order	10.8	11.0 11.6	153.5 347.0 10.1 18.2
	Scenario I	ℓ0-PD CSBR SBR ℓ0LS-CD Sℓ0 IRℓ1 IRLS
	SE	0.9	1.3 2.1	36.7 48.5 3.8	9.4
	TP	29.4	29.3 29.1	0.7	0.8 28.0 27.7
	Order (true: 30) 29.7	29.8 30.2	8.2 20.1 29.8 34.8
	MDLc-SE	3.8	3.5 3.7	686.0 444.9 9.5 114.3
	MDLc-TP	29.5	29.4 29.2	28.6 17.5 28.5 26.4
	MDLc-Order	32.8	32.3 32.1	437.0 449.8 36.5 137.2
	Scenario J	ℓ0-PD CSBR SBR ℓ0LS-CD Sℓ0 IRℓ1 IRLS
	SE	0.3	3.5 5.3	10.3 10.4 2.4	4.3
	TP	9.8	7.3 5.8	0.6	2.6 8.8	9.2
	Order (true: 10)	9.8	8.1 6.9	1.4	5.6 10.0 12.7
	MDLc-SE	2.6	7.7 12.4	176.2 78.6 7.2 69.0
	MDLc-TP	9.7	8.9 8.0	8.6	3.0 8.9	4.1
	MDLc-Order	12.0	15.5 18.5	73.0 74.6 14.9 67.2

[START_REF] Friedman | Multivariate adaptive regression splines[END_REF] DRAFT

This work was carried out in part while C. Soussen was visiting IRCCyN during the academic year 2010-2011 with the financial support of CNRS.

qualitatively good for CSBR and ℓ 0 -PD, the SE and TP scores are very weak. It is hard to discriminate the performance of algorithms because these measures are very weak for all considered algorithms. Moreover, the values of J (λ) found for most algorithms are often similar.

We have also considered an adaptive spline approximation problem generalizing the jump detection problem to the approximation of a signal using piecewise polynomials of degree P = 1 or 2 [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF]. The jump detection problem can indeed be thought of as the approximation with a piecewise constant signal (P = 0). [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] is inspired from the regression spline modeling in [START_REF] Friedman | Multivariate adaptive regression splines[END_REF]. Now, the dictionary atoms are related to the detection of the locations of jumps, changes of slopes and changes of curvatures in the signal y (subdictionaries A 0 , A 1 and A 2 ). The dictionary then takes the form

The generalized version

where each sub-dictionary A p (p ≤ P ) is formed of shifted versions of the one-sided power function i → [max(i, 0)] p . The size of the full dictionary A is approximately m × (P + 1)m.

Hence, it becomes overcomplete as soon as P ≥ 1. We have shown [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF] that SBR is competitive when P = 1 or 2. We have carried out new tests confirming that CSBR and ℓ 0 -PD are more efficient than their competitors in terms of values of J (λ). However, the rate of true positives is low for P ≥ 1 since the location of the change of slopes and of curvatures can hardly be exactly recovered from noisy data.

VI. SOFTWARE

The Matlab implementation of the proposed CSBR and ℓ 0 -PD algorithms is available at www.cran. univ-lorraine.fr/perso/charles.soussen/software.html including programs showing how to call these functions.

VII. CONCLUSION

The choice of a relevant sparse approximation algorithm relies on a trade-off between the desired performance and the computation time one is ready to spend. The proposed algorithms are relatively expensive but very well suited to inverse problems inducing highly correlated dictionaries. A reason is that they have the capacity to escape from local minimizers of J (x; λ) = y -Ax 2 2 + λ x 0 [START_REF] Soussen | From Bernoulli-Gaussian deconvolution to sparse signal restoration[END_REF]. This behavior is in contrast with other classical sparse algorithms.

We have shown the usefulness and efficiency of the two SBR extensions when the level of sparsity is moderate to high, i.e., k/ min(m, n) is lower than 0.1. They remain competitive when k/ min(m, n) ranges between 0.1 and 0.2, and their performance gradually degrade for weaker levels of sparsity, which is an expected behavior for such greedy type algorithms. For a single λ, CSBR is as efficient as SBR, and ℓ 0 -PD improves the SBR and CSBR performance within a larger computation cost. The main benefit [START_REF] Friedman | Multivariate adaptive regression splines[END_REF] DRAFT