About this Journal Submit a Manuscript Table of Contents
Mathematical Problems in Engineering
Volume 2012 (2012), Article ID 478931, 22 pages
http://dx.doi.org/10.1155/2012/478931
Research Article

Sparse Signal Recovery via ECME Thresholding Pursuits

1School of Computer Science and Telecommunication Engineering, Jiangsu University, Zhenjiang 212013, China
2School of Information Science and Technology, Sun Yat-Sen University, Guangzhou 510006, China

Received 17 February 2012; Revised 24 April 2012; Accepted 8 May 2012

Academic Editor: Jung-Fa Tsai

Copyright © 2012 Heping Song and Guoli Wang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The emerging theory of compressive sensing (CS) provides a new sparse signal processing paradigm for reconstructing sparse signals from the undersampled linear measurements. Recently, numerous algorithms have been developed to solve convex optimization problems for CS sparse signal recovery. However, in some certain circumstances, greedy algorithms exhibit superior performance than convex methods. This paper is a followup to the recent paper of Wang and Yin (2010), who refine BP reconstructions via iterative support detection (ISD). The heuristic idea of ISD was applied to greedy algorithms. We developed two approaches for accelerating the ECME iteration. The described algorithms, named ECME thresholding pursuits (EMTP), introduced two greedy strategies that each iteration detects a support set I by thresholding the result of the ECME iteration and estimates the reconstructed signal by solving a truncated least-squares problem on the support set I. Two effective support detection strategies are devised for the sparse signals with components having a fast decaying distribution of nonzero components. The experimental studies are presented to demonstrate that EMTP offers an appealing alternative to state-of-the-art algorithms for sparse signal recovery.

1. Introduction

Sparsity exploiting has recently received a great amount of attention in the applied mathematics and signal processing community. Sparse signal processing algorithms have been developed for various applications. A recent Proceedings of the IEEE special issue on applications of sparse representation and compressive sensing devoted to this topic. Some of the exciting developments addressed in [17]. Compressed sensing, also known as compressive sensing, compressive sampling (CS) [8, 9], has been the subject of significant activity in sparse signal processing where one seeks to recover efficiently a sparse unknown signal vector of dimension 𝑛 via a much smaller number (𝑚) of undersampled linear measurements. For 𝑘-sparse unknown signal 𝑥0𝑛, the sparse signal recovery is intimately related to solving an underdetermined system of linear equations 𝑦=𝐴𝑥0 with sparseness constraint 𝑃0:min𝑥𝑥0s.t.𝐴𝑥=𝑦,(1.1) where 𝑥 is true signal to be recovered, 𝑥0 is used to denote the number of nonzero components of 𝑥, 𝐴𝑚×𝑛 is the measurement matrix, and 𝑦𝑚 is the measurement vector. The key insight in the pioneering work on CS [8, 9] is to replace 0 by 1 for the (𝑃0) problem due to nonconvexity and combinatorial effect. In [10], it is the basis pursuit problem (BP):min𝑥0𝑥01s.t.𝐴𝑥=𝑦.(1.2) Hereafter, numerous researchers developed various computational algorithms for sparse signal recovery [11, 12]. Generally, there are three major classes of computational algorithms: convex relaxation, bayesian inference, and greedy pursuit. Convex relaxation replaces the combinatorial problem (𝑃0) with a convex optimization problem (BP), such as basis pursuit [10], basis pursuit denoising [13], the least absolute shrinkage and selection operator (LASSO) [14], and least angle regression (LARS) [15]. Bayesian inference derives the maximum a posteriori estimator from a sparse prior model, such as sparse bayesian learning (SBL) [16, 17] and bayesian compressive sensing (BCS) [18, 19]. Another popular approach is to use greedy algorithms which iteratively refine a sparse solution by successively selecting one or more elements. These algorithms include matching pursuit (MP) [20], orthogonal matching pursuit (OMP) [2123], subspace pursuit (SP) [24], compressive sampling matching pursuit (CoSaMP) [25], and iterative thresholding algorithms [26, 27].

Iterative hard thresholding (IHT) [27] is a simple and powerful approach for sparse recovery. Recently, an alternative algorithm has been present to alleviate the convergence speed issue of IHT in [28, 29]. Qiu and Dogandzic [28, 29] derived an expectation conditional maximization either (ECME) iteration from a probabilistic framework based on the Gaussian-distributed signals with unknown variance and proposed an acceleration method, termed double overrelaxation (DORE) thresholding scheme, to improve the convergence speed of the ECME algorithm. In addition, Qiu and Dogandzic [28, 29] further proposed an automatic double overrelaxation (ADORE) thresholding method for conditions that the underlying sparsity level is not available. As in study [30], Wang and Yin presented an iterative support detection (ISD) algorithm to refine the failed reconstructions by thresholding the solution of a truncated (BP) problem.

Inspired by the theoretical and empirical evidence of favorable recovery performance of ECME [29] and ISD [30], we combine ECME [29] and ISD [30] to devise novel sparse signal recovery methods, dubbed as ECME thresholding pursuits (EMTP). EMTP detects a support set 𝐼 using an ECME iteration and estimates the reconstructed signal by solving a truncated least-squares problem on the support set 𝐼, and it iterates these two steps for a small number of times. We present two effective support detection strategies (hard thresholding and dynamic thresholding) for the sparse signals with components having a fast-decaying distribution of nonzeros (called fast decaying signals in [30]), which include sparse Gaussian signals, sparse Laplacian signals, and certain power-law decaying signals.

This paper considers the iterative greedy algorithms and abstracts them into two types [31], One-Stage Thresholding (OST) and Two-Stage Thresholding (TST), as discussed in Sections 2.1 and 2.2. Then, we review the initial work of Iterative Support Detection (ISD) in Section 2.3. In Section 3, we describe the proposed approaches. After that, Section 4 details the experimental results. Finally, we conclude this paper in Section 5.

1.1. Notation

We introduce the notation used in this paper.(i)𝑥(𝑡): the algorithms described in this paper are iterative and the reconstructed signal 𝑥 in current iteration 𝑡 is denoted as 𝑥(𝑡). The same convention is used for other vectors and matrices.(ii)𝐼, 𝐴𝐼: index set 𝐼, the matrix 𝐴𝐼 denotes the submatrix of 𝐴 containing only those columns of 𝐴 with indexes in 𝐼. The same convention is used for vectors.(iii)[1,𝑛]𝐼: the complement of set 𝐼 in set {1,2,,𝑛}.(iv)supp(𝑥): the support set of a vector 𝑥, that is, the index set corresponding to the nonzeros of 𝑥, supp(𝑥)={𝑖𝑥𝑖0}.(v)𝐻𝑘(𝑥): the hard thresholding that sets all but the largest in magnitude 𝑘 elements of a vector 𝑥 to zero.(vi)|𝑥|,𝑥𝑝,𝑥𝑇: the absolute value, 𝑝 norm, and transpose of a vector 𝑥, respectively.(vii)𝐴: the Moore-Penrose pseudoinverse of matrix 𝐴𝑚×𝑛. 𝐴=𝐴𝑇(𝐴𝐴𝑇)1 for 𝑚𝑛; 𝐴=(𝐴𝑇𝐴)1𝐴𝑇 for 𝑚𝑛.

2. Related Works

2.1. One-Stage Thresholding

Qiu and Dogandzic derived an expectation conditional maximization either (ECME) iteration based on a probabilistic framework [29] 𝑥(𝑡)=𝐻𝑘𝑥(𝑡1)+𝐴𝑇𝐴𝐴𝑇1𝑦𝐴𝑥(𝑡1).(2.1) Note that ECME iteration reduces to IHT step when the measurement matrix 𝐴 has orthonormal rows (i.e., (𝐴𝐴𝑇) is the identity matrix). These one-stage thresholding algorithms (e.g., IHT [27] and ECME [29]) are guaranteed to recover sparse signals and converge with limited iterations. However, OST is not the empirical choice for practical applications due to slow convergence. To this end, Qiu and Dogandzic proposed an acceleration method, termed double overrelaxation (DORE) thresholding scheme [28, 29], to improve the convergence speed of the ECME algorithm. DORE utilizes two overrelaxation steps: 𝑥1(𝑡)=𝑥(𝑡)+𝛼1𝑥(𝑡)𝑥(𝑡1),𝑥2(𝑡)=𝑥1(𝑡)+𝛼2𝑥1(𝑡)𝑥(𝑡2),(2.2) where 𝛼1,𝛼2 are the line search parameters. Finally, an additional hard thresholding step 𝑥(𝑡)=𝐻𝑘𝑥2(𝑡)(2.3) ensures that the resulting signal is guaranteed to be 𝑘-sparse. In addition, Qiu and Dogandzic further presented an automatic double overrelaxation (ADORE) thresholding method for conditions that the underlying sparsity level is not available.

2.2. Two-Stage Thresholding

The algorithms described in this paper fall into the category of a general family of iterative greedy pursuit algorithms. Following [31], we adopt the name “Two-Stage Thresholding” (TST). Considering a CS recovery problem, the sparse recovery algorithms aim to detect the support set 𝐼 and approximate 𝑦 using 𝑦=𝐴𝐼𝑥𝐼.(2.4) Starting with initial solution 𝑥=0, TST iterates between 2 main steps:

Step 1 (support detection). Detect the support set 𝐼 of the signal 𝑥, that is, select atoms of measurement matrix 𝐴 which have been used to generate 𝑦; in other words, determine active atoms in sparse representation of a signal 𝑥. In some literature, this step also is called basis selection or atom selection.

Step 2 (signal estimation). Update the signal 𝑥 using the least-squares solution on the detected support set 𝐼. 𝑥𝐼=argminz{𝑦𝐴𝐼𝑧22,supp(𝑧)𝐼}, 𝑥[1,𝑛]𝐼=0.

Many algorithms (e.g., Orthogonal Matching Pursuit (OMP) [23], Subapace Pursuit (SP) [24], Compressed Sensing Matching Pursuit (CoSaMP) [25], and Gradient Pursuits (GP) [32]) developed various approaches for Step 1 or Step 2.

2.3. Iterative Support Detection

Considering the failed reconstructions of BP, Wang and Yin [30] proposed an algorithmic framework to refine the BP constructions, called iterative support detection (ISD). ISD alternates between two steps: support detection and signal reconstruction. Initialize the detected support 𝐼= and set the iteration number 𝑡=0; ISD iteratively calls the following steps:(1)signal reconstruction:solve the truncated BP with 𝑇=𝐼𝐶: 𝑥(𝑡)=argmin𝑥𝑥𝑇1s.t.𝑦=𝐴𝑥;(2.5)(2)support detection:detect support set 𝐼 using 𝑥(𝑡) as the reference.

The reliability of ISD relies on the support detection. Wang and Yin devised serval detection strategies for the sparse signals with components having a fast decaying distribution of nonzero components (called fast decaying signals [30]). One of the detection strategies is based on thresholding (we use ISD defined by (2.6) to refer the implementation algorithm in the following context): 𝐼(𝑡)=||𝑥𝑖𝑖(𝑡)||>𝛽𝑡||𝑥max(𝑡)||,𝛽(0,1).(2.6)

3. ECME Thresholding Pursuits

In this section, we describe our proposed approaches, named ECME Thresholding Pursuits (EMTP), that combine OST and TST using the heuristic idea of ISD. The detailed description of the proposed algorithms are presented as follows

Step 1 (initialization). We have the following. Initialize the reconstruction signal 𝑥(0)=0, initialize residual signal 𝑟=𝑦, initialize support set 𝐼(0)=, set the iteration counter 𝑡=1.

Step 2 (support detection). We have the following. Update signal approximation: 𝑥(𝑡)=𝑥(𝑡1)+𝐴𝑟(𝑡1).(3.1) Detect the support set 𝐼: strategy 1: hard thresholding 𝐼(𝑡)=supp(𝐻𝑘(𝑥(𝑡))); strategy 2: dynamic thresholding 𝐼(𝑡)={𝑖|𝑥𝑖(𝑡)|>𝛽𝑡max|𝑥(𝑡)|}.

Step 3 (signal estimation). We have the following.Estimate the signal: 𝑥𝐼(𝑡)(𝑡)=𝐴𝐼(𝑡)𝑥𝑦,(𝑡)[1,𝑛]𝐼(𝑡)=0.(3.2)Update the residual: 𝑟(𝑡)=𝑦𝐴𝑥(𝑡).(3.3)

Step 4 (halting). Check the stopping condition and, if it is not yet time to stop, increment 𝑡=𝑡+1 and return to Step 2. If it is time to stop, the recovered signal 𝑥 has nonzero entries in support set 𝐼(𝑡) and corresponding support vector lies in 𝑥𝐼(𝑡)(𝑡).

We present two thresholding proposals and corresponding algorithm EMTP-𝑘 (shown as Algorithm 1) and EMTP-𝛽 (shown as Algorithm 2), where 𝛽(0,1) is the thresholding parameter.

alg1
Algorithm 1: EMTP-𝑘 algorithm.

alg2
Algorithm 2: EMTP-𝛽 algorithm.

Remarks. (1) EMTP updates the leading elements of 𝑥(𝑡) using a least-squares solution on the support set 𝐼. However, DORE updates all the elements of 𝑥(𝑡) using double overrelaxation steps. Finally, DORE needs a hard thresholding to ensure the new approximation is 𝑘-sparse. EMTP-𝛽 demands no prior knowledge about the underlying sparsity level 𝑘.
(2) EMTP is a special case of TST and utilizes OST in the support detection stage. It is different from OMP, SP, and CoSaMP. When the ECME iteration reduces to IHT step, that is, the reference for support detection 𝑥(𝑡)=𝑥(𝑡1)+𝐴𝑇𝑟(𝑡1) is the gradient descent step for least-squares, OMP, SP and CoSaMP use the correlation between the residual signal and the atoms of the measurement matrix A (i.e., 𝐴𝑇𝑟(𝑡1)) to detect support set. 𝐴𝑇𝑟(𝑡1) is the negative gradient for least squares. It is clear that ECME iteration provides the estimation for the underlying sparse signal. We can employ the heuristic idea of ISD to devise various support detection strategies depending on the underlying sparse signal distribution. The dynamic thresholding method can be performed without the sparsity level 𝑘. EMTP directly detects the support of the underlying sparse signal by referencing the ECME iteration while OMP, SP, and CoSaMP augment the support by picking out the leading values of the negative gradient. OMP each step spots one index into the support, so it requires more iterations than EMTP. SP and CoSaMP spot several indexes into the support, so they need an additional step (i.e., orthogonal projection) to make sure the recovered signal is 𝑘-sparse. Like ISD, EMTP has an appealing self-correction capacity. An EMTP demo is presented in Figure 1.
(3) Like other greedy pursuits such as Subapace Pursuit [24] and CoSaMP [25], EMTP-𝑘 fixes the cardinality of support set 𝐼 and removes previous false detections. However, EMTP-𝛽 refines the support set 𝐼 which is not necessarily increasing and nested over the iterations.
(4) The dynamic thresholding strategy used in EMTP-𝛽 is inspired by ISD [30]. It finds significant nonzeros by comparing a threshold rather than maintaining fixed number (𝑘) of items. It is appealing for the conditions that the underlying sparsity level 𝑘 is not available.
(5) EMTP resembles ISD [30]. EMTP and ISD have the same idea in support detection step with iterative behavior. However, the support detection step is based on different sparse recovery methods, ECME and BP, respectively. EMTP updates the reconstruction signal using a least-squares solution on the detected support set 𝐼. However, ISD iteratively refines the BP solution on the complement of the detected support set 𝐼.
(6) EMTP-𝛽 with large 𝛽 obtains high-quality reconstruction from a small number of measurements. However, because support set 𝐼 grows slowly, EMTP-𝛽 takes a larger number of iterations (to be discussed in Section 4). EMTP-𝑘 wrongly detected elements can often be pruned out in later iterations.
(7) Like ISD (as discussed in [30]), EMTP-𝛽 only performs well for fast decaying signals. It does not work on sparse signals that decay slowly or have no decay at all (e.g., trinary and binary sparse signals). EMTP-𝑘 performs worse for the sparse signals that nonzero elements have similar magnitude [31, 33]. However, we can apply EMTP to non-fast-decaying signals via linear or nonlinear premapping [34, 35].

fig1
Figure 1: An EMTP-𝛽 demo that recovers a Gaussian-distributed sparse vector with 𝑛=200,𝑚=80,𝑘=20,𝛽=0.5.
3.1. Complexity Analysis

Like DORE algorithm, the basic operation is matrix-vector multiplication and sorting operation. Assume that the pseudoinverse matrix 𝐴 is precomputed (with the cost of 𝒪(𝑚2𝑛+𝑚𝑛)). Line 7 in Algorithms 1 and 2 for updating the ECME iteration requires 𝒪(𝑚𝑛). The computation complexity of thresholding operation (line 8 in Algorithms 1 and 2) is 𝒪(𝑛log𝑘). Line 9 involves the partial least-square solver, costing 𝒪(𝑘2𝑚+𝑘𝑚). The step for updating residual requires an additional matrix-vector operation, costing 𝒪(𝑘𝑚). To summarize, one EMTP iteration costs 𝒪(𝑚𝑛+𝑛log𝑘+𝑘2𝑚+2𝑘𝑚), which is significantly less than the complexity of one DORE step (with the cost of 𝒪(2𝑛log𝑘+2𝑚2𝑛+3𝑚𝑛)).

4. Experimental Results

To assess the performance of the proposed approaches, we conduct numerical experiments on computer simulations. All algorithms were implemented and tested in MATLAB v7.6 running on Windows XP with 2.53 GHz Intel Celeron CPU and 2 GB of memory. We compared the proposed approaches with the accelerated ADORE/DORE approaches. The code of ADORE/DORE is available on the authors homepage (http://home.eng.iastate.edu/~ald/DORE.htm). We initialized 𝑥 with zero sparse signal. The search length of ADORE was set to 1. All results were averaged over 100 times Monte Carlo problem instances. We used our unoptimized code (http://cs-notes.googlecode.com/files/EMTP.zip). The least-squares solution 𝑥=argmin𝑦𝐴𝑥22  was implemented using MATLAB pseudoinverse function by x = pinv(A)*y. The MATLAB code was partially adapted from [30, 33, 36]. In all the experiments, the measurement matrix 𝐴 was generated by uniform spherical ensemble (USE), that is, we generate the measurement matrix by sampling each entry independently from the standard norm distribution and then normalize each column to have unit norm. The MATLAB code is given by

A = randn(m, n); A = A./repmat(sqrt(sum(A.^2)), [m 1]).

The underlying 𝑘-sparse vectors were generated by randomly selecting a support set of size 𝑘 and each entry in the support set is sampled uniformly from a specific distribution. The sparse signals were generated in MATLAB by

x = zeros(n,1); support = randperm(n); x(support(1 : k))= v;

and the nonzeros v generated by following code.

The sparse Gaussian signals were generated in MATLAB by

v = randn(k,1);

The sparse Laplacian signals were generated by

z = rand([k, 1]);

v = zeros([k, 1]);

in = z <= 0.5;

ip = z > 0.5;

v(in)=1/lambda*log(2 * z(in));

v(ip)= -1/lambda*log(2 * (1-z(ip))).

The power-law decaying signals were generated by

v = sign(randn(k,1)).* ((1 : k). ^(-1/lambda))’.

The variable lambda controls the rate of decay. We set lambda=10 and lambda=0.5 for sparse Laplacian signals and power-law decaying signals, respectively.

For fair comparison, we stopped iterations once the relative error fell below a certain convergence tolerance or the number of iterations is greater than 100. The convergence tolerance is given by 𝑦𝑟(𝑡)2𝑦2106.(4.1) We empirically evaluate reconstruction performance in terms of signal-to-noise ratio (SNR) and number of iterations. The SNR is defined as SNR(dB)=20log10𝑥0𝑥2𝑥02,(4.2) where 𝑥0 is the true signal and 𝑥 is the recovered signal.

4.1. Influence of Thresholding Parameter

We evaluated the influence of thresholding parameter 𝛽 by varying 𝛽 from 0.4 to 0.7. First, we fixed signal dimension 𝑛=400, the underlying sparsity level 𝑘=20, and varied the number of measurements. The plot of influence of thresholding parameter 𝛽 for sparse Gaussian signals was presented in Figure 2. It is clear that larger thresholding parameter achieves better recovery performance with fewer measurements. However, the cost is more iterations. Large 𝛽 is particularly competitive when the number of measurements is fairly small. It is worth noting that the number of iterations is smaller than the underlying sparsity level 𝑘 for exact recovery, especially for large 𝑘. As shown in Figure 2, EMTP-𝛽 achieves the SNR (dB) around 300, that is, the relative error is almost as low as the double precision. Second, we tested the comparisons by fixing the number of measurements 𝑚=80, as shown in Figure 3, respectively. Finally, we enlarged the dimension of signals (𝑛=2000), as shown in Figures 4 and 5. The latter tests come to similar conclusions as the first test.

fig2
Figure 2: Influence of thresholding parameter 𝛽 for sparse Gaussian signals with 𝑛=400,𝑘=20: comparisons in terms of SNR (dB) and number of iterations.
fig3
Figure 3: Influence of thresholding parameter 𝛽 for sparse Gaussian signals with 𝑛=400, 𝑚=80: comparisons in terms of SNR (dB) and number of iterations.
fig4
Figure 4: Influence of thresholding parameter 𝛽 for sparse Gaussian signals with 𝑛=2000, 𝑘=100: comparisons in terms of SNR (dB) and number of iterations.
fig5
Figure 5: Influence of thresholding parameter 𝛽 for sparse Gaussian signals with 𝑛=2000, 𝑚=400: comparisons in terms of SNR (dB) and number of iterations.

We fixed the ratios of 𝑚/𝑛, 𝑘/𝑚 and present plots cases for the number of iterations as a function of problem size 𝑛. Given sufficient measurements, as shown in Figure 6 with 𝑚/𝑛=0.4, 𝑘/𝑚=0.2, n = 200 : 400 : 2200, EMTP-𝑘 requires significantly fewest iterations and EMTP-𝛽 requires relatively more iterations than DORE. It is clear that the number of iterations for EMTP-𝑘 and DORE is stable to the problem size, and EMTP-𝛽 needs acceptably a few more iterations with increasing the problem size. For another case, as shown in Figure 7 with 𝑚/𝑛=0.3, 𝑘/𝑚=0.3, n = 200 : 400 : 2200, DORE requires more iterations, EMTP-𝑘 needs significantly fewer iterations with increasing the problem size, and EMTP-𝛽 requires relatively stable iterations.

478931.fig.006
Figure 6: Number of iterations as a function of problem size with fixed ratios of 𝑚/𝑛=0.4, 𝑘/𝑚=0.2, n = 200 : 400 : 2200.
478931.fig.007
Figure 7: Number of iterations as a function of problem size with fixed ratios of 𝑚/𝑛=0.3, 𝑘/𝑚=0.3, n = 200 : 400 : 2200.
4.2. Comparisons with the Accelerated OST

We present comparisons in terms of SNR (dB) and number of iterations. We omit the comparisons with 1 minimization methods and some other OST/TST methods, as they have already been investigated in [31, 33]. First, we fixed signal dimension 𝑛=400, the underlying sparsity level 𝑘=20, and varied the number of measurements. The result was plotted in Figure 8. Next, we tested the comparisons by fixing the number of measurements 𝑚=80, as shown in Figure 9. Finally, the test set used larger dimension signals (𝑛=2000). The corresponding results were depicted in Figures 10 and 11. Figures 10, and 11 show the average SNR (dB) (a) and number of iterations (b) for each indeterminacy level 𝑚/𝑛. EMTP achieves significantly larger SNR (dB) than ADORE/DORE. For exact recovery, EMTP obtains the relative error almost as low as the double precision. EMTP-𝛽 constantly needs small number of iterations and EMTP-𝑘 needs the smallest number of iterations for exact recovery. However, ADORE needs hundreds of iterations. Figure 8 was zoomed in for better illustration. Figures 8, 9, 10 and 11 show the average SNR (dB) and number of iterations for each sparsity level 𝑘/𝑚. For larger-dimension signals (𝑛=2000), the results are depicted in Figures 15 and 16. The latter tests come to similar conclusions as the first test. These figures do not fully capture the performance of EMTP-𝛽 since we only set 𝛽=0.5; however EMTP-𝛽 achieved almost superior performance than ADORE/DORE.

fig8
Figure 8: Sparse Gaussian signals with 𝑛=400, 𝑘=20: comparisons in terms of SNR (dB) and number of iterations.
fig9
Figure 9: Sparse Gaussian signals with 𝑛=400, 𝑚=80: comparisons in terms of SNR (dB) and number of iterations.
fig10
Figure 10: Sparse Gaussian signals with 𝑛=2000, 𝑘=100: comparisons in terms of SNR (dB) and number of iterations.
fig11
Figure 11: Sparse Gaussian signals with 𝑛=2000, 𝑚=400: comparisons in terms of SNR (dB) and number of iterations.
4.3. Performance for Other Fast-Decaying Signals

As discussed in [30, 33], thresholding-based sparse recovery methods achieve superior performance only if the nonzero elements of signals have a fast-decaying distribution. We present empirical studies on other fast-decaying signals. Sparse Laplacian signals containing 𝑘=20 nonzero elements for dimension 𝑛=400 were tested and the corresponding result was plotted in Figure 12. For larger dimension signals (𝑛=2000), the results are shown in Figure 13. Power-law decaying signals containing 𝑘=20 nonzero elements for dimension 𝑛=400 were tested and the corresponding result was plotted in Figure 14. For larger-dimension signals (𝑛=2000), the results are depicted in Figure 15. ADORE/DORE were derived from a probabilistic framework based on Gaussian distribution. Surprisingly, they can also work for fast-decaying signals, as depicted in Figures 12, 13, 14, and 15. The conclusions for sparse Gaussian signals also can be generalized for sparse Laplacian signals and power-law decaying signals. Power-law decaying signals for low indeterminacy level 𝑚/𝑛 make an exception that EMTP-𝛽 achieves inferior performance.

fig12
Figure 12: Sparse Laplacian signals with 𝑛=400, 𝑘=20: comparisons in terms of SNR (dB) and number of iterations.
fig13
Figure 13: Sparse Laplacian signals with 𝑛=2000, 𝑘=100: comparisons in terms of SNR (dB) and number of iterations.
fig14
Figure 14: Power-law decaying signals with 𝑛=400, 𝑘=20: comparisons in terms of SNR (dB) and number of iterations.
fig15
Figure 15: Power-law decaying signals with 𝑛=2000, 𝑘=100: comparisons in terms of SNR (dB) and number of iterations.
478931.fig.0016
Figure 16: Probability of exact recovery for DORE, SP, and EMTP as a function of problem sparsity and four problem indeterminacies from the thickest to thinnest lines: 𝛿=𝑚/𝑛=0.05,0.2,0.35,0.5.
4.4. Phase Transitions

Following [31, 33], we present numerical comparisons in terms of phase transitions. For the sparse signal recovery from compressed measurements, the indeterminacy 𝛿=𝑚/𝑛 defines the undersampling of a vector in compressed measurements. Let 𝜌=𝑘/𝑚 be a normalized measure of the sparsity. We fix the problem size 𝑛=400 and test 100 Monte Carlo realizations of each sparse vector at each pair of (𝜌,𝛿). We test a grid of 16×16 linearly spaced (𝜌,𝛿) combinations with 𝜌,𝛿 varying from 0.05 to 0.5 in 16 steps. For each sparsity and indeterminacy pair, then, we find the probability of exact recovery with defining exact recovery when 𝑥0𝑥2/𝑥2<0.01. For each problem indeterminacy, we interpolate the results over all sparsities to find where successful recovery occurs with a probability of 0.5. As a result, we obtain a phase transition plot showing the boundary above which most recoveries fail and below which most recoveries succeed. Figure 16 presents the recovery rates of DORE, SP, and EMTP-𝛽 for Gaussian-distributed sparse vectors as a function of problem sparsity and four problem indeterminacies from the thickest to thinnest lines. In Figure 16, we observe EMTP-𝛽 outperforms SP, and DORE in recovery rates. Figure 17 shows the phase transitions of DORE, SP and EMTP-𝛽 for compressively sensed sparse vectors sampled from Gaussian distribution. Figure 17 indicates that these transitions obey the following hierarchy in recovery performance: EMTP-𝛽> SP > DORE.

478931.fig.0017
Figure 17: Comparison of phase transitions of DORE, SP, and EMTP for Gaussian-distributed sparse vectors.
4.5. Summary

To sum up, we have compared EMTP with ADORE/DORE. EMTP has significant advantages over the accelerated OST in terms of SNR and number of iterations. EMTP can significantly reduce the number of iterations required and achieves significantly higher SNR. For low indeterminacy level 𝑚/𝑛, EMTP requires fewer measurements. EMTP-𝛽 can work with no prior knowledge of the underlying sparsity level 𝑘 yet achieves recovery performance better than ADORE. Furthermore, we compared EMTP with the state-of-the-art greedy algorithm SP in terms of phase transitions. Among all methods, EMTP-𝛽 appeared to be the best.

5. Conclusions and Future Work

In this paper, we propose ECME thresholding pursuits (EMTP) for sparse signal recovery. EMTP detects a support set 𝐼 using the ECME iteration and estimates the reconstructed signal by solving a truncated least-squares problem on the support set 𝐼, and it iterates these two steps for a small number of times. We present two effective support detection strategies (hard thresholding and dynamic thresholding) for the sparse signals with components having a fast-decaying distribution of nonzero components. The experimental studies are presented to demonstrate that EMTP offers an attractive alternative to state-of-the-art algorithms for sparse signal recovery. EMTP can significantly reduce the number of iterations required. We then turn to the problem of reducing the computational cost of each iteration, especially for large scale applications. Future research includes three directions.(1)EMTP requires the precomputation and storage of the pseudoinverse matrix 𝐴. So it is the advisable choice that replacing ECME iteration by IHT step in the support detection stage. IHT does not require the computation and storage of the pseudoinverse matrix 𝐴, which leads to significant advantage for large scale applications.(2)In signal estimation stage, that is, update reconstructed signal by solving least-squares problem, EMTP uses the orthogonal projection, that is, by calculating 𝑥𝐼=𝐴𝐼𝑦 where 𝐴𝐼 is the pseudoinverse of 𝐴𝐼. We will approximate the orthogonal projection efficiently by gradient pursuits [32].(3)This paper devotes efforts to devise computational algorithms which are experimentally reproducible and provides empirical studies as useful guidelines for practical applications. Further, future investigations will test other OST methods as the reference and sophisticated thresholding methods for support detection.

Acknowledgment

This work was supported by the National Science Foundation of China under Grant no. 61074167, 61170126, and the Scientific Research Foundation for Advanced Talents by the Jiangsu University, no. 12JDG050.

References

  1. M. Elad, M. A. T. Figueiredo, and Y. Ma, “On the role of sparse and redundant representations in image processing,” Proceedings of the IEEE, vol. 98, no. 6, pp. 972–982, 2010. View at Publisher · View at Google Scholar · View at Scopus
  2. M. J. Fadili, J. L. Starck, J. Bobin, and Y. Moudden, “Image decomposition and separation using sparse representations: an overview,” Proceedings of the IEEE, vol. 98, no. 6, pp. 983–994, 2010. View at Publisher · View at Google Scholar · View at Scopus
  3. M. D. Plumbley, T. Blumensath, L. Daudet, R. Gribonval, and M. E. Davies, “Sparse representations in audio and music: from coding to source separation,” Proceedings of the IEEE, vol. 98, no. 6, pp. 995–1005, 2010. View at Publisher · View at Google Scholar · View at Scopus
  4. L. C. Potter, E. Ertin, J. T. Parker, and M. Çetin, “Sparsity and compressed sensing in radar imaging,” Proceedings of the IEEE, vol. 98, no. 6, pp. 1006–1020, 2010. View at Publisher · View at Google Scholar · View at Scopus
  5. J. Wright, Y. Ma, J. Mairal, G. Sapiro, T. S. Huang, and S. Yan, “Sparse representation for computer vision and pattern recognition,” Proceedings of the IEEE, vol. 98, no. 6, pp. 1031–1044, 2010. View at Publisher · View at Google Scholar · View at Scopus
  6. A. Y. Yang, M. Gastpar, R. Bajcsy, and S. S. Sastry, “Distributed sensor perception via sparse representation,” Proceedings of the IEEE, vol. 98, no. 6, pp. 1077–1088, 2010. View at Publisher · View at Google Scholar · View at Scopus
  7. R. Robucci, J. D. Gray, L. K. Chiu, J. Romberg, and P. Hasler, “Compressive sensing on a CMOS separable-transform image sensor,” Proceedings of the IEEE, vol. 98, no. 6, pp. 1089–1101, 2010. View at Publisher · View at Google Scholar · View at Scopus
  8. E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489–509, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  9. D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM Review, vol. 43, no. 1, pp. 129–159, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  11. A. M. Bruckstein, D. L. Donoho, and M. Elad, “From sparse solutions of systems of equations to sparse modeling of signals and images,” SIAM Review, vol. 51, no. 1, pp. 34–81, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  12. J. A. Tropp and S. J. Wright, “Computational methods for sparse solution of linear inverse problems,” Proceedings of the IEEE, vol. 98, no. 6, pp. 948–958, 2010. View at Publisher · View at Google Scholar · View at Scopus
  13. B. K. Natarajan, “Sparse approximate solutions to linear systems,” SIAM Journal on Computing, vol. 24, no. 2, pp. 227–234, 1995. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  14. R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society, Series B, vol. 58, no. 1, pp. 267–288, 1996. View at Zentralblatt MATH
  15. B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, “Least angle regression,” The Annals of Statistics, vol. 32, no. 2, pp. 407–499, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  16. M. E. Tipping, “Sparse Bayesian learning and the relevance vector machine,” Journal of Machine Learning Research, vol. 1, no. 3, pp. 211–244, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  17. D. P. Wipf and B. D. Rao, “Sparse Bayesian learning for basis selection,” IEEE Transactions on Signal Processing, vol. 52, no. 8, pp. 2153–2164, 2004. View at Publisher · View at Google Scholar
  18. S. Ji, Y. Xue, and L. Carin, “Bayesian compressive sensing,” IEEE Transactions on Signal Processing, vol. 56, no. 6, pp. 2346–2356, 2008. View at Publisher · View at Google Scholar
  19. S. D. Babacan, R. Molina, and A. K. Katsaggelos, “Bayesian compressive sensing using Laplace priors,” IEEE Transactions on Image Processing, vol. 19, no. 1, pp. 53–63, 2010. View at Publisher · View at Google Scholar
  20. S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” IEEE Transactions on Signal Processing, vol. 41, no. 12, pp. 3397–3415, 1993. View at Publisher · View at Google Scholar · View at Scopus
  21. Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, “Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition,” in Proceedings of the 27th Asilomar Conference on Signals, Systems & Computers, vol. 1, pp. 40–44, Pacific Grove, Calif , USA, November 1993. View at Scopus
  22. J. A. Tropp, “Greed is good: algorithmic results for sparse approximation,” IEEE Transactions on Information Theory, vol. 50, no. 10, pp. 2231–2242, 2004. View at Publisher · View at Google Scholar
  23. J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE Transactions on Information Theory, vol. 53, no. 12, pp. 4655–4666, 2007. View at Publisher · View at Google Scholar
  24. W. Dai and O. Milenkovic, “Subspace pursuit for compressive sensing signal reconstruction,” IEEE Transactions on Information Theory, vol. 55, no. 5, pp. 2230–2249, 2009. View at Publisher · View at Google Scholar
  25. D. Needell and J. A. Tropp, “CoSaMP: iterative signal recovery from incomplete and inaccurate samples,” Applied and Computational Harmonic Analysis, vol. 26, no. 3, pp. 301–321, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  26. I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Communications on Pure and Applied Mathematics, vol. 57, no. 11, pp. 1413–1457, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  27. T. Blumensath and M. E. Davies, “Iterative hard thresholding for compressed sensing,” Applied and Computational Harmonic Analysis, vol. 27, no. 3, pp. 265–274, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  28. K. Qiu and A. Dogandzic, “Double overrelaxation thresholding methods for sparse signal reconstruction,” in Proceedings of the 44th Annual Conference on Information Sciences and Systems (CISS '10), Princeton, NJ, USA, March 2010. View at Publisher · View at Google Scholar · View at Scopus
  29. K. Qiu and A. Dogandzic, “Sparse signal reconstruction via ECME hard thresholding,” IEEE Transactions on Signal Processing, vol. 60, no. 9, pp. 4551–4569, 2012.
  30. Y. Wang and W. Yin, “Sparse signal reconstruction via iterative support detection,” SIAM Journal on Imaging Sciences, vol. 3, no. 3, pp. 462–491, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  31. A. Maleki and D. L. Donoho, “Optimally tuned iterative reconstruction algorithms for compressed sensing,” IEEE Journal on Selected Topics in Signal Processing, vol. 4, no. 2, pp. 330–341, 2010. View at Publisher · View at Google Scholar · View at Scopus
  32. T. Blumensath and M. E. Davies, “Gradient pursuits,” IEEE Transactions on Signal Processing, vol. 56, no. 6, pp. 2370–2382, 2008. View at Publisher · View at Google Scholar
  33. B. Sturm, “A study on sparse vector distributions and recovery from compressed sensing,” http://arxiv.org/abs/1103.6246.
  34. X. Zhang, J. Wen, Y. Han, and J. Villasenor, “An improved compressive sensing reconstruction algorithm using linear/non-linear mapping,” in Proceedings of the Information Theory and Applications Workshop (ITA '11), pp. 146–152, San Diego, Calif, USA, February 2011. View at Publisher · View at Google Scholar · View at Scopus
  35. X. Zhang, Z. Chen, J. Wen, J. Ma, Y. Han, and J. Villasenor, “A compressive sensing reconstruction algorithm for trinary and binary sparse signals using pre-mapping,” in Proceedings of the Data Compression Conference (DCC '11), pp. 203–212, Snowbird, Utah, USA, March 2011. View at Publisher · View at Google Scholar · View at Scopus
  36. V. Stodden, L. Carlin, D. Donoho et al., “SparseLab: seeking sparse solutions to linear systems of equations,” Tech. Rep., Stanford University, 2011.