J. Cent. South Univ. (2015) 22: 3946−3956 DOI: 10.1007/s11771-015-2939-2

A new backtracking-based sparsity adaptive algorithm for distributed compressed sensing XU Yong(徐勇)1, 2, 3 , ZHANG Yu-Jie(张玉洁)1, XING Jing(邢婧)1, 2, 3, LI Hong-wei(李宏伟)1, 3 1. School of Mathematics and Physics, China University of Geosciences, Wuhan 430074, China; 2. Institute of Statistics, Hubei University of Economics, Wuhan 430205, China; 3. Hubei Subsurface Multi-scale Imaging Key Laboratory, China University of Geosciences, Wuhan 430074, China © Central South University Press and Springer-Verlag Berlin Heidelberg 2015 Abstract: A new iterative greedy algorithm based on the backtracking technique was proposed for distributed compressed sensing (DCS) problem. The algorithm applies two mechanisms for precise recovery soft thresholding and cutting. It can reconstruct several compressed signals simultaneously even without any prior information of the sparsity, which makes it a potential candidate for many practical applications, but the numbers of non-zero (significant) coefficients of signals are not available. Numerical experiments are conducted to demonstrate the validity and high performance of the proposed algorithm, as compared to other existing strong DCS algorithms. Key words: distributed compressed sensing; sparsiy; backtracking; soft thresholding

1 Introduction Compressed sensing (CS) [1] is a digital data acquisition theory about signal sensing and compression, which has attracted increasing interests in recent years. The major principle of CS is that sparse or compressible signals with a few significant samples in one basis can be reconstructed almost perfectly from a small number of random projections onto a second basis that is incoherent with the first one [2]. In other words, CS means that certain types of signals can be approximated well by a fraction of samples that contain sufficient information. Assume that x is an N-dimension signal, it is called K-sparse (or compressible) if x can be well approximated by K(<
(1)

where y  R M 1 is called measurement vector with M(<
reduction in digital data acquisition [3]. For recovering x exactly, it requires the measurement matrix with low coherence [1−3]. Gaussian matrices have been proved to satisfy the condition of coherence. A natural and straightforward approach for CS is to translate this recovery problem into l0-norm minimization problem. Unfortunately, solving l0-norm minimization problem is NP-hard since it contains an exhaustive search over all possible solutions, which cannot be used for practical applications [4−5]. Therefore, much research in CS has been done to explore faster recovery algorithms with sustainable reconstruction complexity. DONOHO [1] and CANDES et al [4, 6] indicated that a much easier l1-norm optimization based on linear programming (LP) techniques is equivalent to the l0-norm optimization, if the measurement matrix Φ satisfies the so-called restricted isometry property (RIP) with a constant parameter. Although the l1-norm optimization technique can yield nearly optimal recovery result [4−5], its complexity is still intolerable for some large scale applications, for instance, the complexity of LP algorithms based on interior point methods is O(M2N3/2) [7]. In recent years, several recovery algorithms based on iterative greedy pursuit have been researched widely due to their low complexity for signal recovery. The list

Foundation item: Projects(61203287, 61302138, 11126274) supported by the National Natural Science Foundation of China; Project(2013CFB414) supported by Natural Science Foundation of Hubei Province, China; Project(CUGL130247) supported by the Special Fund for Basic Scientific Research of Central Colleges of China University of Geosciences Received date: 2014−09−05; Accepted date: 2015−01−31 Corresponding author: LI Hong-wei, Professor, PhD; Tel: +86−18627905633; E-mail: [email protected]

J. Cent. South Univ. (2015) 22: 3946−3956

includes orthogonal matching pursuit (OMP) [8], the improved backward optimized orthogonal matching pursuit algorithm (IBOOMP) [9] and the regularized OMP (ROMP) [10]. These greedy pursuit algorithms rest on the idea that the column index accounts for the greatest correlation between the measurement matrix Φ and the residual (the original residual is measurement vector) conducts a good estimation of true support set of signal. According to this strategy, greedy pursuit algorithms explore the support of the signal iteratively. At each iteration, one or more coordinates of the signal are selected for identifying, according to the correlation values between the columns of Φ and the residual produced at previous iteration. If the selected coordinates are believed to be sufficiently reliable, the candidate column indices are added to the current estimate of the support set of the signal x. The iterations repeat until all the coordinates in the correct support set are acquired. Computational complexities of these algorithms are significantly lower than those of LP methods, especially when the signal is very sparse [8−10]. Unfortunately, these algorithms have two obvious drawbacks: 1) they need the sparsity of signal as the halting condition of iteration, i.e., sparsity is required as the prior knowledge, which is not available in many practical applications; 2) they all lack backtracking mechanism for the estimated signal support set at each iteration. Once one or more coordinates are identified to be reliable and added to the candidate list, they are not removed from the list until the algorithm ends. It is obvious that the optimal coordinates obtained at completed iteration may not be contained in the true support set of signal. Therefore, the candidates have to be identified with extreme caution. Subsequently, several new greedy algorithms based on backtracking technique have been proposed, such as compressive sampling matching pursuit (CoSaMP) [11] and the subspace pursuit (SP) [12]. At the k-th iteration of these algorithms, the coordinates obtained at previous iteration are incorporated into the current candidates for reevaluating their reliability. If the measurement vector does not lie in the current estimate for the correct spanning space, one refines the estimate by retaining reliable candidates, discarding the unreliable ones while adding the same number of new candidates during each iteration [12]. The criterion of “reliability” is defined in terms of the order statistics of the inner products of the received signal with the columns of Φ and the subspace projection coefficients. Although these algorithms utilize the backtracking technique, they both need the sparsity of x as prior information, which is unrealistic in many practical applications. In order to overcome the dependence on the sparsity, sparsity adaptive matching pursuit (SAMP) was proposed [13]. SAMP not only releases the limitation of sparsity but also keeps

3947

performance comparable with that of SP, CoSaMP or LP. In the past few years, researchers realize that the CS theory which is mainly designed to exploit intra-signal structures at a single sensor is not sufficient for practical applications, such as wireless sensor network (WSN). Therefore, distributed compressed sensing (DCS), new distributed coding algorithm for multi-signal ensembles that exploit both intra- and inter-signal correlation structures, is proposed [14]. The DCS theory rests on a new concept that we term the joint sparsity of a signal ensemble. In DCS scenario, several sensors measure signals that are individually sparse in some basis and also correlate from sensor to sensor. Each sensor independently encodes its signal by projecting it onto another basis (which is incoherent with the first basis) and then transmits just a number of obtained coefficients to a collection point. Under some right conditions, all of the signals can be jointly reconstructed by decoding at the collection point. Some literatures demonstrate the research for DCS. In Ref. [15], two different joint sparse models (JSM) which can be applied in different scenarios are defined to demonstrate various types of connections among these signals. Based on l1-minimization technique, a single linear program is used to jointly recover sparse signals that satisfy these joint sparse models [16]. Simultaneous orthogonal matching pursuit (SOMP) [17] is proposed for the multi-signal joint sensing and reconstruction (termed as DCS-SOMP). Unfortunately, since DCS-SOMP is inspired by OMP technique, it lacks the backtracking mechanism for the estimated signal support sets. Five greedy algorithms for DCS are proposed in Ref. [18], including simultaneous iterative hard thresholding (SIHT), simultaneous normalized IHT (SNIHT), simultaneous hard thresholding pursuit (SHTP), simultaneous normalized HTP (SNHTP), and simultaneous CoSaMP (SCoSaMP). Inspired by SP algorithm, SUNDMAN et al [19] applied it in DCS scenario to jointly recover sparse signals that satisfy these joint sparse models. Unfortunately, sparsity is needed as the prior information, which is not realistic in practical applications [18−19]. Furthermore, the methods mentioned in Refs. [14−17] cannot deal with the DCS problem based on the model in Ref. [20]. The model can be seen as a generalization over all previously proposed models. To overcome these drawbacks above, SAMP for DCS (referred as DCSSAMP) is proposed [20]. DCSSAMP has provable reconstructing accuracy as that of LP optimization technique, and it demonstrates low complexity comparable to that of OMP methods. In this work, a new iterative greedy algorithm is proposed to solve the DCS problem. This algorithm based on backtracking technique can reconstruct several

J. Cent. South Univ. (2015) 22: 3946−3956

3948

input signals jointly without any prior information of sparsity. Following the same model in Ref. [20], common and individual support-sets are used in this work to characterize jointly sparse signals. Unlike Ref. [15], there is no restriction on the associated signal values. The signal model we use is called mixed supportset model. In order to guarantee the performance of the proposed algorithm, the mechanisms of soft thresholding and cutting are applied at each iteration. For highlighting the advantages of our algorithm, it is compared with DCS-SAMP, joint-OMP and joint-SP [19] through several numerical simulations, in terms of the requirement of measurements, the robustness against sparsity and the numbers of signals. In order to test the validity for practical application, a DCS problem of speech signals is solved by the proposed algorithm and DCS-SAMP respectively. Robustness of our algorithm against noise is tested in these experiments. Recovery of image is also performed by the proposed algorithm.

2 Distributed compressed sensing and mixed support-set model 2.1 Distributed compressed sensing In this work, we consider the same DCS model in Ref. [20], which exploits both intra- and inter-signal correlation structures. That is y j  Ax j ( j  1, 2, , L)

(2)

where each signal x j  R N  L . For the signal sparse in certain transformed domain, there exists a determined sparse basis P by which xj can be sparsely represented by no more than K(<
(3) N×L

represents the ensemble where X=(x1, x2, …, xL)  R of K-sparse signals, Y=(y1, y2, …, yL)  RM×L is the set of measurements and E=(e1, e2, …, eL)  RM×L is the noise term. The goal of DCS is to simultaneously reconstruct X from Y.

2.2 Mixed support-set model In this subsection, three existing joint sparsity models (JSM) that exploit the interrelation between the sparse signals in Eq. (2) are introduced. They are applied in different scenarios. In these models, each signal is itself sparse so that the CS algorithms could be applied to encode and decode each one separately. If there exists no

interrelation between xj, it roughly requires c×K×L (c is a constant) measurements in total to fully represent all the signals and L independent reconstruction procedures. In contrast, taking advantages of the joint sparse models, fewer total measurements are required [14, 17−18]. 2.2.1 JSM-1: Common sparse supports model In this model, all signals are generated from the same sparse set of basis vectors with different coefficients [17]: x j  Ψθ j ( j  1, 2, , L)

(4)

where θj is supported only on the same   {1, 2, …, N} with |Ω|=K and K is the sparsity of xj. This means that all signals are K-sparse with the same coordinates but different coefficients. 2.2.2 JSM-2: Sparse common component + innovations In this model, all signals share a common sparse component while each individual signal carries a sparse innovation component [17], i.e., x j  z  z j ( j  1, 2, , L)

where z  Ψθ z , θ z

0

 K z and z j  Ψθ j , θ j

(5) 0

 K j.

The l0-norm θ 0 only counts the number of nonzero entries in vector θ. Obviously, the Kz-sparse signal z is common to all of the x j (j  1, 2, , L). The signals

zj are the unique components of xj with sparsity Kj. 2.2.3 JSM-3: Sparse common support + sparse independent innovation (mixed support-set) This model extends JSM-2 so that there is no restriction on the associated signal values. In other words, the common component is common only in the support set while nonzero coefficients of the innovation parts and common parts are different among the signals [20]. JSM-3 can be represented as x j  c j  z j  Ψ (θ j  σ j ) ( j  1, 2, , L)

(6)

where the vector θj is common among all signals, supported only on the same   {1, 2,  , N } with   Kz, where Kz is the sparsity of θj; σj is the innovation part coefficient vector with the nonzero support cardinality Kj. One can find that Eq. (6) could be seen as a generalization over Eq. (4) and Eq. (5). Assume that each xj is K-sparsity. If the interrelation among them is ignored, it roughly requires c×K×L measurements to fully represent all of the signals and L independent reconstruction procedures. In contrast, applying the mixed support-set model, the DCS method requires only L

c  ( K c   K j ) measurements and one reconstruction j 1

procedure to recover the sparse signals simultaneously [20]. In this work, a DCS algorithm is proposed based on

J. Cent. South Univ. (2015) 22: 3946−3956

JSM-3.

3 Proposed algorithm for DCS 3.1 Algorithm description In this section, the proposed algorithm is discussed in detail. The algorithm is based on backtracking technique, not requiring sparsity as prior information. Moreover, the mechanism of soft thresholding is adapted for preselecting coordinates. For pruning the candidates, the cutting mechanism is performed for the proposed algorithm. For convenience and without loss of generality, each signal xj is assumed to be K-sparse in time domain. The recovery process is divided into several iterations. At each iteration, some coordinates in the true support set which contains the column indices accounting for s (the number s is determined by the soft threshold) greatest amounts of residual energy across all signals are identified (Step 1). Then, the indices are incorporated into the estimate of support set obtained at last iteration, generating a new candidate list (Step 2). Under the right condition, the candidate list is refined by the cutting mechanism and a new residual is built based on the cutting result, otherwise the new residual is constructed directly based on the candidate list (Step 3). The algorithm iterates this procedure until all the coordinates in the real support set are included in the estimated support set. Finally, estimates of the signals can be obtained by solving a least-square problem (Step 4). For clarity, some notions we propose are listed as follows: apr, bpr: Fuzzy threshold parameters;

3949

Ω={i:, (v )i  thr  max(v )}. Step 2: Ck  Ck 1   , I  length(Ck )  1, F= max( ACk Y

row 

Step 3: If r

2

, I ), r  Y  AF AF Y .

 rk 1

2

 pru  Y

2

Ck  F , rk  r . end

rk  Y  ACk ACk Y , k  k  1. end

Step 4: xˆ j (Ck )  ACk y j , xˆ j ({1, 2, , N }  Ck)=0. Inspired by Ref. [21], apr and bpr are used to build soft threshold parameter thr for controlling the number of selected coordinates in Step 1. It is vital to select the values of apr and bpr. Too small values of them may decrease the recovery precision while too large may result in the increase of performing expense. On the other hand, when we consider the stopping factor, too large value may decrease the recovery precision. The cutting factor pru mainly impact the performing time of algorithm. Obviously, there must be a balance between precision and performing expense for all the parameters. 3.2 Theoretical analysis This subsection describes theoretical analysis for the thresholding mechanism in our algorithm. First of all, the concept of restricted isometry property (RIP) is introduced. Definition 1 [12]: A matrix A  R M  N is set to satisfy (RIP) parameters (K, δ) for K≤M, 0≤δ≤1, if for all index sets I  {1, 2, , N } such that |I|≤K and for all I x  R , one has 2 2

2 2

Cutting factor;

(1   ) x

stop:

Stopping factor for iteration;

k:

Iteration index;

A+:

Pseudo-inverse of matrix A;

where |I| is the cardinality of I. We define δK, the RIP constant, as the infimum of all parameters δ for which the RIP holds, i.e.,

AI:

A submatrix of A which is composed of the columns of A with indices i  I ; A vector composed of the entries of vector x with indices i  I ; Residue in the k-th iteration;

x(I): rk: A

row 

: Taking ∞-norm of matrix A row by row;

Max(x, b): Set of indices corresponding to the b largest amplitude components of x. The whole procedure of our algorithm is described as follows: Input: Measurement matrix A, measurement Y, apr, bpr, pru and stop. Initialization: xˆ j =0 (j=1, 2, …, L), r0=Y, k=0, C0=  Step 1: While rk 1 2  stop  Y 2

v  AT rk 1

row 

, thr  apr  (bpr  apr )  rand(1),

 AI x

 (1   ) x

2 2

pru:

 K : inf{ : (1   ) x

2 2

 AI x

2 2

I

 I  K ,x  R }

(7)

2

 (1   ) x 2 ,

(8)

For all index sets Γ and all A for which the RIP holds s≥|Γ|, the following inequalities hold trivially [22]: (1    ) x ( ) A* y

2

2

 A* A x ( )

 (1    ) y

2

2

 (1    ) x ( ) 2 (9) (10)



where x  R N , y  R . Suppose Γ and Λ are two disjoint sets. From Ref. [23], if the RIP holds A with     , then A* A x

2

    x

where x  R



2

(11)

. Another conclusion demonstrated in

J. Cent. South Univ. (2015) 22: 3946−3956

3950

Ref. [23] is for two integers K1 and K2,

From Eq. (15) and Eq. (20), we have

 K1   K 2

(12)

is established if K1≤K2. Following the discussion above, the rationality of the thresholding mechanism of our algorithm can be explored. Denote C as the final list in our algorithm and q=|C|. One may wonder how to ensure the correct indices of the support set of xj (j=1, 2, …, L) being found in the iterations. Corollary 1: Denote a index set Λ={i1, i2, …, is}, let s  N and s≤q, the RIP holds A with    q  s  (1  2) q

, then

2

(21)

2

Following Eq. (12) and Eq. (21), we have

 sq x j

2

 (1   s  q ) x j (C)   s  q x j (C \ C) 2

2

(22)

If the entries of |xj| are arranged in descending order, i.e., x j (n1 )  x j (n2 )    x j (ns )  x j (ns 1 )   

(23)

x j (nq )

sup

(14)

Am , y j

m{1,2,, N }, m  i ,i

where y j  Ax j  AC x j (C ). Proof: Suppose   C  , and from Eq. (11), we have  AC* y j

 AC* AC x j (C )

2

| AC* AC x j (C ) 

2

(15) Denote C  {n1 , n2 , , ns } as the index corresponding to the s largest entries of |xj| and

x j (n1 )  x j (n2 )    x j (ns ) 

sup

n{1,2,, N }, nC

set

x j ( n)

(16) where |xj| is a vector composed of the absolute values of ~ entries of xj. Obviously, C is a subset of C. From Eq. (14), *

2



AC* y j

AC* A C \C

2



AC* AC x j (C )

x j (C \ C)

AC* AC \ C x j (C \ C )

2



2

|

AC* AC x j (C ) 

AC* AC x j (C ) 2

2

(17)

2

 (1   s ) x j (C )

2

  q x j (C \ C )

(19)

2

2

 (1   s ) x j (C)   q x j (C \ C) 2

2

(20)

(25)

~ From the definition of C and Eq. (25), it is obtained: x j (C)

2

s xj q



(26)

2

Note that 2 x j (C)  x j (C \ C) 2

2 2

 x j (C )

2 2

 xj

2

(27)

2

then x j (C \ C )

2

 1

s xj q

(28)

2

Following Eq. (22), Eq. (26) and Eq. (28), we have

 sq x j

2

 (1   s  q )

s xj q

2

  sq 1 

s xj q

2

(29)

This means s s s   sq   sq   sq 1  q q q

(30)

Note that s s  1  2 q q

then *

x 2j (ns )    x 2j (nq ))

(18)

2

On the other hand, from Eq. (11),

AC* AC \C x j (C \ C)

q( x 2j (n1 )  x 2j (n2 )    x 2j (ns ))  s( x 2j (n1 )   



where C \ C is composed of the entries of C that are not ~ in C . From Eq. (9), we get

AC* AC x j (C)

(24)

Actually, there are s(q−s) items on each side of Eq. (24) while items on the left are greater than or equal to that on the right. Based on Eq. (24), it is obtained:

Ai1 , y j  Ai2 , y j    Ais , y j 

2

x 2j (ns  2 )    x 2j (nq ))

(13)

if

A y j

 (1   s ) x j (C)   q x j (C \ C)

(q  s)( x 2j (n1 )  x 2j (n2 )    x 2j (ns ))  s( x 2j (ns 1 ) 

 C  

A y j

2

then

s

A* y j

 sq x j

(31)

so, s  (1  2) s  q q

(32)

J. Cent. South Univ. (2015) 22: 3946−3956

3951

that is

 sq 

s (1  2) q

(33)

Obviously, Eq. (33) is contradictory with the hypothesis in Corollary 1. Hence,   C  . The proof is completed. Corollary 1 reveals that our algorithm based on thresholding mechanism enables at least one correct index in the support set of xj being found in the first iteration. As y j  Ax j  AC x j (C ), yj can be represented as  x j (Ck )  y j  ACk , AC Ck   x j (C  Ck )   





(34)

where Ck is the final list in the k-th iteration. For the following iterations, there are C−Ck correct indices to be found. For yj, the residue in the k-th iteration is approximately equal to AC Ck x j (C  Ck ). In Corollary 1, this residue and C−Ck are used to replace yj and C respectively while the RIP for A is taken with    q 1  1/(1  2) q , which means at least one correct index of the support set of xj can be found in the following iteration. Thus, we ensure our algorithm based on thresholding mechanism can find correct indices in each iteration.

are assumed to be sparse with sparsity Kz and Kj, respectively. Kj is assumed to be equal for all signals and the measurement matrix A is taken as a M×N random Gaussian matrix.

4.1 Performance of proposed algorithm In this subsection, the performance of the proposed algorithm is demonstrated intuitively. 200 trials are run to test the reconstruction ability of the proposed algorithm with N=200, M=90, Kz=15, Kj=2 and L=3. In each trial, different realizations of the measurement matrix A and X are used. In order to evaluate the reconstructions, the errors (E) between xj and corresponding estimates xˆ j (j=1, 2, …, L) are defined in Eq. (35). The larger the errors are, the better the performance is:

E( j) 

xj

2

x j  x j

( j  1, 2, , L)

(35)

2

Figure 1 demonstrates that, for all the three signals, their errors are kept at a high level over all runs. Actually, all the errors are roughly in interval [1×1014, 17×1014], which shows a high performance of our algorithm.

4 Numerical experiments In this section, the performance of the proposed algorithm is compared with that of DCS-SAMP, JointOMP and Joint-SP, in terms of the requirement for measurements, the robustness against sparsity and the numbers of signals. Note that all the four algorithms can be applied in the model JSM-3, however, Joint-SP and Joint-OMP require the sparsity of signals as the prior information. Algorithms proposed in Refs. [17] and [18] are only applicable to JSM-1, so they are not discussed in this section. A DCS problem of speech signals, the experiment for investigating the robustness against noise and demonstration on a sparse image exhibit the application potential of our algorithm in the end of this section. As discussed in previous sections, a balance between precision and performing expense for all parameters of our algorithm should be found. By a large number of tests, a right combination is found: apr=0.8, bpr=1, stop=6×10−4, and pre=10−2. They are fixed in the experiments below. To evaluate the performance of the proposed algorithm, L input signals whose nonzero coefficients are random Gaussian values are generated. All the signals satisfy the condition for JSM-3. The common components cj and independent innovation components zj

Fig. 1 Errors for all signals in each run

4.2 Comparisons among proposed algorithm, DCSSAMP, Joint-OMP and Joint-SP 4.2.1 Requirement for measurements for four algorithms The probability of success recovery vs the number of measurements is investigated. DCS-SAMP, JointOMP, Joint-SP and the proposed algorithm are used respectively to recover signals with N=200, Kz=15, Kj=6, and L=5, while the number of measurements M is increased from 55 to 100. 200 trials are performed to test the stability of the four algorithms against measurements. The probability of success is represented by the frequency of success in these trials. In a trial, “success” defined for each E(j) (j=1, 2, …, 5) is more than 105. The experiment result is shown in Fig. 2. From Fig. 2, for all the four algorithms, the

3952

J. Cent. South Univ. (2015) 22: 3946−3956

Fig. 2 Probability of success vs number of measurements for four algorithms

Fig. 3 Probability of success vs sparsity of common components of signals for four algorithms

probability of success increases with the increase of M. The probability curves of Joint-OMP and Joint-SP are very close to those of our algorithm over all the values of M. No matter Joint-SP, Joint-OMP or our algorithm, the probability curve increases slowly when it has reached 0.9. However, our algorithm only requires 70 measurements to achieve the probability 0.9 while the number is 75 for Joint-SP and Joint-OMP. In case of M >75, Joint-SP and Joint-OMP are slightly better than our algorithm. Among the four algorithms, DCS-SAMP is the worst for its high demand for M. Actually, its success probability roughly stays at zero when M is no more than 75. The largest increase for the probability curve is obtained by increasing M from 75 to 85. For reaching the probability of 0.9, DCS-SAMP requires 90 measurements in this experiment, which is significantly more than the requirements for other algorithms. As discussed above, there is no significant difference between the performance of our algorithm and that of Joint-SP and Joint-OMP. However, their performances are obviously better than those of DCS-SAMP, in terms of the requirements for measurements. 4.2.2 Robustness against sparsity for four algorithms Considering that our algorithm is based on JSM-3, the relationship between increased sparsity and the performance of algorithm is explored. In this experiment, DCS-SAMP, Joint-SP and Joint-OMP are still compared with our algorithm. Firstly, only the sparsity of common components varies for trials to explore the performance of the four algorithms. 200 trials are performed with N=200, M=90, Kj=2 and L=5, while Kz varies from 25 to 33. The experiment result is shown in Fig. 3. Of particular interest is the sparsity level at which the recovery probability drops far below 1, i.e., the critical sparsity which leads to significant errors in the reconstruction algorithms when they are exceeded. The reason why our analysis begins with Kz=25 lies in the fact that all the four algorithms perform perfectly with

Kz≤25. Their significant differences of performance against sparsity are demonstrated when Kz>27. Actually, when the sparsity of common components increases from 27 to 33, the probability of success of Joint-SP and JointOMP drops rapidly. On the other hand, the trials reveal that the critical sparsity of our algorithm and DCSSAMP algorithm by far exceed that of Joint-SP and Joint-OMP in this experiment. The reconstruction capability of our algorithm is comparable to that of DCSSAMP. The reason why Joint-SP and Joint-OMP are sensitive to the sparsity of common components may lie in the fact that the algorithms utilize the strategy of jointly finding the common support-set to recover each signal iteratively [19]. When the sparsity of common support-set comes to a certain level, at which the algorithms cannot successfully find the common support-set, the curves of probability of success drop far below 1. Secondly, what happens when the sparsity of innovation components is increased for the four algorithms is considered. It must be pointed out that, compared with Kz, the whole sparsity of X is more sensitive to Kj. For all signals xj (j=1, 2, …, L), if each xj is sparse with common sparsity K z and innovation L

sparsity Kj, the sparsity for JSM-3 is K z   K j [20]. j 1

Considering Kj is assumed to be equal for all signals, the whole sparsity for JSM-3 is Kz+L×Kj. This means that a small Kj may lead to a large sparsity for JSM-3, which is not conducive to reconstruction. In this experiment, 200 trials are performed with N=200, M=90, Kz=15, and L=5, while Kj varies from 2 to 9. The experiment result is shown in Fig. 4. Thanks to the fixed sparsity of common components, the probability curves for Joint-SP and Joint-OMP decrease slightly in this experiment. The slight decrease is attributed to the increased sparsity of innovation components, which prevents signals from being

J. Cent. South Univ. (2015) 22: 3946−3956

3953

Fig. 4 Probability of success vs sparsity of innovation components of signals for four algorithms

Fig. 5 Probability of success vs number of signals for two algorithms

recovered successfully. The probability curve of our algorithm has similar trend to that of Joint-SP and JointOMP, however the curve of our algorithm is uniformly lower than that of Joint-SP and Joint-OMP. When the sparsity of innovation components varies from 8 to 9, the probability of our algorithm drops from 0.9 to 0.8, which is not an inspiring result. Among the four algorithms, DCS-SAMP is the worst. The probability of DCS-SAMP stays above 0.9 when the sparsity of innovation components varies from 2 to 6. The largest decrease is obtained by increasing the innovation sparsity from 6 to 8, as a result, the probability drops below 0.1. When the innovation sparsity comes to 9, DCS-SAMP is completely ineffective. Following the discussion above, our algorithm is robust against both sparsity of common components and innovation components in this experiment. 4.2.3 Robustness against number of signals for four algorithms What must be pointed out first is that we do not discuss this problem for Joint-SP and Joint-OMP, since the procedures of Joint-SP and Joint-OMP are recovering signals one by one essentially [19]. Therefore, the increased number of signals has no obvious effect on the performances of Joint-SP and Joint-OMP. In this experiment, L varies while other quantities are fixed with N=200, M=100, Kz=15 and Kj=5. 200 trials are performed for observation. Figure 5 demonstrates what we are interested in. From Fig. 5, no matter our algorithm or DCSSAMP, the probability of success decreases with the increase of L. The reason lies in the fact that, for either our algorithms or DCS-SAMP, signals are recovered jointly, instead of one by one. As discussed in last subsection, the sparsity of signals is Kz+L×Kj, and more signals bring the increase of sparsity for joint reconstruction. The probabilities of success for two algorithms are very close when L≤8. As we can see, in the case of L≤8, their probabilities of success can stay

around 0.9. L=8 is the watershed for the two algorithms. For our algorithm, the probability of success decreases slightly with the increase of L. The probability stays above 0.8, even L has reached 11. Unfortunately, when L is increased from 8 to 11, the probability of success of DCS-SAMP drops from 0.88 to 0 rapidly, which is a shocking decrease. The significant decline mainly occurs when L varies from 8 to 9. We can draw a conclusion that our algorithm is more robust than DCS-SAMP. On the other hand, taking into account that Joint-SP and Joint-OMP require sparisty as prior information, our algorithm is the most practical and robust among the four algorithms in this experiment. 4.3 Applications of our algorithm 4.3.1 A DCS problem for speech signals In practical applications, the sparsity of signals is unknown, for example in radar imaging [24]. Furthermore, not all signals are sparse in time domain. For some practical signals that are sparse in transformed domain, our algorithm can also recover them with tolerable errors. In Fig. 6, there are three speech signal fragments of length N=1000, which are compressed through a generated random Gaussian matrix of size 500×1000. Discrete cosine transform (DCT) is applied to obtain sparse representation of each signal (Fig. 7). Considering that these representations are sparse, our algorithm is employed to recover them. The recovery result is shown in Fig. 8. It is obvious that all the sparse representations are reconstructed properly. Finally, the inverse discrete cosine transform (IDCT) is applied to convert the reconstructed data back to time domain, and recover the three original speech signals. The final result is demonstrated in Fig. 9. To evaluate the reconstruction results, 100 trials are performed. Signal-to-reconstruction-noise-ratio (SRNR, Rstr) is used to evaluate the reconstruction errors, which is defined as

J. Cent. South Univ. (2015) 22: 3946−3956

3954

Fig. 6 Three speech signal fragments: (a) Speech signal 1; (b) Speech signal 2; (c) Speech signal 3 2  xj  2 Rstr (j )  E 10lg  x j  xˆ j 

  (j  1, 2, 3) 2  2

Fig. 7 DCT coefficients of oringinal signals: (a) Speech signal 1; (b) Speech signal 2; (c) Speech signal 3

(36)

where the calculation of mathematical expectation is replaced by an average over the number of trials. The larger the SRNR is, the better the performance is. The obtained SRNR values for the original signals are respectively 36.5272, 39.7656 and 36.6348, which are tolerable for recovery in practical applications. Also, DCS-SAMP is tried to recover these signals, but the recovery completely fails. Except that, since the sparsity is not available in practical applications, Joint-SP and Joint-OMP are not considered in this experiment. 4.3.2 Robustness against noise for our algorithm Except the unknown sparsity, the robustness against noise is another factor that determines the applications perspective of algorithms. To evaluate the robustness

Fig. 8 Recovery of DCT coefficients of original signals: (a) Speech signal 1; (b) Speech signal 2; (c) Speech signal 3

J. Cent. South Univ. (2015) 22: 3946−3956

Fig. 9 Recovery of original speech signals: (a) Speech signal 1; (b) Speech signal 2; (c) Speech signal 3

against noise for our algorithm, M-dimensional vector ej (j=1, 2, …, L) is added to each measurement vector yj (j=1, 2, …, L), respectively, where ej is an independent identically distributed white Gaussian noise, and the symbol SNR is applied to characterize the signal to noise rate between yj and additive noise ej. In this experiment, for M=50 and M=100, with N=200, Kz=15, Kj=2 and L=3, the SNR (Rsn) varies from 1 to 35. 100 trials are conducted to investigate the SRNR for each signal. Figure 10 shows that the SRNR curves with M=50 are uniformly lower than those with M=100 for three signals. One can attribute to the fact that small number of measurements is not conducive to recovery. The larger the number of measurements is, the more information is contained for recovery. No matter M=50 or M=100,

Fig. 10 SRNR vs SNR for three signals

3955

SRNR is increased almost linearly with the increase of SNR for each signal. For M=100, when Rsn=15, the SRNR of each signal is very close to 20, which is an acceptable value in noisy environment. However, to obtain SRNR of 20, noise with Rsn=20 is required to be added to measurements for M=50. 4.3.3 Demonstration on a sparse image In order to evaluate the recovery performance of our algorithm in a more realistic case, recovery of the 256×256 image “cameraman” is demonstrated (Fig. 11). The recovery is performed using 8×8 blocks. The purpose of such processing is breaking the recovery problem into a number of smaller and hence simpler problems. The image “cameraman” is first preprocessed such that each 8×8 block is K-sparse in the 2D Haar Wavelet basis B, where K=13, i.e. for each block, only 13 largest magnitude wavelet coefficients are kept. Note that, in this case, the signal is not itself sparse, but has a sparse representation in a basis B. Hence, for each block, the measurement is obtained via the sensing matrix AB. For each block, M(=32) observations are taken, where the entries of A are randomly drawn from the Gaussian distribution with mean 0 and variance 1. Denoting xj and yj as the j-th block and corresponding measurement, the CS model is y j  ABx j  As j

(37)

where s j  Bx j is the decomposition vector of the j-th block. Instead of recovering all decomposition vectors independently, our algorithm is applied to jointly recover them. Typically, these decomposition vectors could be jointly recovered four by four. After recovering these decomposition vectors, all blocks are reconstructed and then combined to obtain the whole image. Figure 12 shows the recovery result that provides a peak signalto-noise ratio (PSNR) value of 28.3. By comparing Fig. 12 with Fig. 11, the recovery performance of our algorithm is acceptable.

Fig. 11 Test image cameraman

J. Cent. South Univ. (2015) 22: 3946−3956

3956 [7]

[8]

[9]

[10]

[11]

[12]

Fig. 12 Recovery result of cameraman by our algorithm (Rstr=28.3)

[13]

5 Conclusions [14]

A new DCS algorithm with backing tracking strategy is derived for joint signal recovery. The algorithm utilizes mechanisms such as soft thresholding and cutting for precise recovery. Extensive simulations confirm its validity and high performance, as compared with other three existing strong DCS algorithms. The highlights of the algorithm are its capability of recovering several signals simultaneously without sparsity as prior information and its stability against the number of measurements and the sparsity of signals. The capability of solving DCS problem with unknown sparsity and the robustness against noise exhibit its application perspectives. Further research is required to develop a computationally efficient version of the proposed algorithm in noisy environment as well as to better analytically characterize the performance of the algorithm.

[18]

References

[19]

[1] [2]

[3]

[4]

[5]

[6]

DONOHO D L. Compressed sensing [J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289−1306. LIU Jing, HAN Chong-zhao, HU Yu. A novel compressed sensing based method for reconstructing sparse signal using space-time data in airborne radar system [C]// Chinese Control Conference. Xi’an, China, 2013: 4826−4831. RAUHUT H, SCHNASS K, VANDERGHEYNST P. Compressed sensing and redundant dictionaries [J]. IEEE Transactions on Information Theory, 2008, 54(5): 2210−2219. CANDES E J, ROMBERG J, TAO T. Stable signal recovery from incomplete and inaccurate measurements [J]. Communications on Pure and Applied Mathematics, 2006, 59(8): 1207−1223. CANDES E J, ROMBERG J, TAO T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information [J]. IEEE Transactions on Information Theory, 2006, 52(2): 489−509. CANDES E J, TAO T. Decoding by linear programming [J]. IEEE Transactions on Information Theory, 2005, 51(12): 4203−4215.

[15]

[16]

[17]

[20]

[21]

[22]

[23]

[24]

NESTEROV Y, NEMIROVSKII A. Interior-point polynomial algorithms in convex programming [M]. Philadelphia: SIAM, 1994: 20−89. TROPP J A, GILBERT A C. Signal recovery from random measurements via orthogonal matching pursuit [J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655−4666. FANG Hong, ZHANG Quan-bing, WEI Sui. Image reconstruction based on improved backward optimized orthogonal matching pursuit algorithm [J]. Journal of South China University of Technology: Natural Science, 2008, 36(8): 23−27. (in Chinese) NEEDELL D, VERSHYNIN R. Signal recovery from inaccurate and incomplete measurements via regularized orthogonal matching pursuit [J]. IEEE Journal of Selected Topics in Signal Process, 2010, 4(2): 310−316. NEEDELL D, TROPP J A. CoSaMP: Iterative signal recovery from incomplete and inaccurate samples [J]. Applied and Computational Harmonic Analysis, 2008, 26(3): 301−321. DAI Wei, MILENKOVIC O. Subspace pursuit for compressive sensing signal reconstruction [J]. IEEE Transactions on Information Theory, 2009, 55(5): 2230−2249. DO T T, LU G, NGUYEN N, TRAN T D. Sparsity adaptive matching pursuit algorithm for practical compressed sensing [C]// 42nd Asilomar Conference on Signals, Systems and Computers. Pacific Grove, United States, 2008: 581−587. BARON D, DUARTE M F, WAKIN M B, SARVOTHAM S, BARANIUK R G. Distributed compressed sensing [J]. IEEE Transaction on Information Theory. http://dsp.rice.edu/publications/ distributed-compressive-sensing. DUARTE M F, SARVOTHAM S, WAKIN M B, BARON D, BARANIUK R G. Joint sparsity models for distributed compressed sensing [C]// Online Proceedings of the Workshop on Signal Processing with Adaptive Sparse Structured Representations. Rennes, France, 2005. http://www.researchgate.net/publication/220043736_ Joint_sparsity_models_for_distributed_compressed_sensing. BARON D, DUARTE M F, SARVOTHAM S, WAKIN M B, BARANIUK R G. An information-theoretic approach to distributed compressed sensing [C]// 43rd Allerton Conference on Communication, Control, and Computing. Monticello, United States, 2005. http://dsp.rice.edu/publications/information-theoretic-approachdistributed-compressed-sensing. DUARTE M F, SARVOTHAM S, BARON D, WAKIN M B, BARANIUK R G. Distributed compressed sensing of jointly sparse signals [C]// 39th Asilomar Conference on Signals, Systems and Somputers. Pacific Grove, United States, 2005: 1537−1541. BLANCHARD J, CERMAK M, HANLE D, JING Yi-rong. Greedy algorithms for joint sparse recovery [J]. IEEE Transactions on Signal Processing, 2014, 62(7): 1694−1704. SUNDMAN D, CHATTERJEE S, SKOGLUND M. Greedy pursuits for compressed sensing of jointly sparse signals [C]// European Signal Processing Conference. Barcelona, Spain, 2011: 368−372. WANG Qun, LIU Zhi-wen. A robust and efficient algorithm for distributed compressed sensing [J]. Computers and Electrical Engineering, 2011, 37(6): 916−926. GAN Wei, XU Lu-ping, ZHANG Hua, SU Zhe. Greedy adaptive recovery-algorithm for compressed sensing [J]. Journal of Xidian University, 2012, 39(3): 54−62. BLUMENSATH T, DAVIES M E. Iterative hard thresholding for compressed sensing [J]. Applied and Computational Harmonic Analysis, 2009, 27(3): 265−274. DAI Wei, OLGICA M. Subspace pursuit for compressive sensing: Closing the gap between performance and complexity [EB/OL]. [2008−05−10]. http: // www.dsp. rice.edu/cs. LIU Zhen, WEI Xi-zhang, LI Xiang. Low sidelobe robust imaging in random frequency-hopping wideband radar based on compressed sensing [J]. Journal of Central South University, 2013, 20(3): 702−714. (Edited by FANG Jing-hua)