Charles Soussen

Member, IEEE Jérôme Idier
email: jerome.idier@irccyn.ec-nantes.fr.

David Brie

Junbo Duan
email: jduan@tulane.edu.

From Bernoulli-Gaussian deconvolution to sparse signal restoration

Formulated as a least square problem under an ℓ 0 constraint, sparse signal restoration is a discrete optimization problem, known to be NP complete. Classical algorithms include, by increasing cost and efficiency, Matching Pursuit (MP), Orthogonal Matching Pursuit (OMP), Orthogonal Least Squares (OLS), stepwise regression algorithms and the exhaustive search. We revisit the Single Most Likely Replacement (SMLR) algorithm, developed in the mid-80's for Bernoulli-Gaussian signal restoration. We show that the formulation of sparse signal restoration as a limit case of Bernoulli-Gaussian signal restoration leads to an ℓ 0 -penalized least square minimization problem, to which SMLR can be straightforwardly adapted.

The resulting algorithm, called Single Best Replacement (SBR), can be interpreted as a forward-backward extension of OLS sharing similarities with stepwise regression algorithms. Some structural properties of SBR are put forward. A fast and stable implementation is proposed. The approach is illustrated on two inverse problems involving highly correlated dictionaries. We show that SBR is very competitive with popular sparse algorithms in terms of trade-off between accuracy and computation time.

Index Terms

Sparse signal estimation; inverse problems; Bernoulli-Gaussian signal restoration; SMLR algorithm; mixed ℓ 2 -ℓ 0 criterion minimization; Orthogonal Least Squares; stepwise regression algorithms.

I. INTRODUCTION

Sparse signal restoration arises in inverse problems such as Fourier synthesis, mono-and multidimensional deconvolution, and statistical regression. It consists in the decomposition of a signal y as a linear combination of a limited number of elements from a dictionary A. While formally very similar, sparse signal restoration has to be distinguished from sparse signal approximation. In sparse signal restoration, the choice of the dictionary is imposed by the inverse problem at hand whereas in sparse approximation, the dictionary has to be chosen according to its ability to represent the data with a limited number of coefficients.

Sparse signal restoration can be formulated as the minimization of the squared error y -Ax 2 (where • refers to the Euclidean norm) under the constraint that the ℓ 0 pseudo-norm of x, defined as the number of non-zero entries in x, is small. This problem is often referred to as subset selection because it consists in selecting a subset of columns of A. This yields a discrete problem (since there are a finite number of possible subsets) which is known to be NP-complete [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF]. In this paper, we focus on "difficult" situations in which some of the columns of A are highly correlated, the unknown weight vector x is only approximately sparse, and/or the data are noisy. To address subset selection in a fast and sub-optimal manner, two approaches can be distinguished.

The first one, which has been the most popular in the last decade, approximates the subset selection problem by a continuous optimization problem, convex or not, that is easier to solve [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF][START_REF] Rao | Subset selection in noise based on diversity measure minimization[END_REF][START_REF] Candès | Enhancing sparsity by reweighted ℓ1 minimization[END_REF][START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF][START_REF] Wipf | Iterative reweighted ℓ1 and ℓ2 methods for finding sparse solutions[END_REF]. In particular, the ℓ 1 relaxation of the ℓ 0 -norm has been increasingly investigated [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF], leading to the LASSO optimization problem.

The second approach addresses the exact subset selection problem using either iterative thresholding [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF][START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF][START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF][START_REF] Blumensath | Normalized iterative hard thresholding: Guaranteed stability and performance[END_REF] or greedy search algorithms. The latter gradually increase or decrease by one the set of active columns.

The simplest greedy algorithms are Matching Pursuit (MP) [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] and the improved version Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF]. Both are referred to as forward greedy algorithms since they start from the empty active set and then gradually increase it by one element. In contrast, the backward algorithm of Couvreur and Bresler [START_REF] Couvreur | On the optimallity of the backward greedy algorithm for the subset selection problem[END_REF] starts from a complete active set which is gradually decreased by one element. It is, however, only valid for undercomplete dictionaries. Forward-backward algorithms (also known as stepwise regression algorithms) in which insertions and removals of dictionary elements are both allowed, are known to yield better recovery performance since an early wrong selection can be counteracted by its further removal from the active set [START_REF] Efroymson | Multiple regression analysis[END_REF][START_REF] Berk | Forward and backward stepping in variable selection[END_REF][START_REF] Broersen | Subset regression with stepwise directed search[END_REF][START_REF] Miller | Subset Selection in Regression[END_REF]. In contrast, the insertion of a wrong element is irreversible when using forward algorithms. We refer the reader to [START_REF] Miller | Subset Selection in Regression[END_REF]Chapter 3] for an overview of the forward-backward algorithms in subset selection.

The choice of the algorithm depends on the amount of time available and on the structure of matrix A. In favorable cases, the sub-optimal search algorithms belonging to the first or the second approach provide solutions having the same support as the exhaustive search solution. Specifically, if the unknown signal is highly sparse and if the correlation between any pair of columns of A is low, the ℓ 1 -norm approximation provides optimal solutions [START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF]. But when fast algorithms are unsatisfactory, it is relevant to consider slower algorithms being more accurate and remaining very fast compared to the exhaustive search. The Orthogonal Least Squares algorithm (OLS) [START_REF] Chen | Orthogonal least squares methods and their application to non-linear system identification[END_REF] which is sometimes confused with OMP [START_REF] Blumensath | On the difference between Orthogonal Matching Pursuit and Orthogonal Least Squares[END_REF], falls into this category. Both OLS and OMP share the same structure, the difference being that at each iteration, OLS solves as many least square problems as there are non-active columns while OMP only performs one linear inversion. In this paper, we derive a forward-backward extension of OLS allowing an insertion or a removal per iteration, each iteration requiring to solve n least square problems, where n is the size of x.

The proposed forward-backward extension of OLS can be viewed as a new member of the family of stepwise regression algorithms. The latter family traces back to 1960 [START_REF] Efroymson | Multiple regression analysis[END_REF], and other popular algorithms were proposed in the 1980's [START_REF] Miller | Subset Selection in Regression[END_REF] and more recently [START_REF] Haugland | A Bidirectional Greedy Heuristic for the Subspace Selection Problem[END_REF]. Note that forward-backward extensions of OMP have also been proposed [START_REF] Zhang | Adaptive forward-backward greedy algorithm for learning sparse representations[END_REF][START_REF] Herzet | Bayesian pursuit algorithms[END_REF]. In contrast with the other stepwise regression algorithms, our approach relies on a bi-objective formulation in order to handle the trade-off between low residual and low cardinality. This formulation reads as the minimization of the ℓ 0 -penalized least square cost function y -Ax 2 +λ x 0 . Then, we design a heuristic algorithm to minimize this cost function in a suboptimal way. While the other forward-backward strategies [15-17, 21, 22] aim at handling the same trade-off, most of them are not expressed as optimization algorithms, but rather as empirical schemes without any connexion with an objective function. Moreover, some of them involve discrete search parameters that control variable selection or de-selection [START_REF] Efroymson | Multiple regression analysis[END_REF][START_REF] Berk | Forward and backward stepping in variable selection[END_REF][START_REF] Zhang | Adaptive forward-backward greedy algorithm for learning sparse representations[END_REF] while others do not involve any parameter [START_REF] Broersen | Subset regression with stepwise directed search[END_REF][START_REF] Haugland | A Bidirectional Greedy Heuristic for the Subspace Selection Problem[END_REF].

An exception can be made for Broersen's algorithm [START_REF] Broersen | Subset regression with stepwise directed search[END_REF] since it aims at minimizing y -Ax 2 + λ x 0 for a specific λ value corresponding to Mallows' C p statistic. However, it is only valid for undercomplete problems. On the contrary, our proposed algorithm is general and valid for any λ value. It does not necessitate to tune any other parameters (e.g., stopping parameters).

Our starting point is the Single Most Likely Replacement (SMLR) algorithm which proved to be a very efficient tool for the deconvolution of a Bernoulli-Gaussian signal [START_REF] Kormylo | Maximum-likelihood detection and estimation of Bernoulli-Gaussian processes[END_REF][START_REF] Mendel | Optimal Seismic Deconvolution[END_REF][START_REF] Goussard | A new algorithm for iterative deconvolution of sparse spike trains[END_REF][START_REF] Champagnat | Unsupervised deconvolution of sparse spike trains using stochastic approximation[END_REF]. We show that sparse signal restoration can be seen as a limit case of maximum a posteriori (MAP) Bernoulli-Gaussian restoration which results in an adaptation of SMLR to subset selection. The paper is organized as follows. In Section II, we introduce the Bernoulli-Gaussian model and the Bayesian framework from which we formulate the sparse signal restoration problem. In Section III, we adapt SMLR resulting in the socalled Single Best Replacement (SBR) algorithm. In Section IV, we propose a fast and stable SBR implementation. Finally, Sections V and VI illustrate the method on the sparse spike deconvolution with a Gaussian impulse response and on the joint detection of discontinuities at different orders in a signal.

II. SPARSE SIGNAL ESTIMATION USING A LIMIT BERNOULLI-GAUSSIAN MODEL

A. Preliminary definitions and working assumptions

Given an observation vector y ∈ Ê m and a dictionary A = [a 1 , . . . , a n ] ∈ Ê m×n , a subset selection algorithm aims at computing a weight vector x ∈ Ê n yielding an accurate approximation y ≈ Ax. The columns a i corresponding to the non-zero weights x i are referred to as the active (or selected) columns.

Throughout this paper, no assumption is made on the size of A: m can be either smaller or larger than n. A is assumed to satisfy the unique representation property (URP): any min(m, n) columns of A are linearly independent. This assumption is usual when m n; it is stronger than the full rank assumption [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using FOCUSS: A re-weighted minimum norm algorithm[END_REF]. When m n, it amounts to the full rank assumption. Although URP was originally introduced to guarantee uniqueness of sparse solutions [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using FOCUSS: A re-weighted minimum norm algorithm[END_REF], we use this assumption to propose a valid algorithm. It can actually be relaxed provided that the search strategy guarantees that the selected columns are linearly independent (see Section VI-C for details).

The support of a vector x ∈ Ê n is the set S(x) ⊆ {1, . . . , n} defined by i ∈ S(x) if and only if

x i = 0. We denote by Q ⊆ {1, . . . , n} the active set and by q ∈ {0, 1} n the related vector defined by

q i = 1 if and only if i ∈ Q. When Card[Q] min(m, n), let A Q be the submatrix of size m × Card[Q]
formed of the active columns of A. We define the least square solution and the related squared error:

x Q arg min S(x)⊆Q {E(x) = y -Ax 2 } (1) E Q E(x Q ) = y -Ax Q 2 .
(2)

B. Bayesian formulation of sparse signal restoration

We consider the restoration of a sparse signal x from a linear observation y = Ax + n, where n stands for the observation noise. An acknowledged probabilistic model dedicated to sparse signals is the Bernoulli-Gaussian (BG) model [START_REF] Kormylo | Maximum-likelihood detection and estimation of Bernoulli-Gaussian processes[END_REF][START_REF] Mendel | Optimal Seismic Deconvolution[END_REF][START_REF] Champagnat | Unsupervised deconvolution of sparse spike trains using stochastic approximation[END_REF]. For such model, deterministic optimization algorithms [START_REF] Champagnat | Unsupervised deconvolution of sparse spike trains using stochastic approximation[END_REF] and Markov chain Monte Carlo techniques [START_REF] Cheng | Simultaneous wavelet estimation and deconvolution of reflection seismic signals[END_REF] are used to compute the MAP and the posterior mean, respectively. Hereafter, we define the BG model and then consider its estimation in the joint MAP sense.

A BG process can be defined using a Bernoulli random vector q ∈ {0, 1} n coding for the support and a Gaussian random vector r ∼ N (0, σ 2 x I n ), with I n the identity matrix of size n. Each sample x i of x is modeled as x i = q i r i [START_REF] Kormylo | Maximum-likelihood detection and estimation of Bernoulli-Gaussian processes[END_REF][START_REF] Mendel | Optimal Seismic Deconvolution[END_REF]. The Bernoulli parameter ρ = Pr(q i = 1) is the probability of presence of signal and σ 2

x controls the variance of the nonzero amplitudes x i = r i . The Bayesian formulation consists in inferring x = (q, r) knowing y. The MAP estimator can be obtained by maximizing the marginal likelihood l(q | y) [START_REF] Champagnat | Unsupervised deconvolution of sparse spike trains using stochastic approximation[END_REF] or the joint likelihood l(q, r | y) [START_REF] Mendel | Optimal Seismic Deconvolution[END_REF][START_REF] Goussard | A new algorithm for iterative deconvolution of sparse spike trains[END_REF]. Following [START_REF] Mendel | Optimal Seismic Deconvolution[END_REF] and assuming a Gaussian white noise n ∼ N (0, σ 2 n I m ), independent from x, Bayes' rule leads to:

L(q, r) -2σ 2 n log[l(q, r | y)] = y -A∆ q r 2 + σ 2 n σ 2 x r 2 + λ q 0 + c (3)
where λ = 2σ 2 n log(1/ρ -1), ∆ q is the diagonal matrix of size n whose diagonal elements are q i (x reads x = ∆ q r), and c is a constant. Now, a signal x is sparse if some entries x i are equal to 0. Since this definition does not impose constraints on the range of the non zero amplitudes, we choose to use a limit Bernoulli-Gaussian model in which the amplitude variance σ 2

x is set to infinity. Note that a parallel limit development was done, independently from our work, in the conference paper [START_REF] Herzet | Bayesian pursuit algorithms[END_REF]. In Appendix A, we show that the minimization of L w.r.t. x = (q, r) rereads:

min x∈Ê n {J (x; λ) = y -Ax 2 + λ x 0 }. ( 4
)
This formulation is close to that obtained in the Bayesian subset selection literature [18, Chapter 7] using an alternative Bernoulli-Gaussian model. In the latter model, the Gaussian prior relies on R Q r instead of r, with R Q the Cholesky factor of the Gram matrix A t Q A Q . This leads to a cost function of the form (4), the difference being that λ depends on the amplitude variance σ 2

x and tends to infinity as σ 2

x tends to infinity [START_REF] George | Calibration and empirical Bayes variable selection[END_REF][START_REF] Chipman | The practical implementation of Bayesian model selection[END_REF].

Remark 1 (Noise-free case) The Bayesian development above is valid for noisy data. In the noise-free case, we define the sparse solution as the limit of arg min x J (x; λ) when λ tends towards 0. According to classical results in optimization [START_REF] Nocedal | Numerical optimization[END_REF]Chapter 17], if {λ k } is a sequence decreasing towards 0 and x k is an exact global minimizer of J (x; λ k ), then every limit point of the sequence {x k } is a solution of arg min x x 0 s.t. y -Ax 2 is minimal. In Appendix B, we derive a more precise result: "the set of minimizers of J (x; λ) is constant when λ is close enough to 0 (λ = 0). It is equal to the set of sparsest solutions to y = Ax in the overcomplete case, and to the unconstrained least-squares solution in the undercomplete case."

In the following, we focus on the minimization problem (4). The hyperparameter λ is fixed. It controls the level of sparsity of the desired solution. The algorithm that will be developed relies on an efficient search of the support of x. The search strategy is based on the definition of a neighborhood relationship between two supports: two supports are neighbors if one is nested inside the other and the largest support has one more element.

III. SINGLE BEST REPLACEMENT ALGORITHM

We propose to adapt the SMLR algorithm to the minimization of the mixed ℓ 2 -ℓ 0 cost function J (x; λ) defined in (4). To clearly distinguish SMLR which specifically aims at minimizing (3), the adapted algorithm will be termed as Single Best Replacement (SBR).

A. Principle of SMLR and main notations

SMLR [START_REF] Kormylo | Maximum-likelihood detection and estimation of Bernoulli-Gaussian processes[END_REF] is a deterministic coordinatewise ascent algorithm to maximize likelihood functions of the form l(q | y) (marginal MAP estimation) or l(q, r | y) (joint MAP estimation). In the latter case, it is easy to check from (3) that given q, the minimizer of L(q, r) w.r.t. r has a closed form expression r = r(q). Consequently, the joint MAP estimation reduces to the minimization of L(q, r(q)) w.r.t. q.

At each SMLR iteration, all the possible single replacements of the support q (set q i = 1q i while keeping the other q j , j = i unchanged) are tested, then the replacement yielding the maximal decrease of L(q, r(q)) is chosen. This task is repeated until no single replacement can decrease L(q, r(q)) anymore.

The number of possible supports q being finite and SMLR being a descent algorithm, it terminates after a finite number of iterations.

Before adapting SMLR, let us introduce some useful notations. We denote by Q•i a single replacement, i.e., an insertion or removal into/from the active set Q:

Q • i    Q ∪ {i} if i / ∈ Q, Q\{i} otherwise.
When Card[Q] min(m, n), we define the cost function

J Q (λ) E Q + λCard[Q] (5)
involving the squared error E Q defined in [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. By definition of J (x

Q ; λ) = E Q + λ x Q 0 , J Q (λ) coincides with J (x Q ; λ) when the support of x Q is equal to Q.
Although it aims at minimizing J (x; λ), the proposed SBR algorithm involves the computation of

J Q (λ) rather than J (x Q ; λ).
We make this choice because J Q (λ) can be computed and updated more efficiently, the computation of x Q being no longer necessary. In subsection III-C, we show that for noisy data, the replacement of J (x Q ; λ) by J Q (λ) has a negligible effect.

B. The Single Best Replacement algorithm

SMLR can be seen as an exploration strategy for discrete optimization rather than an algorithm specific to a posterior likelihood function. Here, we use this strategy to minimize J (x; λ). We rename the algorithm Single Best Replacement to remove any statistical connotation.

SBR works as follows. Consider the current support Q. The n single replacements Q • i are tested, i.e., we compute the squared errors E Q•i and we memorize the values of J Q•i (λ). If the minimum of

J Q•i (λ) is lower than J Q (λ)
, then we select the index yielding this minimum value:

ℓ ∈ arg min i∈{1,...,n} J Q•i (λ). (6)
The next SBR iterate is thus defined as

Q ′ = Q • ℓ. This task is repeated until J Q (λ) cannot decrease
anymore. By default, we use the initial empty support. The algorithm is summarized in Table I.

C. Case where some active amplitudes are zero

We show that this case almost surely never arises when the data y are corrupted with "non degenerate" noise.

Theorem 1 Let y = y 0 + n where y 0 ∈ Ê m is fixed and n is an absolute continuous random vector, i.e., admitting a probability density w.r.t. the Lebesgue measure. Then, when Input: A, y, λ and support Q1 (Card[Q1] min(m, n))

Card[Q] min(m, n), the probability that x Q 0 < Card[Q] is equal to 0. Proof: Let k = Card[Q] and t Q be the minimizer of y-A Q t 2 over Ê k . t Q reads t Q = V Q y where matrix V Q = (A t Q A Q ) -1 A t Q is of size k × m, and x Q 0 = t Q 0 k. Denoting by v 1 , . . . , v k ∈ Ê m the row vectors of V Q , t Q 0 < k if
Step 1: Set j = 1.

Step 2: For i ∈ {1, . . . , n}, compute JQ j •i(λ).

Compute ℓ using [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF].

If J Q j •ℓ (λ) < JQ j (λ), Set Qj+1 = Qj • ℓ.
else, Terminate SBR.

End if.

Set j = j + 1 and go to Step 2.

Output: support Qj = SBR(Q1; λ) the inner product). Because A Q is full rank, V Q is full rank and then ∀i, v i = 0. Denoting by H ⊥ (v i )
the hyperplane of Ê m which is orthogonal to v i , we have

x Q 0 < k ⇐⇒ y ∈ k i=1 H ⊥ (v i ). (7)
Because the set i H ⊥ (v i ) has a Lebesgue measure equal to zero and the random vector y admits a probability density, the probability of event ( 7) is zero.

Theorem 1 implies that when dealing with real noisy data, it is almost sure that all active coefficients

x i are non-zero. Hence, each SBR iterate Q almost surely satisfies J (x Q ; λ) = J Q (λ).
In any case, SBR can be applied without restriction and the properties stated below (e.g., termination after a finite number of iterations) remain valid when an SBR iterate satisfies

x Q 0 < Card[Q].

D. Properties of SBR

Proposition 1 Under the assumptions of Theorem 1, each SBR iterate x Q is almost surely a local minimizer of J (x; λ). In particular, the SBR output satisfies this property.

Proof:

Let x = x Q be an SBR iterate. According to Theorem 1, the support S(x) = Q almost surely. Setting ε = min i∈Q |x i | > 0, it is easy to check that if x ′ ∈ Ê n satisfies x ′ -x < ε, then S(x ′ ) ⊇ S(x) = Q, thus x ′ 0 x 0 . Assume that x ′ satisfies x ′ -x < ε. • If S(x ′ ) = Q, then, by definition of x = x Q , we have E(x ′ ) E(x). Thus, J (x ′ ; λ) J (x; λ). June 16, 2011 DRAFT • Otherwise, J (x ′ ; λ) = E(x ′ ) + λ x ′ 0 E(x ′ ) + λ( x 0 + 1)
. By continuity of E, there exists

a neighborhood V(x) of x such that if x ′ ∈ V(x), |E(x ′ ) -E(x)| < λ. Thus, if x ′ ∈ V(x), x ′ -x < ε and S(x ′ ) ⊃ Q, then J (x ′ ; λ) > E(x) + λ x 0 = J (x; λ). Finally, if x ′ ∈ V(x) and x ′ -x < ε, then J (x ′ ; λ) J (x; λ).
Termination: Because SBR is a descent algorithm, a support Q cannot be explored twice and SBR terminates after a finite number of iterations. We emphasize that no stopping condition is needed unlike many algorithms which require to set a maximum number of iterations and/or a threshold on the squared error variation (CoSaMP, Subspace Pursuit, Iterative Hard Thresholding, Iterative Reweighted ℓ 1 ).

OLS as a special case: When λ = 0, 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][START_REF] Chen | Fast orthogonal least squares algorithm for efficient subset model selection[END_REF].

The removal operation never occurs because it yields an increase of the squared error

J Q (0) = E Q .
Empty solutions: We characterize the λ-values for which SBR yields an empty solution.

Remark 2 SBR(∅; λ) yields the empty set if and only if

λ λ max max i ( a i , y 2 / a i 2 ).
This result directly follows from checking that any insertion trial fails, i.e., ∀i, E {i} + λ E ∅ . It allows us to design an automatic procedure which sets a number of λ-values adaptively to the data in order to compute SBR solutions at different sparsity levels (see Section VI-D).

Relation between SBR and SMLR:

The main difference between both algorithms is that SMLR involves the inversion of a matrix of the form

A t Q A Q + αI Card[Q] whereas SBR computes the inverse of A t Q A Q .
In the case of SMLR, the term αI Card[Q] acts as a regularization on the amplitude values. It avoids instabilities when A Q is ill conditioned at the price of handling the additional hyperparameter α. On the contrary, instabilities may occur while using SBR. In the next section, we focus on this issue and propose a stable implementation.

IV. IMPLEMENTATION ISSUES

Given the current support Q, an SBR iteration consists in computing the squared error

E Q ′ for any replacement Q ′ = Q • i, leading to the computation of J Q ′ (λ) = E Q ′ + λCard[Q ′ ].
G Q A t Q A Q when Q is modified by one element. The Cholesky factorization takes the form G Q = L Q L t Q where L Q is a lower triangular matrix of size k = Card[Q]. Also, L Q is better conditioned than G Q ,
improving the stability of matrix inversion.

We now give the main updating equations. Full detailed derivation can be found in Appendix C.

A. Efficient strategy based on the Cholesky factorization

The replacement tests only rely on the current matrix L Q and do not require its update.

1) Single replacement tests: An insertion test Q ′ = Q ∪ {i} takes the form:

J Q ′ (λ) -J Q (λ) = λ - l t Q,i L -1 Q A t Q y -a t i y 2 a i 2 -l Q,i 2 (8) with l Q,i = L -1 Q A t Q a i .
This computation mainly requires a triangular system inversion (computation of l Q,i in O(k 2 ) elementary operations) up to the pre-computation of L -1 Q (A t Q y) at the beginning of the current SBR iteration.

According to [START_REF] Miller | Subset Selection in Regression[END_REF][START_REF] Reeves | An efficient implementation of the backward greedy algorithm for sparse signal reconstruction[END_REF], a removal test

Q ′ = Q\{i} reads J Q ′ (λ)-J Q (λ) = x Q (i) 2 /γ i -λ where x Q (i)
is the ith element in vector x Q and γ i is the diagonal element of G -1 Q corresponding to the position of a i in A Q . The overall removal tests mainly amount to the inversion of the triangular matrix L Q (in

O(k 3 ) operations) as the computation of γ i for all i and of G -1 Q A t Q y (i.e., the values of x Q (i)) from L -1 Q are both in O(k 2 ).
Note that insertion and removal tests can be easily done in parallel. In Matlab, this parallel implementation leads to a significant save of computation time due to the SIMD capabilities of Matlab.

2) Updating the Cholesky factorization:

The update of L Q can be easily done in the insertion case by adding the new column a i at the last position in A Q∪{i} . The new matrix L Q ′ is a 2 × 2 block matrix whose upper left block is L Q (see Appendix C). The removal case requires more care since a removal breaks the triangular structure of L Q . The update can be done by performing either a series of Givens planar rotations [START_REF] Haugland | A Bidirectional Greedy Heuristic for the Subspace Selection Problem[END_REF] or a positive rank 1 Cholesky update [START_REF] Gill | Methods for modifying matrix factorizations[END_REF]. We describe the latter strategy in Appendix C. The Cholesky factorization update is in O(k 2 ) in the insertion case and in O((k -

I) 2 )
in the removal case where I denotes the position of the column to be removed in A Q .

B. Reduced search

Additionally, we propose an acceleration of SBR yielding the same iterates with a reduced search.

We notice that a column removal Q ′ = Q\{i} yields an increase of the squared error and a decrease of the penalty equal to λ. Hence, the maximum decrease of J Q (λ) which can be expected is λ. The acceleration of SBR consists in testing insertions first. If any insertion leads to

J Q (λ) -J Q ′ (λ) > λ,
then removals are not worth being tested. Otherwise, the removals have to be tested as stated in Table I.

We have implemented this acceleration systematically.

C. Memory requirements and computation burden

The actual implementation may vary depending on the size and the structure of matrix A. We briefly describe the main possible implementations.

When the size of A is relatively small, the computation and storage of the Gram matrix A t A prior to any SBR iteration (storage of n 2 scalar elements) avoids to recompute the vectors A t Q a i which are needed when the insertion of a i into the active set is tested. The storage of the other quantities (mainly L Q ) that are being updated amounts to O(k 2 ) scalar elements and a replacement test costs O(k 2 ) elementary operations in average.

When A is larger, the storage of A t A is no longer possible, thus A t Q a i must be recomputed for any SBR iteration. This computation costs km elementary operations and now represents the most important part of an insertion test. When the dictionary has some specific structure, this limitation can be alleviated, enabling a fast implementation even for large n. For instance, if a large number of pairs of columns of A are orthogonal to each other, A t A can be stored as a sparse array. Also, finite impulse response deconvolution problems enable a fast implementation since A t A is then a Toeplitz matrix (save northwest and/or south-east submatrices, depending on the boundary conditions). The knowledge of the autocorrelation of the impulse response is sufficient to describe most of the Gram matrix.

All these variants have been implemented 1 . In the following, we analyze the behavior of SBR for two difficult problems involving highly correlated dictionaries: the deconvolution of a sparse signal with a Gaussian impulse response (Section V) and the joint detection of discontinuities at different orders in a signal (Section VI).

V. DECONVOLUTION OF A SPARSE SIGNAL WITH A GAUSSIAN IMPULSE RESPONSE

This is a typical problem for which SMLR was introduced [START_REF] Champagnat | Unsupervised deconvolution of sparse spike trains using stochastic approximation[END_REF]. It affords us to study the ability of SBR to perform an exact recovery in a simple noise-free case (separation of two Gaussian signals) and to test SBR in a noisy case (estimation of a larger number of Gaussians) and compare it with other algorithms.

For simulated problems, we denote by x ⋆ the exact sparse signal, the data reading y = Ax ⋆ + n. The dictionary columns are always normalized: a i 2 = 1. The signal to noise ratio (SNR) is defined by SNR = 10 log(P y /P n ), where P y = Ax ⋆ 2 /m is the average power of the noise-free data and P n is the variance of the noise process n.

1 Matlab codes provided by the authors can be downloaded at http://ieeexplore.org. In our Matlab implementation, the insertion and removal tests are done in parallel.

d = 20 0 0 2 ⋆ 2 ⋆ 2 ⋆ 2 ⋆ 2 ⋆ 2 ⋆ d = 13 0 1 3 4 5 2 ⋆ 2 ⋆ 2 ⋆ d = 6 0 1 1 3 5 6 8 2 ⋆

A. Dictionary and simulated data

The impulse response h is a Gaussian signal of standard deviation σ, sampled on a regular grid at integer locations. It is approximated by a finite impulse response of length 6σ by thresholding the smallest values, allowing for fast implementation even for large size problems (see subsection IV-C). The deconvolution problem leads to a Toeplitz matrix A whose columns are obtained by shifting the signal h. The dimension of A is chosen to have any Gaussian feature resulting from the convolution h * x ⋆ belonging to the observation window {1, . . . , m}. This implies that A is slightly undercomplete (m > n).

B. Separation of two close Gaussian features

We first analyze the ability of SBR to separate two Gaussian features ( x ⋆ 0 = 2) from noise-free data.

The centers of both Gaussian features lay at a relative distance d (expressed as a number of samples) and their weights x ⋆ i are set to 1. We analyze the SBR outputs for decreasing λ-values by computing their cardinality and testing whether they coincide with the true support S(x ⋆ ). Table II shows the results obtained for a problem of size 300 × 270 (σ = 5) with distances equal to d = 20, 13, and 6 samples.

It is noticeable that the exact recovery always occurs provided that λ is sufficiently small. This result remains true even for smaller distances (from d = 2). When the Gaussian features strongly overlap, i.e., for d 13, the size of the output support first increases while λ decreases, and then removals start to occur, enabling the exact recovery for lower λ's.

C. Behavior of SBR for noisy data

We consider a more realistic simulation in which the data are of larger size (m = 3000 samples) and noisy. The impulse response h is of size 301 (σ = 50) yielding a matrix A of size 3000 × 2700, and the SNR is set to 20 dB. composed of 17 spikes that are uniformly located in {1, . . . , n}. The non-zero amplitudes x ⋆ i are drawn according to an i.i.d. Laplacian distribution. Let us remark that the limit Bernoulli-Gaussian model is not a proper probabilistic model so that one cannot use it to design simulated data. We choose a Laplacian distribution since the non-zero amplitudes are more heterogeneous than with a Gaussian distribution with finite variance. On Fig. 1(b-d), we display the SBR results for three λ-values. For large λ's, only the main Gaussian features are found. When λ decreases, the smaller features are being recovered together with spurious features. Removals occur for λ 0.8 yielding approximations that are more accurate than those obtained with OLS and for the same cardinality (the residual y -Ax 2 is lower) while when λ > 0.8, the SBR output coincides with the OLS solution of same cardinality. Note that the theoretical value of λ obtained from (3) is equal to 0.3 yielding a support of cardinality 18. The residual is slightly lower than that obtained with λ = 0.5. The exact support of x ⋆ is never found because the data are noisy 1. SBR(λ = 0.5) is run first yielding a support of cardinality k sbr =18. Then, we run OLS(k sbr ), OMP(k sbr ), homotopy for ℓ1 regression [START_REF] Donoho | About SparseLab[END_REF], and IRℓ1(λ) [START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF]. The ℓ1 result is the homotopy iterate of cardinality k sbr yielding the least value of J (x; λ). and the neighboring columns of A are highly correlated. In such difficult case, one needs to perform a wider exploration of the discrete set {0, 1} n by introducing moves that are more complex than single replacements. Such extensions were already proposed in the case of SMLR. One can for instance shift an existing spike x i forwards of backwards [START_REF] Chi | Improved maximum-likelihood detection and estimation of Bernoulli-Gaussian processes[END_REF] or update a block of neighboring amplitudes jointly (e.g.,

x i and x i+1 ) [START_REF] Allain | Efficient binary reconstruction for non-destructive evaluation using gammagraphy[END_REF]. Various search strategies are also reported in [START_REF] Miller | Subset Selection in Regression[END_REF]Chapter 3].

D. Comparison of SBR with other sparse algorithms

We compared SBR with classical and recent sparse algorithms: OMP, OLS, CoSaMP [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF], Subspace Pursuit [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF], Iterative Hard Thresholding (IHT) [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF][START_REF] Blumensath | Normalized iterative hard thresholding: Guaranteed stability and performance[END_REF], ℓ 1 regression [START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF] and Iterative Reweighted ℓ 1 (IRℓ 1 ) [START_REF] Candès | Enhancing sparsity by reweighted ℓ1 minimization[END_REF][START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF]. A general trend is that thresholding algorithms perform poorly when the dictionary columns are strongly correlated. CoSaMP and Subspace Pursuit yield the worst results: they stop after a very few iterations as the squared error increases from one iteration to the next. On the contrary, IHT guarantees that the squared error decreases but the convergence is very slow and the results remain poor in comparison with SBR. In the simulation of Fig. 1(c), SBR performs 12 iterations (only insertions) leading to a support of cardinality 12. Meanwhile, the number of iterations of IHT before convergence is huge: both versions of IHT presented in [START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF] require at least 10,000 iterations to converge, leading to an overall computation time (22 and 384 seconds) that is much larger than the SBR computation time (3 seconds). Fig. 2 is a synthetic view of the performance of SBR, OLS, OMP, ℓ 1 regression, and IRℓ 1 for a given sparsity level λ. The computation time and the value of J (x; λ) are shown on the horizontal and vertical axes, respectively. This enables us to define several categories of algorithms depending on their locations on the 2D plane: the outputs of fast algorithms (OMP and ℓ 1 ) lay in the upper left region whereas slower but more efficient algorithms (OLS, SBR, and IRℓ 1 ) yield points laying in the lower right region. We chose not to represent the outputs of thresholding algorithms since they yield poorer performance, i.e., points located either in the upper right (IHT) or upper left (CoSaMP, Subspace Pursuit) regions. In details, we observed that ℓ 1 regression tends to overestimate the support cardinality and to place several spikes at very close locations. We used Donoho's homotopy implementation [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] and found that it requires many iterations: homotopy runs during 200 iterations before reaching a support of cardinality 18 when processing the data of Fig. 1 (we recall that homotopy starts from the empty set and performs a single support replacement per iteration). The performance of ℓ 1 regression fluctuates around that of OMP depending on the trials and the sparsity level. Regarding IRℓ 1 , we used the Adaptive LASSO implementation from Zou [START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF] since it is dedicated to the minimization of J (x; λ). We stopped the algorithm when two successive ℓ 1 iterates share the same support. For the simulation of Fig. 1, IRℓ 1 and SBR yield comparable results in that one algorithm does not outperform the other for all λ values, but IRℓ 1 generally performs slightly better (Fig. 2). We designed other simulations in which the nonzero weights x ⋆ i are spread over a wider interval. In this case, SBR most often yields the best approximations. Fig. 2 is representative of the empirical results obtained while performing many trials. Obviously, the figure may significantly change depending on several factors among which the λ-value and the tuning parameters of IRℓ 1 . The goal is definitely not to conclude that an algorithm always outperforms the others but rather to sketch a classification of groups of algorithms according to the trade-off between accuracy and computation time.

VI. JOINT DETECTION OF DISCONTINUITIES AT DIFFERENT ORDERS IN A SIGNAL

We now consider another challenging problem: the joint detection of discontinuities at different orders in a signal [START_REF] Smith | Nonparametric regression using Bayesian variable selection[END_REF][START_REF] Vetterli | Sampling signals with finite rate of innovation[END_REF]. We process both simulated and real data and compare the performance of SBR with respect to OMP, Bayesian OMP (BOMP) which is an OMP based forward-backward algorithm [START_REF] Herzet | Bayesian pursuit algorithms[END_REF], OLS, ℓ 1 regression [START_REF] Donoho | Fast solution of ℓ1-norm minimization problems when the solution may be sparse[END_REF], and IRℓ 1 [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]. Firstly, we formulate the detection of discontinuities at a single order as a spline approximation problem. Then, we take advantage of this formulation to introduce the joint detection problem.

Signal a 0 i i -1 i i + 1 Signal a 1 i i -1 i i + 1 Signal a 2 i i -1 i i + 1
Fig. 3. Signals a p i related to the pth order discontinuities at location i. a 0 i is the Heaviside step function, a 1 i is the ramp function, and a 2 i is the one-sided quadratic function. Each signal is equal to 1 at location i and its support is equal to {i, . . . , m}.

A. Approximation of a spline of degree p

Following [START_REF] Smith | Nonparametric regression using Bayesian variable selection[END_REF], we introduce the dictionary A p of size m × (mp) formed of shifted versions of the one-sided power function k → [max(k, 0)] p for all possible shifts (see Fig. 3) and we address the sparse approximation of y by the piecewise polynomial A p x p (actually, we impose as initial condition that the spline function is equal to 0 for k 0). It consists in the detection of the discontinuity locations (also referred to as knots in the spline approximation literature) and the estimation of their amplitudes:

x p i codes for the amplitude of a jump at location i (p = 0), the change of slope at location i (p = 1), etc. Here, the notion of sparsity is related to the number of discontinuity locations.

B. Piecewise polynomial approximation

We formulate the joint detection of discontinuities of orders p = 0, . . . , P by appending the elementary dictionaries A p in a global dictionary A = [A 0 , . . . , A P ]. The product Ax yields a sum of piecewise polynomials of degree lower than P with a limited number of pieces. The dictionary A is overcomplete since it is of size m × s, with s = (P + 1)(m -P/2) > m for P 1. Moreover, any column a p i of A p overlaps all other columns a q j because their respective supports are the intervals {i, . . . , m} and {j, . . . , m}. The discontinuity detection problem is difficult as most algorithms are very likely to position wrong discontinuities in their first iterations. For example, when approximating a signal with two discontinuities at distinct locations i and j, greedy algorithms start to position a first (wrong) discontinuity in between i and j, and forward greedy algorithms cannot remove it.

C. Adaptation of SBR

The above defined dictionary does not satisfy the unique representation property. Indeed, it is easy to check that the difference between two discrete ramps at locations i and i + 1 yields the discrete Heaviside function at location i: a 1 ia 1 i+1 = a 0 i . We thus need to slightly modify SBR in order to ensure that only full rank matrices A Q are explored. The modification is based on the following proposition which gives a sufficient condition for full rankness of A Q .

Proposition 2 Let n i denote the number of columns a p i , p ∈ {0, . . . , P } which are active for sample i. Let us define the binary condition C(i):

• if n i = 0, C(i) 1; • if n i 1, C(i) n i+j = 0, j = 1, . . . , n i -1 If Q satisfies ∀i, C(i) = 1, then A Q is full rank.
Proposition 2 is proved in Appendix D. Basically, it states that we can allow several discontinuities to be active at the same location i, but then, the next samples i+1, . . . , i+n i -1 must not host any discontinuity.

This condition ensures that there are at most n i discontinuities in the interval {i, . . . , i + n i -1} of length n i . The SBR adaptation consists in testing an insertion only when the new support

Q ′ = Q ∪ {(i, p)}
satisfies the above condition.

D. Numerical simulations

We first set P = 1 leading to the piecewise affine approximation problem. The noise-free data y = Ax ⋆ of Fig. 4(a) are of size m = 1000 with x ⋆ 0 = 18 discontinuities. According to Remark 2, we compute the value λ max above which the SBR output is the empty set, and we run SBR with λ j = λ max 10 -j/2

for j = 0, . . . , 20. For the least λ-value, SBR yields an exact recovery (see Fig. 4(a)). For comparison purpose, we also run 27 iterations of OMP and OLS. The "ℓ 2 -ℓ 0 " curves represented on Fig. 4(b) express the squared residual y -Ax 2 versus the cardinality x 0 for each algorithm (we plot the first 27 iterates of OMP and OLS and for all j, we plot the output of SBR(λ j ) after full convergence of SBR).

Whatever the cardinality, SBR yields the least residual. For noisy data, the "ℓ 2 -ℓ 0 " curve corresponding to SBR still lays below the OMP and OLS curves for most cardinalities. In the next paragraph, we also consider the Bayesian OMP, ℓ 1 regression, and IRℓ 1 algorithms for further comparisons.

E. AFM data processing

In Atomic Force Microscopy (AFM), a force curve measures the interatomic forces exerting between a probe associated to a cantilever and a nano-object. Specifically, the recorded signal z → y(z) shows the force evolution versus the probe-sample distance z, expressed in nanometers. Researching discontinuities (location, order, and amplitude) in a force curve is a challenging task because they are used to provide a precise characterization of the physico-chemical properties of the nano-object (topography, energy of adhesion, etc.) [START_REF] Butt | Force measurements with the atomic force microscope: Technique, interpretation and applications[END_REF].

The data displayed on Fig. 5(a) are related to a bacterial cell Shewanella putrefaciens laying in aqueous solution, interacting with the tip of the AFM probe [START_REF] Gaboriaud | Spatially resolved force spectroscopy of bacterial surfaces using force-volume imaging[END_REF]. A retraction force curve is recorded by positioning the tip in contact with the bacterial cell, and then gradually retracting the tip from the sample until it We search for the discontinuities of orders 0, 1, and 2. Similar to the processing of simulated data, we run SBR with 14 λ-values and we run OLS and OMP until iteration 41. For each algorithm, we plot the "ℓ 2 -ℓ 0 " curve and the curve displaying the time of reconstruction versus the cardinality (Figs. 5(b,c)). These figures show that the performance of SBR is at least equal and sometimes better than that of OLS.

Both algorithms yield results that are far more accurate than OMP at the price of a larger computation time.

Fig. 6 displays the approximations yielded by the three algorithms together with the BOMP, ℓ 1 , and IRℓ 1 approximations. For the largest value λ 1 , SBR runs during 6 iterations (4 insertions and 2 removals) yielding a support of cardinality 2. SBR performs better than other algorithms (Figs. 6(a-f)). Although IRℓ 1 yields the most accurate approximation, it relies on 4 dictionary columns leading to a larger value of J (x; λ 1 ). We observed the same behavior for the lowest value λ 2 (subfigures (g-l)). Again, SBR yields the least value of J (x; λ 2 ) among all algorithms. Moreover, SBR provides a very precise localization of both first order discontinuities (subfigure (a)) which are crucial information for the physical interpretation of the data. On the contrary, all other algorithms fail for the highest sparsity level, and some do not even succeed for the lowest. Specifically, OLS accurately locates both first order discontinuities when 5 iterations have been performed (the desired discontinuities are the first and the last ones among the 5) while OMP fails even after 5 iterations. LASSO and BOMP yield very poor approximations for the highest sparsity level and approximations with many dictionary columns for the lowest sparsity level. In terms of value of the cost function J (x; λ), BOMP and LASSO fluctuate around OMP but they are far outperformed by OLS, SBR, and IRℓ 1 .

VII. CONCLUSION

A. Discussion

We performed comparisons for two problems involving highly correlated dictionary columns. SBR is at least as accurate as OLS and sometimes more accurate, with a slightly larger cost of computation.

We also considered sparse algorithms that are slower than OLS. SBR was found to be very competitive in terms of trade-off between accuracy and computation time. Although OLS based forward-backward algorithms yield a relatively large computational cost per iteration, we have noticed that for correlated dictionaries, the number of SBR iterations (i.e., of elementary modifications of the support) is much lower than the number of support modifications performed by several other algorithms. Typically, IHT and IRℓ 1 can often be more expensive than SBR. Additionally, SBR terminates within a finite number of iterations, thus it does not require to tune any empirical stopping parameter. The limitation of SBR in terms of speed arises when the dictionary A is unstructured and the size of A is too large to store A t A.

The inner products a t i a j must then be recomputed for each iteration, which is relatively burdensome. In the recent literature, it is often acknowledged that the cost function J (x; λ) has a large number of local minimizers therefore discouraging its direct optimization [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]. Many authors thus choose to minimize an approximate cost function in which the ℓ 0 norm |x i | 0 is replaced with a nonconvex continuous function ϕ(x i ). However, when the range of values of the (expected) nonzero amplitudes x i = 0 is wide, it is difficult to find a good approximation ϕ(x i ) of |x i | 0 for all x i . Selecting an appropriate ϕ function generally relies on the introduction of a degree of freedom whose tuning is not obvious [START_REF] Candès | Enhancing sparsity by reweighted ℓ1 minimization[END_REF][START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed ℓ 0 norm[END_REF]. For instance, the IRℓ 1 algorithm can be interpreted as an approximate ℓ 2 -ℓ 0 minimization method where the ℓ 0 norm is replaced with ϕ(x i ; ε) = 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]. The parameter ε controls the "degree of nonconvexity" of the surrogate function ϕ ( 2 ). Although J (x; λ) has a large number of local minima, we have found that SBR is often as accurate as algorithms based on the nonconvex approximation of J . Moreover, SBR is simple to use. The good behavior of SBR is somehow related to the result of Proposition 1 which states that any SBR iterate is almost surely a local minimizer of J . We conclude that SBR is actually capable to "skip" local minima with a large cost J (x; λ).

B. Perspectives

In the proposed approach, the main difficulty relies in the choice of the λ-value. If a specific cardinality or approximation residual is desired, one can resort to a trial and error procedure in which a number of λ-values are tried until the desired approximation level is found. In [START_REF] Duan | A continuation approach to estimate a solution path of mixed L2-L0 minimization problems[END_REF], we sketched a continuation version in which a series of SBR solutions are computed for decreasing levels of sparsity λ, and the λ-values are recursively computed. This continuation version is showing promising results and will be the subject of a future extended contribution. A similar perspective was actually proposed by Zhang to generalize his FoBa algorithm in a path-following algorithm (see the discussion section in [START_REF] Zhang | Adaptive forward-backward greedy algorithm for learning sparse representations[END_REF]).

Another important perspective is to investigate whether SBR can guarantee exact recovery in the noise-free case under some conditions on matrix A and on the unknown sparse signal x ⋆ . According to Remark 1, we will study the behavior of SBR when λ → 0. In the simulations done in Sections V and VI, we observed that SBR is able to perform exact recoveries provided that λ is sufficiently small. This promising result is a first step towards a more general theoretical study.

APPENDIX A DETAILED DEVELOPMENT OF LIMIT BG SIGNAL RESTORATION

Consider the Bernoulli-Gaussian model x = (q, r) introduced in Section II-B and the joint MAP formulation (3) involving the cost function L(q, r). Given q, let us split r into two subvectors u and t indexed by the null and non-null entries of q, respectively. Since r 2 = t 2 and A∆ q r = A Q t do not depend on u, we have min u L(q, t, u) = L(q, t, 0). Thus, the joint MAP estimation problem reduces to the minimization of L(q, t, 0) w.r.t. (q, t). In the limit case σ 2

x → ∞, this problem rereads: min q,t {L(q, t, 0) = y -A Q t 2 + λ q 0 }.

The equivalence between ( 9) and (4) directly follows from the change of variable x = {q, t} where q and t are the support and non-zero amplitudes of x.

APPENDIX B PROOF OF REMARK 1

The proof of the result stated in Remark 1 is based on the two following lemmas.

Lemma 1 For λ > 0, any minimizer of J (x; λ) takes the form x Q with Card[Q] min(m, n).

Proof of lemma: According to the URP assumption, any min(m, n) columns of A yield an unconstrained minimizer of y -Ax 2 . Let x LS be such minimizer, with x LS 0 min(m, n), and let u be a minimizer of J (x; λ). J (u; λ) J (x LS ; λ) implies that u 0 x LS 0 + (E(x LS ) -E(u))/λ

x LS 0 min(m, n).

We denote by Q the support of u. The related least-square solution x Q obviously satisfies E(x Q ) E(u) and x Q 0 Card[Q] = u 0 , thus J (x Q ; λ) J (u; λ). Since u is a minimizer of J (x; λ), we have J (x Q ; λ) = J (u; λ) hence E(x Q ) = E(u). Because of the URP assumption, the least-squares minimizer over Q is unique, thus u = x Q .

Lemma 2 There exists λ min > 0 such that for 0 < λ λ min , the minimizers of J (x; λ) are unconstrained minimizers of y -Ax 2 .

Proof of lemma: When λ tends towards 0, we have for all Q, J (x

Q ; λ) = E Q + λ x Q 0 → E Q .
In particular, J (x QLS ; λ) → E QLS with x QLS an unconstrained minimizer of y -Ax 2 yielded by a subset Q LS of cardinality min(m, n). Because the number of possible subsets Q is finite and for all Q, E Q E QLS , there exists λ min > 0 such that for 0 < λ λ min , the subsets Q ⋆ minimizing J (x Q ; λ)

satisfy E Q ⋆ = E QLS . Consequently, the minimizers of J (x; λ) are unconstrained least-squares solutions according to Lemma 1.

Proof of Remark 1:

The proof directly follows from the application of Lemma 2. We denote by X λ the set of minimizers of J (x; λ).

In the undercomplete case, there is a unique unconstrained least-square minimizer x LS . Thus, X λ = {x LS } for 0 λ λ min .

In the overcomplete case, we denote by X ⋆ the set of sparsest solutions to y = Ax. To show that X λ = X ⋆ for 0 < λ λ min , we consider x ∈ X ⋆ and x ′ ∈ X λ . According to Lemma 2, x ′ satisfies y = Ax ′ , then J (x ′ ; λ) = λ x ′ 0 . By definition of X ⋆ , we have y = Ax and J (x; λ) = λ x 0 J (x ′ ; λ).

Because x ′ ∈ X λ is a minimizer of J , we deduce that x ′ 0 = x 0 , then x ′ ∈ X ⋆ and x ∈ X λ . We have proved that X λ = X ⋆ for 0 < λ λ min .

APPENDIX C UPDATE OF THE CHOLESKY FACTORIZATION

At each SBR iteration, n linear systems of the form t Q G -1 Q A t Q y must be solved, the corresponding squared errors reading

E Q = y -A Q t Q 2 = y 2 -y t A Q t Q .
Using the Cholesky factorization

G Q = L Q L t Q , t Q rereads t Q = L -t Q L -1 Q A t Q y, thus E Q = y 2 -L -1 Q A t Q y 2 . ( 10
)
Insertion of a new column after the existing columns: Including a new column leads to

A Q ′ = [A Q , a i ].
Thus, the new Gram matrix reads as a 2 × 2 block matrix:

G Q ′ =   G Q A t Q a i (A t Q a i ) t a i 2  
and the Cholesky factor of G Q ′ can be straightforwardly updated:

L Q ′ =   L Q 0 l t Q,i a i 2 -l Q,i 2   (11)
with l Q,i = L -1 Q A t Q a i . The update (8) of J Q (λ) = E Q + λCard[Q] directly follows from ( 10) and [START_REF] Blumensath | Normalized iterative hard thresholding: Guaranteed stability and performance[END_REF]. Removal of an arbitrary column: When removing a column a i , updating L Q remains possible although more complex. This idea was developed by Ge et al. [START_REF] Ge | Enhanced sampling schemes for MCMC based blind Bernoulli-Gaussian deconvolution[END_REF] who update the Cholesky factorization of matrix

G -1
Q . We adapt it to the direct (simpler) factorization of G Q . Let I be the position of a i in A Q (with 1 I Card[Q]). L Q can be written in a block matrix form:

L Q =      Λ 0 0 b t d 0 C e F      ( 12
)
where the lowercase characters refer to the scalar (d) and vector quantities (b, e) appearing in the Ith row and in the Ith column. The computation of G Q = L Q L t Q and the removal of the Ith row and the Ith column in G Q lead to

G Q ′ =   Λ 0 C F     Λ t C t 0 F t   +   0 e   0 e t .
By identification with G Q ′ = L Q ′ L t Q ′ and because the Cholesky factorization is unique, L Q ′ necessarily reads:

L Q ′ =   Λ 0 C X   , ( 13
)
where X is a lower triangular matrix satisfying XX t = F F t + ee t . The problem of computing X from F and e is classical; it is known as a positive rank 1 Cholesky update and there exists a stable algorithm in O(f 2 ) operations, where f = Card[Q] -I is the size of F [START_REF] Gill | Methods for modifying matrix factorizations[END_REF].

  Our implementation is inspired by the fast implementation of the homotopy algorithm for ℓ 1 regression [3, 34]. It consists in maintaining the Cholesky factorization of the Gram matrix

5 Fig. 1 .

 51 Fig. 1. Gaussian deconvolution results. Problem of size 3000 × 2700 (σ = 50). (a) Generated data, with 17 Gaussian features and with SNR = 20 dB. The exact locations x ⋆ are labeled o. (b,c,d) SBR outputs and data approximations with empirical settings of λ. The estimated amplitudes x are shown with vertical spikes. The SBR outputs (supports) are of size 5, 12, and 18, respectively. The computation time always remains below 3 seconds (Matlab implementation).

Fig. 2 .

 2 Fig. 2. Comparison of sparse algorithms in terms of trade-off between accuracy (J (x; λ)) and CPU time for the deconvolution problem of Fig. 1. SBR(λ = 0.5) is run first yielding a support of cardinality k sbr =18. Then, we run OLS(k sbr ), OMP(k sbr ),

Fig. 4 .

 4 Fig. 4. Joint detection of discontinuities of orders 0 and 1. The dictionary is of size 1000 × 1999 and the data signal includes 18 discontinuities. The true and estimated discontinuity locations are represented with unfilled black and filled gray labels. The shape of the labels (circular or triangular) indicates the discontinuity order. The dashed gray and solid black curves represent the data signal y and its approximation Ax for the least λ-value. (a) Approximation from noise-free data. The recovery is exact. (b) "ℓ2-ℓ0" curves showing the squared residual versus the cardinality for the SBR, OLS, and OMP solutions. (c,d) Similar results for noisy data (SNR = 20 dB).

Fig. 5 .

 5 Fig. 5. Joint detection of discontinuities of orders 0, 1, and 2 (problem of size 2167 × 6498). (a) Experimental AFM data showing the force evolution versus the probe-sample distance z. (b) Squared residual versus cardinality for the SBR, OLS, and OMP solutions. (c) Time of reconstruction versus cardinality.

7 Fig. 6 .

 76 Fig.6. AFM data processing: joint detection of discontinuities at orders 0, 1, and 2. The estimated discontinuities x are represented with vertical spikes and with a label indicating the discontinuity order. (a) SBR output of cardinality 2: 4 insertions and 2 removals have been done (λ1 = 120). (b-f) OLS and OMP outputs after 2 iterations, BOMP and IRℓ1[START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF] outputs for λ = λ1, homotopy iterate (LASSO) leading to the minimal value of J (x; λ1). (g-l) Same simulation with a lower λ-value (λ2 = 8.5). The SBR output is of cardinality 5 (7 insertions and 2 removals).

TABLE I SBR

 I ALGORITHM. BY DEFAULT, Q1 = ∅.

TABLE II SEPARATION

 II OF TWO GAUSSIAN FEATURES FROM NOISE-FREE DATA WITH SBR. d STANDS FOR THE DISTANCE BETWEEN THE GAUSSIAN FEATURES. WE DISPLAY THE SIZE OF THE SUPPORT OBTAINED FOR A SEQUENCE OF DECREASING λ-VALUES λ0 > λ1 > . . . > λ7. THE LABEL ⋆ INDICATES AN EXACT RECOVERY FOR A SUPPORT OF CARDINALITY 2.

	λ	λ0	λ1	λ2	λ3	λ4	λ5	λ6	λ λ7

June 16, 2011 DRAFT

In the comparisons with SBR, we set ε = 0 following[START_REF] Zou | The adaptive Lasso and its oracle properties[END_REF].June 16, 2011 DRAFT

This work was carried out in part while C. Soussen was visiting

APPENDIX D

PROOF OF PROPOSITION 2 Let us first introduce some notations specific to the piecewise polynomial dictionary problem. Consider a subset Q of columns a p i and let i -= min{i | n i > 0} denote the lowest location of an active entry (we recall that n i denotes the number of active columns for sample i). Up to a reordering of the columns

where A i -gathers the n i -active columns a p i such that i = i - and A i -gathers the remaining active columns (with i > i -). The following lemma is a key element to prove Proposition 2.

Proof: Let I = n i -denote the number of discontinuities at location i -and let 0 p 1 < p 2 < . . . < p I denote their orders, sorted in the ascending order. Suppose that there exist two families of scalars

Let us show that all µ-values are then equal to 0.

Rewriting the first I nonzero equations in this system and because Q satisfies the condition of Proposition 2, we have, for all k ∈ {i -, . . . , i -+ I -1}, I j=1 µ pj i -(k + i --1) pj = 0. Hence, the polynomial F (X) = I j=1 µ pj i -X pj has I positive roots. Because any non-zero polynomial formed of I monomials of different degree has at most I -1 positive roots [47, p. 76], F is the zero polynomial, thus all scalars µ pj iare 0. We deduce from [START_REF] Couvreur | On the optimallity of the backward greedy algorithm for the subset selection problem[END_REF] and from the full rankness of A i -that µ p i = 0 for all (i, p). We have shown that the column vectors of A Q are linearly independent, i.e., that A Q is full rank.

The proof of Proposition 2 directly results from the recursive application of Lemma 3. Starting from the empty set, all the indices, sorted by decreasing order, are successively included.

ACKNOWLEDGMENT

The authors would like to thank Dr. Grégory Francius from LCPME (UMR CNRS 7564, Nancy, France) for providing them with real AFM data.