Sparse Signal Recovery via ECME Thresholding Pursuits

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 [1–7]. 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 𝑥₀∈ℝ^𝑛, the sparse signal recovery is intimately related to solving an underdetermined system of linear equations 𝑦=𝐴𝑥₀ 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 ℓ₀ by ℓ₁ for the (𝑃_ℓ0) problem due to nonconvexity and combinatorial effect. In [10], it is the basis pursuit problem (BP):min_𝑥0‖‖𝑥₀‖‖_ℓ1s.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) [21–23], 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 𝐴∈ℝ^𝑚×𝑛. 𝐴^†=𝐴^𝑇(𝐴𝐴^𝑇)⁻¹ for 𝑚≤𝑛; 𝐴^†=(𝐴^𝑇𝐴)⁻¹𝐴^𝑇 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).(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(𝑡)_2(𝑡)=𝑥^(𝑡)+𝛼₁𝑥^(𝑡)−𝑥^(𝑡−1),=𝑥^1(𝑡)+𝛼₂𝑥^1(𝑡)−𝑥^(𝑡−2),(2.2) where 𝛼₁,𝛼₂ 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 𝐼. 𝑥_𝐼=argmin_z{‖𝑦−𝐴_𝐼𝑧‖²²,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_𝑥‖‖𝑥_𝑇‖‖₁s.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,  initialize residual signal 𝑟=𝑦,  initialize support set 𝐼⁽⁰⁾=∅,  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.

Input: Measurement matrix 𝐴, measurements y, sparsity level 𝑘 Output: The reconstructed signal 𝑥 (1)  Initialization: (2)  𝑡=1 //iteration number (3)  𝑥(0)=𝟎 //initial signal (4)  𝑟(0)=𝑦 //initial residual (5)  𝐼(0)=∅ //initial support set (6)  while  halting criterion false do (7) 𝑥(𝑡)=𝐻𝑘(𝑥(𝑡−1)+𝐴†𝑟(𝑡−1)) (8) 𝐼(𝑡)={𝑖∶𝑥𝑖(𝑡)≠0} (9) 𝑥𝐼(𝑡)(𝑡)=𝐴†𝐼(𝑡)𝑦 (10) 𝑥(𝑡)[1,𝑛]⧵𝐼(𝑡)=0 (11) 𝑟(𝑡)=𝑦−𝐴𝑥(𝑡) (12) 𝑡=𝑡+1 (13)  end while (14)  return  𝑥

Algorithm 1 

EMTP-𝑘 algorithm.

Input: Measurement matrix 𝐴, measurements 𝑦, thresholding parameter 𝛽 Output: The reconstructed signal 𝑥 (1)  Initialization: (2)  𝑡=1 //iteration number (3)  𝑥(0)=𝟎 //initial signal (4)  𝑟(0)=𝑦 //initial residual (5)  𝐼(0)=∅ //initial support set (6)  while  halting criterion false do (7) 𝑥(𝑡)=𝑥(𝑡−1)+𝐴†𝑟(𝑡−1) (8) 𝐼(𝑡)={𝑖∶|𝑥𝑖(𝑡)|>𝛽𝑡max|𝑥|} (9) 𝑥𝐼(𝑡)(𝑡)=𝐴†𝐼(𝑡)𝑦 (10) 𝑥(𝑡)[1,𝑛]⧵𝐼(𝑡)=0 (11) 𝑟(𝑡)=𝑦−𝐴𝑥(𝑡) (12) 𝑡=𝑡+1 (13)  end while (14)  return  𝑥

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].

(a) The 1st iteration

(b) The 2nd iteration

(c) The 3rd iteration

(d) The 4th iteration

(e) The 5th iteration

(f) The 6th iteration

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 𝒪(𝑚²𝑛+𝑚𝑛)). 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 𝒪(𝑘²𝑚+𝑘𝑚). The step for updating residual requires an additional matrix-vector operation, costing 𝒪(𝑘𝑚). To summarize, one EMTP iteration costs 𝒪(𝑚𝑛+𝑛log𝑘+𝑘²𝑚+2𝑘𝑚), which is significantly less than the complexity of one DORE step (with the cost of 𝒪(2𝑛log𝑘+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‖𝑦−𝐴𝑥‖²²  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 ‖‖𝑦−𝑟^(𝑡)‖‖₂‖𝑦‖₂≤10⁻⁶.(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)=−20log₁₀‖‖𝑥₀−𝑥‖‖₂‖‖𝑥₀‖‖₂,(4.2) where 𝑥₀ 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.

(a)

(b)

Figure 2 

Influence of thresholding parameter 𝛽 for sparse Gaussian signals with 𝑛=400,𝑘=20: comparisons in terms of SNR (dB) and number of iterations.

(a)

(b)

Figure 3 

Influence of thresholding parameter 𝛽 for sparse Gaussian signals with 𝑛=400, 𝑚=80: comparisons in terms of SNR (dB) and number of iterations.

(a)

(b)

Figure 4 

Influence of thresholding parameter 𝛽 for sparse Gaussian signals with 𝑛=2000, 𝑘=100: comparisons in terms of SNR (dB) and number of iterations.

(a)

(b)

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.

Figure 6 

Number of iterations as a function of problem size with fixed ratios of 𝑚/𝑛=0.4, 𝑘/𝑚=0.2, n = 200 : 400 : 2200.

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 ℓ₁ 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.

(a)

(b)

(c)

Figure 8 

Sparse Gaussian signals with 𝑛=400, 𝑘=20: comparisons in terms of SNR (dB) and number of iterations.

(a)

(b)

Figure 9 

Sparse Gaussian signals with 𝑛=400, 𝑚=80: comparisons in terms of SNR (dB) and number of iterations.

(a)

(b)

Figure 10 

Sparse Gaussian signals with 𝑛=2000, 𝑘=100: comparisons in terms of SNR (dB) and number of iterations.

(a)

(b)

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.

(a)

(b)

Figure 12 

Sparse Laplacian signals with 𝑛=400, 𝑘=20: comparisons in terms of SNR (dB) and number of iterations.

(a)

(b)

Figure 13 

Sparse Laplacian signals with 𝑛=2000, 𝑘=100: comparisons in terms of SNR (dB) and number of iterations.

(a)

(b)

Figure 14 

Power-law decaying signals with 𝑛=400, 𝑘=20: comparisons in terms of SNR (dB) and number of iterations.

(a)

(b)

Figure 15 

Power-law decaying signals with 𝑛=2000, 𝑘=100: comparisons in terms of SNR (dB) and number of iterations.

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.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.

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.