Improved iterative shrinkage-thresholding for sparse signal recovery via Laplace mixtures models

In this paper, we propose a new method for support detection and estimation of sparse and approximately sparse signals from compressed measurements. Using a double Laplace mixture model as the parametric representation of the signal coefficients, the problem is formulated as a weighted ℓ1 minimization. Then, we introduce a new family of iterative shrinkage-thresholding algorithms based on double Laplace mixture models. They preserve the computational simplicity of classical ones and improve iterative estimation by incorporating soft support detection. In particular, at each iteration, by learning the components that are likely to be nonzero from the current MAP signal estimate, the shrinkage-thresholding step is adaptively tuned and optimized. Unlike other adaptive methods, we are able to prove, under suitable conditions, the convergence of the proposed methods to a local minimum of the weighted ℓ1 minimization. Moreover, we also provide an upper bound on the reconstruction error. Finally, we show through numerical experiments that the proposed methods outperform classical shrinkage-thresholding in terms of rate of convergence, accuracy, and of sparsity-undersampling trade-off.


Introduction
In this paper, we consider the standard compressed sensing (CS) setting [1], where we are interested in recovering high-dimensional signals x ∈ R n from few linear measurements y = Ax + η, where A ∈ R m×n , m n, and η is a Gaussian i.i.d. noise. The problem is underdetermined and has infinitely many solutions. However, much interest has been focused on finding the sparsest solution, i.e., the one with the smallest number of nonzero components [2]. This involves the minimization of 0 pseudonorm [3], which is NP-hard.
A practical alternative is to use the 1 regularization, leading to the basis pursuit (BP, [4]) problem in the absence of noise, or the least absolute shrinkage and selection operator (Lasso, [5]) in the presence of noise. They can be efficiently solved by iterative shrinkagethresholding algorithms (ISTA, [6][7][8]) that are generally *Correspondence: chiara.ravazzi@ieiit.cnr.it 1 National Research Council of Italy, IEIIT-CNR, c/o Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy Full list of author information is available at the end of the article first-order methods followed by a shrinkage-thresholding step. Due to its implementation simplicity and suitability for high-dimensional problems, a large effort has been spent to improve their speed of convergence [9][10][11][12], asymptotic performance in the large system limit [13,14], and ease of use [15].
In a Bayesian framework, 1 minimization is equivalent to a maximum a posteriori (MAP) estimate [16] modeling the signal coefficients using a Laplace prior, in the sense that we need to solve the same optimization problem. Although the Laplace probability density function does not provide a relevant generative model for sparse or compressible signals [17], the non-differentiability at zero of the cost function leads to select a sparse solution, providing empirical success of 1 regularization.
However, 1 minimization alone does not fully exploit signal sparsity. In fact, in some cases, a support estimate [18] could be employed to reduce the number of measurements needed for good reconstruction via BP or Lasso, e.g., by combining support detection with weighted or truncated 1 minimization [19]. The idea of combining support information and signal estimation has appeared in CS literature with several assumptions [20][21][22][23][24][25][26][27]. For example, in [25], the authors employ as prior information an estimate T of the support of the signal and propose a truncated 1 minimization problem. Another piece of literature [28] considers a weighted 1 minimization with weights w i = −log p i where p i is the probability that x i = 0.
In this paper, we propose an iterative soft support detection and estimation method for CS. It is worth remarking that in our setting prior information on the support T or p i is not available. The fundamental idea is to combine the good geometric properties of the 1 cost function associated to the Laplacian prior with a good generative model for sparse and compressible vectors [17]. For this purpose, we use a Laplace mixture model as the parametric representation of the prior distribution of the signal coefficients. Because of the partial symmetry of the signal sparsity, we know that each coefficient should have one out of only two distributions: a Laplace with small variance with high probability and a Laplace with large variance with low probability. We show empirically that this model fits better with the distribution of the Haar wavelet coefficients in test images. Then, we cast the estimation problem as a weighted 1 minimization method that incorporates the parametric representation of the signal.
We show that the proposed framework is able to improve a number of existing methods based on shrinkage-thresholding: by estimating the distribution of the components that are likely to be nonzero from signal estimates at each iteration (support detection), the shrinkage-thresholding step is tuned and optimized, thereby yielding better estimation. As opposed to other adaptive methods [10], we are able to prove, under suitable conditions, the convergence of the proposed tuned method. Moreover, we derive an upper bound on the reconstruction error. We apply this method to several algorithms, showing by numerical simulation that it improves recovery in terms of both speed of convergence and sparsity-undersampling trade-off, while preserving the implementation simplicity.
Compared to the literature on reconstruction methods that combine iterative support detection and weighted 1 minimization, the identification of the support is not nested or incremental over time as in [29][30][31]. Moreover, the choice of weights in 1 minimization is based on the Bayesian rules and a probabilistic model and not on greedy rules as in [19,32]. This feature also marks the difference with respect to reweighted 1 / 2 minimization, where the weights are chosen with the aim of approximating the τ norm with τ ∈ (0, 1] .

Outline
The paper is organized as follows. In Section 2, the basic CS theory and the classical methods based on 1 minimization for sparse recovery are reviewed. The proposed parametric model for sparse or highly compressible signals is described in Section 3 and compared with the related literature. In Section 4, the estimation problem based on Laplace mixture models is introduced and recast as a weighted 1 minimization problem. Then, in Section 5, the proposed approach is used to improve a number of existing methods based on shrinkagethresholding. Numerical experiments are presented in Section 6 and some concluding remarks (Section 7) complete the paper. The theoretical results are rigorously proved in Appendices 1, 2, 3, and 4.

Notation
We conclude this introduction with some notation. We denote column vectors with small letters and matrices with capital letters. If x ∈ R n , we denote its jth element with x j and, given S ∈[ n] := {1, . . . , n}, by x| S the subvector of x corresponding to the indexes in S. The support set of x is defined by supp(x) = {i ∈[ n] : x i = 0}, and we use x 0 = |supp(x)|. Finally, the symbol x with no subscript is to be understood as the Euclidean norm of the vector x. We denote as r(x) the nonincreasing rearrangement of x We denote with s = {x ∈ R n : |supp(x)| ≤ s} and define It should be checked that Given a matrix A, A T denotes its transpose.

Sparse signal recovery from compressed measurements
Compressive sensing aims to recover a sparse signal x ∈ R n from m ≤ n random projections of the form where y ∈ R m is the observation vector, A ∈ R m×n is the measurement matrix, and η is an additive noise. For example, in the transform domain compressive sig-nal reconstruction [1], A = , where ∈ R n×n is the sparsifying basis (i.e., multiplying by corresponds to performing inverse transform), the entries of x is the transform coefficient vector that has k nonzero entries, and ∈ R m×n is the sensing matrix, whose rows are incoherent with the columns of .
Conventional reconstruction methods involve 1 regularization [4]. In particular, it has been shown that, in the absence of noise and under suitable assumptions on the matrix A, the basis pursuit problem (BP) can exactly recover a k-sparse signal (i.e., with a number of nonzero coefficients not larger than k) from m = O(k log(n/k)) measurements with high probability [33,34]. In the presence of noise, one of the most popular convex relaxation methods is the least absolute shrinkage and selection operator (Lasso, [5]), which requires to solve the following unconstrained problem where λ is a positive regularization parameter. Even when the vector is not exactly sparse, under compressibility assumptions of the signals to be recovered, the 1 regularization provides estimates with a controlled error [33]. More formally, we recall the following definition [17].
If x ∈ R n is not exactly sparse but compressible, the support is intended as the set of significant components supp( k (x)).
One drawback of classical 1 minimization is that it fails to penalize the coefficients in different ways. In this paper, we propose a new family of methods that incorporate two tasks: iterative support detection and signal recovery.

A Bayesian view
In a Bayesian framework, if the noise in (1) is white Gaussian (we suppose for simplicity unitary standard deviation), BP and Lasso may be interpreted as a Bayesian MAP estimate [16]. In fact, imposing the 1 norm as penalty in the cost function is equivalent to modeling the signal coefficients x i as independent and identically distributed as a Laplace distribution, namely where, using the Bayes rule, f x|y (x|y) is given by in case of BP, and for Lasso by and Z is a normalization factor so that f * x|y (x|y)dx = 1. Despite the good geometric properties of the 1 cost function associated to such prior, that allow to select a sparse solution, the Laplace prior does not provide a relevant generative model for sparse or compressible signals [17]. In fact, if x n ∈ R n is distributed as i.i.d. with respect to Laplace distribution with scale parameter λ, then for any sequence k n such that lim n→∞ k n /n = κ ∈[ 0, 1], it holds almost surely Therefore, the vectors generated from i.i.d. Laplace distribution are not compressible since we cannot have both κ and ε small at the same time.
Moreover, for a large class of real signals that have highly non-Gaussian statistics, the Laplace model does not provide a good fit to the empirical probability density function. We show this with a simple experiment (experiment 1) Fig. 1. We calculate a single vertical wavelet subband coefficients of several real images of size 256 × 256 pixels, and we compute for each image the best fitting of the double Laplace density function obtained by maximizing the likelihood of the data under that assumption. In Fig. 2, the empirical density function of the Haar wavelet coefficients using 256 bins is shown in the log domain (solid line) and the dashed line corresponds to the best fitting instance of the Laplace density function. It should be noticed that the Laplace density captures peaks at zero but is less accurate along the tails.
In [35], Lasso is proved to provide a robust estimation that is invariant to the signal prior. In sharp contrast, the Bayesian Lasso is able to provide an estimator with minimum mean squared error by incorporating the signal model in the estimation problem [14,36], but the assumption that the signal prior is known in advance is not reasonable in most practical cases. Hence, it becomes crucial to incorporate in the recovery procedure new tools for adaptively learning sparsity models. Other models have been proposed for compressible signals [37][38][39] using more accurate probability density functions than double Laplace distribution. However, two issues generally appear when an accurate but complex signal prior is used: (1) it can be hard to estimate the model parameters and (2) the optimal estimators may not have simple closed form solution and their computation may require high computational work [40]. In fact, although the double Laplace prior is not the most accurate model, this is an especially convenient assumption since the MAP estimator has a simple and closed form [41].
Our goal is to use a compressible distribution as parametric representation of the signal coefficients, able to combine support detection and estimation, and to

Proposed approach: two-component Laplace mixture for support detection
We consider a two-state mixture model as a prior that describes our knowledge about the sparsity of the solution to (1). Because of the partial symmetry of the signal sparsity, we consider the case in which x is a random variable with components of the form where u i are identically and independently distributed (i.i.d.) as Laplace(0, α), v i are i.i.d. according to Laplace(0, β) and z i are i.i.d. Bernoulli random variables with probability mass function f (z i = 1) = 1 − p, with p < < 1/2, α ≈ 0, and β > > 0, in order to ensure that we have few large coefficients. We thus consider the conditional distribution of the data: let = (α, β) where in presence of noise, and Z is a normalization factor so that f (x|y; )dx = 1. This mixture model is completely described by three parameters: the sparsity ratio p < < 1/2, α that is expected to be small, and β > α if the signal is sparse. It should be noticed that vectors generated from this distribution are typically compressible according to Definition 1.
Proposition 1 Let x n ∈ R n be i.i.d. with respect to ( 4). Then, for any sequence k n such that lim n→∞ k n /n = κ ∈ [ 0, 1], it holds almost surely where t is the unique solution of The proof is a consequence of proposition 1 in [17] and is deferred to Appendix 1.
In Fig. 3 (experiment 2), we compare the compressibility parameters (κ, ε) of the Laplace distribution and of 2-LMM distribution with α = 0.1, β = 10 for several values of p. It should be noted that Laplace distribution is not a compressible distribution (we can not have κ and ε small at the same time), whereas 2-LMM distribution are compressible if parameter p is sufficiently small, as we can have both κ and ε small at the same time.
We now compute the empirical probability density function of the Haar wavelet coefficients of several images and the best fitting of the mixture of two Laplace density functions computed by maximizing the likelihood of the data. The computation has been carried out via expectation maximization algorithm [16]. In Fig. 4, we show the results for several images. In order to compare the two parametric representations of sparsity, in Table 1, the Kullback-Leibler divergence of the best fitting probability models and the empirical probability density function are computed for the two models. It can be noticed that a single Laplace is a poor model for the Haar For convenience, we consider p fixed as a guess of the degree of signal's sparsity, whereas = (α, β) will be unknown. The choice to keep p fixed does not entail a significant restriction to our analysis.

Proposition 2
The following optimization problems are equivalent Table 1 The Kullback-Leibler divergence of the best fitting probability models and the empirical probability density function are computed for the two models and Given y, A, instead of solving optimization problem in (7), we consider the minimization of the following modified cost function: where K = pn , and Compared to (7), the optimization problem in (10) 1. Introduces the parameter, which is a regularization term used to avoid singularities whenever one of the Laplace components collapses onto a specific data point since we expect that α ≈ 0, since we seek a sparse solution; this fact will be clear later; 2. Introduces the constraint π ∈ n−K , which enforces a sparse solution.
Similar computations can be carried out for the case with noise, leading to min min x∈R n min It should be noted that there is not a closed form solution to problems (10) and (12).
However, partial minimizations of V = V (x, π, α, β, ) with respect to just one of the variables have simple representation. More precisely, we have the following expressions (see Lemma 3 in Appendix 3).

Proposition 3 Let
and In the following section, we present several iterative algorithms to approximately solve these optimization problems.

Iterative shrinkage/thresholding algorithms
The literature describes a large number of approaches to address minimization of (2) and (3). Popular iterative methods belong to the class of iterative shrinkagethresholding algorithms. These methods can be understood as a special proximal forward backward iterative scheme [42] and are appealing as they have lower computational complexity per iteration and lower storage requirements than interior-point methods. In fact, these types of recursions are a modification of the gradient method to solve a least square problem, where the dominant computational effort lies in a relatively cheap matrixvector multiplication involving A and A and the only difference consists in the application of a shrinkage/soft thresholding operator, which promotes sparsity of the estimate at each iteration.
More precisely, let τ (t) t∈N be a sequence in (0, ∞) such that inf t∈N τ (t) > 0 and sup t∈N τ (t) < 2 A −2 2 , and let u (t) t∈N be a sequence in R n . Then, for every t ∈ N let where η S γ is a thresholding function to be applied elementwise defined as The simplest form, known as iterative shrinkagethresholding algorithms (ISTA, [6]), considers u (t) = 0 and τ (t) = τ < 2 A −2 2 for all t ∈ N. This algorithm is guaranteed to converge to a minimizer of the Lasso [6]. Moreover, as shown in [43], if A fulfills the so-called finite basis injectivity condition, the convergence is linear. However, the factor determining the speed within the class of linearly convergent algorithms depends on local well conditioning of the matrix A, meaning that ISTA can converge arbitrarily slowly in some sense, which is also often observed in practice.
In order to speed up ISTA, alternative algorithms have exploited preconditioning techniques or adaptivity, combining a decreasing thresholding strategy with adaptive descent parameter. However, the lack of a model-based thresholding policy makes this algorithm very sensitive to the signal statistics and the accuracy is not always guaranteed. In [13], the thresholding and descent parameters are optimally tuned in terms of phase transitions, i.e., they maximize the number of nonzeros at which the algorithm can successfully operate. However, preconditioning can be very expensive and there is no proof of convergence for adaptive methods.
Finally, other variations update the next iterate using not only the previous estimation, but two or more previously computed iterates. Among all the proposed techniques with a significantly better rate of convergence and phase transitions, we recall (a) fast iterative shrinkagethresholding algorithm (FISTA, [9]) obtained by (14) choosing and (b) approximate message passing (AMP, [14]) with threshold recursion proposed in [44] In this section, we show how to adapt these numerical methods to solve the weighted minimization problem via 2-LMM.

2-LMM-tuned iterative shrinkage-thresholding
Let us consider the problem of minimizing (11). Since information about the locations of the nonzero coefficients of the original signal is not available a priori, the task of selecting the parameters α, β, and π is performed iteratively. We propose an alternating method for the minimization of (11), inspired by the EM algorithm [16]. The pseudocode of the algorithm is reported in Algorithm 1. The strategy can be summarized as follows.
1. Let t := 0 and set an initial estimate K for the sparsity level, p = K/n, a small value α (0) ≈ 0 (e.g., α (0) = 0.1), the initial configuration π (0) = 1, and (0) = 1. Since π (0) = 1, β (0) can be arbitrary since it is not used in the first step of the algorithm. 2. Given the observed data y and current parameters π i , α, andβ, a new estimation x (t+1) of the signal is obtained by moving in a minimizing direction of weighted Lasso 3. The posterior distribution of the signal coefficients is evaluated and thresholded by keeping its n − K biggest elements and setting the others to zero. It is worth remarking that this step differs from the E-step of a classical EM algorithms as a thresholding operator σ n−K is applied in order to promote the sparsity in the probability vector π. 4. Given the probabilities, we use them to re-estimate the mixture parameters α (t) and β (t) . 5. Set t := t + 1 and iterate until the stopping criteria is satisfied, e.g., until the estimate stops changing

Relation to prior literature
As already observed, Algorithm 1 belongs to the more general class of methods for weighted 1 norm minimization [45][46][47] (see (18)). Common strategies for iterative reweighting 1 minimization (IRL1, [45]) that have been explored in literature re-compute weights at every iteration using the estimate at the previous iteration ω where χ and are appropriate positive constants. In Algorithm 1, the weights ω (t) i are chosen to jointly fit the signal prior and, consequently, depend on all components of the signal and not exclusively on the value x (t) i . Our strategy is also related to threshold-ISD [19] that incorporates support detection in the weighted 1 minimization and runs as fast as the basis Computation of 1 -weights: Gradient/Thresholding step: Posterior distribution evaluation: Regularization parameter: Parameters estimation: pursuit. Given a support estimate, the estimation is performed by solving a truncated basis pursuit problem. Also in [48], an iterative algorithm, called WSPGL1, is designed to solve a sequence of weighted LASSO using a support estimate, derived from the data, and updated at every iteration. Compared to threshold-ISD and WSPGL1, 2-LMM-tuned iterative shrinkage-thresholding does not use binary weights and is more flexible. Moreover, in threshold-ISD, like CoSaMP, the identification of the support is based on greedy rules and not chosen to optimally fit the prior distribution of the signal. A prior estimation based on EM was incorporated within the AMP framework also in [49] where a Gaussian mixture model is used as the parametric representation of the signal. The key difference in our approach is that model fitting is used to estimate the support and to adaptively select the best thresholding function with the minimum mean square error. The necessity of selecting the best thresholding function is also proposed in parametric SURE AMP [50] where a class of parametric denoising functions is used to adaptively choose the best-in-class denoiser. However, at each iteration, parametric SURE AMP needs to solve a linear system and the number of parameters affects heavily both performance and complexity.

Convergence analysis
Under suitable conditions, we are able to guarantee the convergence of the iterates produced by Algorithm I and discuss sufficient condition for optimality. (12) if it satisfies the following relation , then it is a local minimizer of (12).
The proof can be obtained with similar techniques, devised in [51], and we omit the proof for brevity. This result provides a necessary condition for optimality and shows that, being the function in (12) not convex, τ -stationarity points are only local minima. The next theorem ensures that also the sequence (ζ (t) ) converges to a limit point which is also a τ -stationary point of (12) of the algorithm and, from Theorem 1, a local minimum for (12). Moreover, in Theorem 3, we derive an upper bound on the reconstruction error.

Definition 3
Let A be an m×n matrix and let 1 ≤ s ≤ n be an integer. The matrix A is said to satisfy the s-restricted isometry property with restricted isometry constant δ s ∈ (0, 1) if, for every x with |supp(x)| ≤ s, it holds Theorem 3 (2LMM-ISTA: upper bound on the error) Suppose that A is an m×n sampling matrix with restricted isometry constant δ 2K . Let x ∞ be the output of Algorithm 1 with It should be noticed that the result in Theorem 3 implies that the mean square error MSE = e 2 /n = c + λ In this sense, we have provided conditions verifiable a posteriori for convergence in a neighborhood of the solution. This is a common feature in shrinkage-thresholding methods. In [52], it is shown that, in the absence of noise, if certain conditions are satisfied, the error provided by Lasso is O(λ), where λ is the regularization parameter. Since ISTA and FISTA converge to a minimum of the Lasso, we argue that the same estimate holds also for the error between the provided estimations and the true signal.
The proof of Theorems 2 and 3 are postponed to Appendices 3 and 4 and are obtained using arguments of variational analysis and analysis of τ -stationary point of (12), respectively. bound, obtained using Theorem 3, is 5.143 · 10 −10 . Using error bounds in [52], we are able to guarantee that the solutions provided by ISTA and FISTA are accurate with an error only proportional to λ = 10 −3 .

Numerical results, experiments, and discussion
In this section, we compare several first-order methods with their versions augmented by the support estimation, in terms of convergence times and empirical probability of reconstruction in the absence and in the presence of noise. It is worth remarking that this does not represent a challenge among all first-order methods for compressed sensing. Our aim is to show that the combination of support detection and estimation using an iterative reweighted first-order method can improve several iterative shrinkage methods. In other terms, we want to show that, given a specific algorithm for CS, the speed and the performance can be improved via its 2-LMM counterpart. Moreover, in order to show that the choice of the weights is important to obtain fast algorithms and good performance, we have employed an iterative shrinkage method for iterative reweighted 1 minimization algorithm (IRL1). In [45], IRL1 requires to solve at each step a weighted 1 minimization. This algorithm has computational complexity which is not comparable with the iterative shrinkage/thresholding algorithms since each iteration has complexity of order O(n 3 ). We employ a shrinkage-thresholding method for IRL1 in the spirit of [53] and show that the performance are not as good as in the proposed methods. What we mean as IRL1 is summarized in Algorithm 2. Computation of 1 -weights: Gradient/Thresholding step:

Rate of convergence
As a first experiment, we consider Bernoulli-Gaussian signals [54]. More precisely, the signal to be recovered has length n = 560 with k = 50 nonzero elements drawn from a N(0, 4), respectively. The sensing matrix A with m = 280 rows is sampled from the Gaussian ensemble with zero mean and variance 1/m. We fix λ = 10 −3 and τ = 0.19, and the mixture parameters are initialized α (0) = 1, π (0) = 1, K = k + 10, and p = K/n. In Fig. 5, we compare the convergence rate of ISTA, FISTA, IRL1, and AMP with the corresponding methods with 2-LMM-tuning (2-LMM-ISTA, 2-LMM-FISTA, and 2-LMM-AMP). In particular, the mean square error (MSE) of the iterates MSE(t) = x (t+1) − x 2 /n averaged over 50 instances is depicted as a function of the iteration number.
A few comments are in order. The sparsity problem that ISTA, FISTA, and AMP are intended to approximately solve is the same (basis pursuit or Lasso problem). However, the convergence results are different for these iterative algorithms. More precisely, in the absence of noise, • ISTA and FISTA, under certain conditions, are analytically proved to converge to a minimum of Lasso. This solution is shown to provide only an approximation of the sparse solution x which is controlled by Lasso parameter λ. More precisely, x − x 2 ≤ Cλ where C ∈ R and perfect reconstruction is not guaranteed.
• AMP instead is not proved to converge in a deterministic sense. In [14], only the average case performance analysis is carried out. The authors exploit the randomness of A and instead of calculating the limit of x t − x 2 , they show the convergence in the mean square sense E x t − x 2 → 0.
In Fig. 5 the accuracy of the solution provided by ISTA, FISTA, and AMP are different. The difference of AMP has been already explained. The difference between ISTA and FISTA is due to the fact that λ = 0.005 for ISTA (to speed up the algorithm) and λ = 0.001 for FISTA. As already observed, we are not interested in a challenge among all first-order methods for CS. Our aim is to show that the combination of support detection and estimation using an iterative reweighted first-order method can improve a series of iterative shrinkage methods. More precisely, given a specific algorithm for CS, the speed can be improved via its 2-LMM counterpart.
It should be noted that the proposed algorithms are much faster than classical iterative shrinkagethresholding methods: there is about a 81,37, and 35% of reduction in the number of iterations needed for the convergence of 2-LMM-ISTA, 2-LMM-FISTA, and 2-LMM-AMP, respectively.

Effect of the prior distribution
We now show a second experiment: we fix n = 512 and take the fraction of the nonzero coefficients fixed to ρ = k/n and we study the effect of the nonzero coefficients distribution on the empirical probability of reconstruction for different values of k ∈[ 1, 250]. More precisely, x i ∼ (1−ρ)δ 0 (x i )+ρg(x i ) where g is a probability distribution function and δ 0 is the Dirac delta function. In Table 1, the acronyms of the considered distributions are summarized (see also [55]). Figures 6,7,8, and 9 (left) show the empirical recovery success rate, averaged over 50 experiments, as a function of the signal sparsity for different signal priors (see Table 2). For all recovery algorithms, the convergence tolerance has been fixed to 10 −4 . In this case, the elements of matrix A with m = 350 are sampled from a normal distribution with variance 1/m. We have fixed a total number of iterations equal to 1000. The algorithm parameters have been initialized as in Table 3.
It should be noticed that for ISTA, IRL1, and 2-LMM-ISTA, we have chosen λ = 0.005, instead of λ = 0.001 as in FISTA and 2-LMM-FISTA, in order to speed up the convergence. The convergence of ISTA and IRL1 are extremely slow, and before 1000 iterations, we get only an approximation with an error of order 10 −3 for sparsity larger than 50 in most of the cases.
It should be noticed that the 2-LMM-tuning improves the performance of iterative shrinkage-thresholding methods in terms of sparsity-undersampling trade-off. For example for 5P1, it turns out that the signal recovery with 2-LMM-tuning is possible with 30%, 63% sparsity level higher than FISTA, and AMP, respectively. show the average running times (CPU times in seconds) of the algorithms computed over the successful experiments, and the error bar represents the standard deviation of uncertainty for different signal priors. These graphs demonstrate the benefit of 2-LMM-tuning for iterative shrinkage/thresholding methods. Not only 2-LMM-tuning shows better performance in the reconstruction but it also runs much faster than traditional methods. Despite the additional per iteration computational cost needed to update the mixture parameters, the gain of the 2-LMM-tuning ranges from 2 to over 6 times, depending on the signal prior. The algorithm efficiency can be attributed to the simple form of the model used as parametric representation of the signal and the improved runtime performance comes from the effective denoising so that fewer iterations are required to converge.

Reconstruction in imperfect scenarios
In this section, we compare the first-order methods with their versions augmented by the support estimation for recovery of signals in imperfect scenarios where the signal is not exactly sparse or the measurements are noisy.

Not exactly sparse signals
In this experiment, we investigate the performance of first-order methods with 2-LMM-tuning for signals that are not strictly sparse. We consider signals of the form x = x 0 + ξ where x 0 is drawn i.i.d from the ensemble of Bernoulli-Gaussian signals (see G1 in Table 2) of length n = 512 with sparsity level k = 56 and ξ ∈ R n is a vector whose components are distributed as N 0, σ 2 with σ = 0.01. Here, k can be interpreted as the compressibility level of the signal x . The sensing matrix A with m ∈ [160, 360] rows is sampled from the Gaussian ensemble with zero mean and variance 1/m. Then, the reconstruction is performed using measurements y = Ax . The first-order methods are compared with their versions augmented by 2-LMM-tuning: the parameters have been initialized as follows: λ = 10 −3 , τ = 0.1, α (0) = 1, π (0) = 1, K = k + 10, and p = K/n. Figure 10 shows the MSE achieved by 3000 iterations of the algorithms as a function of measurements used for the reconstruction. As we can see, the algorithms with 2-LMM-tuning have similar reconstruction performance and outperform significantly their traditional counterpart (for example 2-LMM-AMP needs 62% of measurements required by AMP to reach a similar accuracy of MSE = 10 −4 ). In this case the improved performance comes from the effective denoising so that fewer iterations are required to achieve a better accuracy.

Reconstruction with noisy measurements
We now fix n = 512 and we study the performance of the proposed methods in scenarios with inaccurate measurements according to (1). In this case, x is a random Bernoulli-Uniform signal (see U1 in Table 2) with sparsity degree k = 56 and the noise η is white Gaussian noise with standard deviation σ = 0.01. In this case, the parameters are set as follows: λ = 10 −3 , τ = 0.1, α (0) = 1, π (0) = 1, K = k + 10, and p = K/n. The MSE achieved by 3000 iterations is depicted as a function of the number of measurements used in the reconstruction. Also in this case, the best results are obtained by methods with 2-LMM-tuning. The efficiency of the proposed algorithms allows to reduce the number of measurements required to achieve a satisfactory level of accuracy. As can be noticed from Fig. 10  In Fig. 11, we show that the proposed methods are robust against noise. More precisely, the mean square error, averaged over 50 runs, and obtained after 3000 iterations, is depicted as a function of signal-to-noise ratio (SNR), defined as follows As to be expected, as the SNR increases the MSE goes to zero. Moreover, the MSE of the proposed algorithms are smaller than those obtained via classical iterative thresholding algorithms. As already observed, the MSE of ISTA is very high compared to the other methods. This is due to the fact that the algorithm is very slow and the number of iterations are not enough to reach a good recovery error.
In our setup, we have considered only Gaussian Noise and the robustness against multiplicative noise is out of our scope. This would require a drastic modification of the proposed algorithms and will be subject of future research. For example, when the underlying sparse or approximately sparse signal is a vector of nonnegative intensities whose measurements are corrupted by Poisson noise, standard CS techniques cannot be applied directly, as the Poisson noise is signal-dependent [56,57]. In this case, the rationale of our method can be adapted combining the use of mixtures models with exponential distribution instead Laplace distribution with penalized negative Poisson log-likelihood objective function with nonnegativity constraints. We refer to [58] for more details on the model and on the implementations of related iterative thresholding algorithms.
If the multiplicative noise is due to hardware's amplification and is not signal-dependent, we can model the measurements as follows  (N (0, σ 2 ))), where D is a diagonal matrix of noise and σ is the parameter governing the amplitude of decalibration. To address this problem, the most standard existing approach is the blind calibration for compressed sensing [59]. More precisely, the sparse regularization is exploited considering A as an inaccurate estimate of the true measurement system A = DA and y = Ax +( A−A)x +η ≈ Ax +ε+η with ε an estimate of the magnitude of this added noise A − A x . We refer to [59] for more sophisticated approaches of blind supervised calibration and adaptations of classical methods that perform both calibration and reconstruction.

Comparison with structured sparsity-based Bayesian compressive sensing
Many authors have recently developed structured sparsity-based Bayesian compressive sensing methods in order to deal with different signals arising in several applications and adaptively explore the statistical structure of nonzero pattern. We refer the interested reader to the repository http://people.ee.duke.edu/~lcarin/BCS. html for an introduction to Bayesian compressive sensing (BCS) methods and to structured sparsity-based Bayesian compressive sensing. For example, [60] proposes a spatio-temporal sparse Bayesian method to recover multichannel signals simultaneously, not only exploiting temporal correlation within each channel signal but also exploiting inter-channel correlations among different signals. This method has been shown to provide several advantages in applications in brain computer interface and electroencephalographybased driver's drowsiness estimation in terms of measurements for reconstruction and computational load. In [61], using a new bilinear time-frequency representation, a redesigned BCS approach is developed for the problem of spectrum estimation of multiple frequency-hopping signals, arising in various communication and radar applications in the context of multiple-input multiple-output (MIMO) operations in the presence of random missing observations. Another example of structured sparsitybased Bayesian compressive sensing comes from the context of reconstruction of signals and images that are sparse in the wavelet basis [62] or in DCT basis with applications to image compression. More precisely, in [62], the statistical structure of the wavelet coefficients is exploited explicitly using a tree-structured Bayesian compressive sensing approach. This tree structure assumption shows several advantages in terms of number of measurements required for reconstruction.
It is worth remarking that in our approach, we do not use any prior information on the structure of the sparsity pattern and we expect that structured sparsity-based methods outperform our approach. A detailed comparison and an ad hoc adaptation of our approach to all specific frameworks mentioned above is out of the scope of this paper. However, in this section, we propose a numerical comparison of 2-LMM-tuning FISTA with tree-structured wavelet-based Bayesian compressive sensing (WBCS). The implementation of WBCS algorithm used for comparison are implemented via a hierarchical Bayesian framework, with the tree structure incorporated naturally in the prior setting. See TS-BCS for wavelet via MCMC in the repository http:// people.ee.duke.edu/~lcarin/BCS.html for a detailed description of the code.
For the comparison, we consider the setting in [62] with a signal of length n = 512 that are sparse in the Haar wavelet basis and whose coefficients are not independent as in classical compressed sensing framework. Specifically, under the wavelet basis, if a parent node in a wavelet Fig. 11 Analysis of robustness against noise: the mean square error is depicted as a function of the SNR tree is zero or close to zero, with a very large probability, its children nodes are also zero or close to zero. We refer to [62] for details on the generation of the signal. The sparsity of the considered signal is k = 63. In this case, the parameters are set as follows: λ = 10 −3 , τ = 0.05, α (0) = 1, π (0) = 1, K = k + 10, and p = K/n. In Fig. 12 the reconstruction error achieved by 10000 iterations is depicted as a function of the number of measurements used in the reconstruction. It is worth remarking that WBCS explores the statistical structure of the wavelet coefficients to reduce the number of CS measurements and goes beyond simply assuming that the data are compressible in a wavelet basis. As to be expected, since more a priori knowledge is employed, WBCS shows better performance in terms of reconstruction accuracy. However, the gap is not large and the 2-LMM-FISTA tuning is able to learn the sparsity model and, as soon as the number of measurements is larger than 200, we obtain a good reconstruction accuracy, of order 10 −4 for 2-LMM-FISTA and of order 10 −5 for WBCS.

Deblurring images
In order to show the effectiveness of the 2-LMM-tuning, in this section, we repeat the same experiment proposed in [9] for deblurring two test images (Lena and cameraman). In [9], it has been shown that FISTA significantly outperforms ISTA and other first-order methods in terms of the number of iterations required to achieve a given accuracy. For this reason, we compare the performance of FISTA with our proposed algorithm 2-LMM-FISTA.
In the considered setting, both images have equal size 256 × 256 and all pixels of the original images are scaled into the range between 0 and 1. A Gaussian blur of size 9× 9 and standard deviation 4 are applied to both images and an additive zero mean white Gaussian noise with standard deviation 10 −4 is added. The original and observed images are given in Figs. 13 and 14, respectively. We then test FISTA and 2-LMM-FISTA for recovery, where y represents the (vectorized) observed image, and A = RW , where R is the matrix representing the blur operator and W is the inverse of a three-stage Haar wavelet transform. The regularization parameter is fixed as in [9] λ = 2·10 −5 , and the blurred image is used as initial condition. For 2-LMM-ITA, the parameters are set as follows: α (0) = 1, π (0) = 1, K = 10000, and p = K/n.
In Fig. 15, the evolution of the error in dB is depicted as a function of the number of iterations. In particular, the images produced by 2-LMM-FISTA exhibit better quality than those obtained by using the classical version of FISTA. In Fig. 13 and 14, the reconstructions obtained by the competing methods are shown for Lena and cameraman, respectively. As can be seen in Figs. 13 and 14, 2-LMM-FISTA achieves significantly better visual quality, as the amount of noise is minimized and visual artifacts are greatly reduced. This is also reflected by the reconstruction PSNR, which is significantly higher for 2-LMM-FISTA.

Conclusions
In this paper, we proposed a new method to perform both support detection and sparse signal estimation in compressed sensing. Combining MAP estimation with the parametric representation of the signal with a Laplace mixture model, we formulated the problem of  reconstruction as a reweighted 1 minimization. Our contribution includes theoretical derivation of necessary and sufficient conditions for reconstruction in the absence of noise. Then, 2-LMM-tuning has been proposed to improve the performance of several iterative shrinkagethresholding algorithms. Iterative procedures have been designed by combining EM algorithm with classical iterative thresholding methods. Numerical simulations show that these new algorithms are faster than classical ones and outperform them in terms of phase transitions. Topics of our current research is to use similar technique based on Laplace mixture models for robust compressed sensing, where measurements are corrupted by outliers (see [63] and reference therein).

Appendix 1: proof of Proposition 1
The proof of Proposition 1 is a direct consequence of a more general result on compressible prior result, formally stated in Proposition 1 in [17]. Lemma 1 (Proposition 1.1 in [17]) Suppose x n ∈ R n is i.i.d. with respect to a distribution p(x). Denote p(x) := 0 for x < 0, and p(x) := p(x) + p(−x) for x ≥ 0 as the probability density function of |X n |, and F(t) := P(|X| ≤ t) as its cumulative density function. Assume that F is continuous and strictly increasing on some interval [ a, b], with F(a) = 0 and F(b) = 1, where 0 ≤ a ≤ b ≤ ∞. For any 0 < κ ≤ 1, define the following function is also well defined for κ = 0, and for any sequence k n such that lim n→∞ k n /n = κ ∈[ 0, 1] the following holds almost surely Proof of Proposition 1. For 2-LMM, we have Then, let t = t(κ) = F −1 (1 − κ), i.e., be the solution We now compute Then, the assertion is proved by applying Lemma 1 and we obtain
Proof of Proposition 2 . Let us consider J(x, π; ) and minimize with respect to π i by taking x all the other variable fixed. By imposing the constraint and the minimizing value is given by From last equality and from Proposition 4, we conclude the thesis.

Appendix 3: proof of Theorem 2
In this section, we prove rigorously Theorem 2, which guarantees the convergence of 2-LMM-ISTA to a limit point. We start from the following preliminary results. Let V : R n × n−K × R × R × R → R be the function defined in (12) V (x, π, α, β, ) where H : [0, 1] → R is the natural entropy function . and J is defined in (11).
Let a be the vector satisfying given by is the nonincreasing rearrangement of a. It should be noticed that arg min being γ π 0 = (n − K)γ just a constant for π ∈ n−K . The minimum of (23) can be calculated by minimizing with respect to each π i individually and V (x, π, α, β, ) + γ π 0 = g(π i ) + C, where C is independent of π i and To derive the minimum, we distinguish two cases, π i = 0 and π i = 0. In the first case, the element-wise cost is (ignoring the constant terms) 0. In the second case, the minimum cost (again ignoring the constant terms) is attained for π i = a i if π i = 0. Comparing the cost for both cases, i.e, g(a i ) < 0, we obtain By definition of γ , we get arg min From this result and the fact that we conclude that π is the minimizing value of V for fixed x, α, β, .

Lemma 4 Define the surrogate functional
then a, π, α, β, ) with Proof By developing the least squares in (12) is straightforward to show that where χ is a function independent of x. By differentiating the function with respect x, we obtain the thesis. (22) is not increasing along the iterates

Proposition 5 The function V defined in
Proof From Lemma 3 and 4, it should be noticed that, for each time t ∈ N, we have Since V (x, π, α, β, ) is an increasing function in and being (t+1) = min{ (t) , θ (t+1) } ≤ (t) by definition, we obtain (t) and therefore, using π (t+1) = arg min see Lemma 3, we get It should be noticed that for all x V S (x, x, π, α, β, ) = V (x, π, α, β, ) and V S (x, a, π, α, β, ) ≥ V (x, x, π, α, β, ) for all a = x. Then, we have The following lemma implies that these algorithms converge numerically when the number of iterations goes to infinity.
We conclude that for all ∈ N T t=1 1 2τ By letting → ∞, we obtain that the series is convergent and, we obtain the necessary condition 0 ≤ 1 2τ as t → ∞ and from inequality in (27) the assertion is proved.
(f) Follows by noticing that the series is telescopic (g) Is a direct consequence bound in (26)

Lemma 6 The sequence (x (t) ) t∈N is bounded
Proof We now prove that both α (t) and β (t) must be upper bounded. From (25), there exists a constant C = V x (1) , π (1) , α (1) , β (1) , (1) λ + n log 2 such that zero when → ∞, there exists t 0 ∈ N and a constant χ ∈ R such that if t > t 0 then Iterating the argument and letting → ∞, and we conclude that the sequence (x (t ) ) ∈N is bounded and so is (β (t ) ) from which we get the contradiction. We conclude that (α (t) ) t∈N and (β (t) ) t∈N are both upper bounded by a constant χ > 0 and so (x (t) ) t∈N : (1 − π i )β (t) ≤ χn Proposition 6 Any accumulation point is a τ -stationary point of (12) of the algorithm and satisfies the equalities in (19a)-(19e).
(a) There exists a sequence q q∈N such that 0 = π t q i = a t q i → π i = 0. This means that there exists q 0 ∈ N such that ∀q > q 0 we have a (t q ) i < a t q j , ∀j ∈ supp π and, by letting q → ∞, (b) There exists 0 ∈ N such that ∀ > 0 π (t ) i = 0, from which a j , ∀j ∈ supp π and by letting → ∞ we get (30).