- Email: [email protected]
Similar sensing matrix pursuit: An efficient reconstruction algorithm to cope with deterministic sensing matrix
Signal Processing 95 (2014) 101–110
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Similar sensing matrix pursuit: An efficient reconstruction algorithm to cope with deterministic sensing matrix Jing Liu a,n, Mahendra Mallick b, ChongZhao Han a, XiangHua Yao a, Feng Lian a a b
School of Electronics and Information Engineering, Xi'an Jiaotong University, Xi'an, PR China Anacortes, WA, USA
a r t i c l e i n f o
abstract
Article history: Received 24 September 2012 Received in revised form 1 August 2013 Accepted 19 August 2013 Available online 8 September 2013
Deterministic sensing matrices are useful, because in practice, the sampler has to be a deterministic matrix. It is quite challenging to design a deterministic sensing matrix with low coherence. In this paper, we consider a more general condition, when the deterministic sensing matrix has high coherence and does not satisfy the restricted isometry property (RIP). A novel algorithm, called the similar sensing matrix pursuit (SSMP), is proposed to reconstruct a K-sparse signal, based on the original deterministic sensing matrix. The proposed algorithm consists of off-line and online processing. The goal of the off-line processing is to construct a similar compact sensing matrix containing as much information as possible from the original sensing matrix. The similar compact sensing matrix has low coherence, which guarantees a perfect reconstruction of the sparse vector with high probability. The online processing begins when measurements arrive, and consists of rough and refined estimation processes. Results from our simulation show that the proposed algorithm obtains much better performance while coping with a deterministic sensing matrix with high coherence compared with the subspace pursuit (SP) and basis pursuit (BP) algorithms. & 2013 Elsevier B.V. All rights reserved.
Keywords: Compressed sensing Deterministic sensing matrix Reconstruction algorithm Similar sensing matrix pursuit algorithm
1. Introduction Compressed sensing has received considerable attention recently, and has been applied successfully in diverse fields, e.g. image processing [1,2], underwater acoustic communication [3], wireless communication [4] and radar systems [5–9]. The central goal of compressed sensing is to capture attributes of a signal using very few measurements. In most work to date, this broader objective is exemplified by the important special case in which a K-sparse vector x A RN (with N large) is to be reconstructed from a small number M of linear measurements with K oM oN. K-sparse signals are the signals that can be represented by K significant coefficients over an N-dimensional basis. This can be compactly described via y ¼ Φx þe: n
ð1Þ
Corresponding author. E-mail address: [email protected] (J. Liu).
0165-1684/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sigpro.2013.08.009
Here, y A RM denotes a measurement vector, Φ represents an M N sensing matrix, and e is an M 1 noise vector. The two fundamental questions in compressed sensing are as follows: how to construct suitable sensing matrix Φ, and how to recover the K-sparse vector x from the measurement vector y efficiently. Tables A1 and A2 in Appendix A list notations for variables used in the paper. In early work of compressed sensing, the entries of the sensing matrix Φ are generated by an independent and identically distributed (i.i.d) Gaussian or Bernoulli process, or from random Fourier ensembles [10–12]. In general, the exact solution to the above second question is shown to be an NP-hard problem [13,14]. If the number of samples (M) exceeds the lower bound of M 4OðK log ðN=KÞÞ, l1 minimization (e.g. BP algorithm) can be performed instead of the exact l0 minimization with the same solution for almost all the possible inputs [14]. An alternative approach to sparse signal recovery is based on the idea of iterative greedy pursuit, and tries to approximate the solution to
102
J. Liu et al. / Signal Processing 95 (2014) 101–110
l0 minimization directly. The greedy algorithms include matching pursuit (MP) [15], orthogonal matching pursuit (OMP) [16], regularized OMP (ROMP) [17], stagewise OMP (StOMP) [18], SP [19], compressive sampling matching pursuit (CoSaMP) [20] and backtracking-based matching pursuit (BAOMP) [21], etc. The reconstruction complexity of these approximate algorithms is significantly lower than that of BP algorithm. In compressed sensing, one of the well-studied conditions on the sensing matrix, which guarantees stable recovery for a number of reconstruction algorithms, is the RIP [13,14]. If a sensing matrix Φ whose column vectors have unit norm and satisfies ð1δK Þ‖x‖22 r‖Φx‖22 rð1 þ δK Þ‖x‖22
ð2Þ
for all possible K-sparse vectors with restricted isometry constant (RIC) δK , then Φ is said to obey K-RIP with δK . The RIP with suitable constant δK guarantees perfect reconstruction [13,14], but it is very hard to check whether a sensing matrix satisfies RIP or not. Coherence, the maximal correlation between two columns in a sensing matrix, is also a well-known performance measure for sensing matrices. For a matrix Φ with columns φ1 ; φ2 ; …; φN , the coherence of Φ is defined as μðΦÞ ¼
max
1 r i;j r N and
jφTi φj j : i a j J φ i J J φj J
ð3Þ
Coherence plays a central role in the sensing matrix construction, because small coherence implies the RIP [22]. In this paper, we are interested in deterministic sensing matrices. Deterministic sensing matrices are useful because in practice, the sampler has to be a deterministic matrix. Although random matrices perform quite well on the average, there is no guarantee that a specific realization works. For the deterministic approaches, the Vandermond matrices seem to be good options, since any K columns of a K N Vandermond matrix are linearly independent. However, when N increases, the constant δK rapidly approaches 1 and some of the K K submatrices become ill-conditioned [23]. A connection between the coding theory and sensing matrices is established in [24] where second order Reed– Muller codes are used to construct bipolar matrices. However, they lack a guarantee on the RIP order. In [25], the authors propose a series of deterministic sensing matrices, the binary, bipolar, and ternary compressed sensing matrices which satisfy the RIP condition. The key concept of coherence is extended to pairs of orthonormal bases. This enables a new choice of the sensing matrices: one simply selects an orthonormal basis that is incoherent with the sparsity basis, and obtains measurements by selecting a subset of the coefficients of the signal in the chosen basis [26]. This approach has successful applications in radar systems [9,27], where an additional sensing matrix H is introduced and the received signal is compressed further by making nonadaptive, linear projections of the direct data sampled at the Nyquist frequency. However, neither of these algorithms mention the hardware implementation of the additional sensing matrix, which is very complex and expensive. In practice, it is challenging to design a deterministic sensing matrix having low coherence. In this paper, we
consider a more general condition when the deterministic sensing matrix has high coherence and does not satisfy the RIP condition. A novel algorithm, called the SSMP algorithm, is proposed to reconstruct the K-sparse signal based on the original deterministic sensing matrix. The proposed algorithm consists of two parts: the off-line processing and the online processing. The goal of the off-line processing is to construct a similar compact sensing matrix with low coherence, which contains as much information as possible from the original sensing matrix. The online processing begins when the measurements arrive, which consists of a rough estimation process and a refined estimation process. In the rough estimation process, an SP algorithm is used to find a rough estimate of the true support set, which contains the indices of the columns that contribute to the original sparse vector. Three kinds of structures of the estimated support set are considered, and three individual refined estimation processes are carried out under these three conditions. We observe from simulation results that the proposed algorithm obtains much better performance when coping with the deterministic sensing matrix with high coherence compared with the SP and BP algorithms. The paper is organized as follows. Section 2 introduces the proposed similar sensing matrix pursuit algorithm. Section 3 presents simulation results, and Section 4 summarizes conclusions. 2. Similar sensing matrix pursuit (SSMP) algorithm Recently, algorithms used to cope with deterministic sensing matrices focus on designing a sensing matrix which satisfies the RIP condition (or with low coherence). However, in practice it is challenging to design a deterministic sensing matrix having a very small restricted isometry constant (coherence). In this paper, we consider a more general condition when a deterministic sensing matrix has high coherence and does not satisfy the RIP condition. We concentrate in developing a novel reconstruction algorithm rather than building a deterministic sensing matrix with low coherence. A novel algorithm called the SSMP is proposed to reconstruct the K-sparse signal, based on the original deterministic sensing matrix. This section introduces the proposed SSMP algorithm. First, the key component of the proposed algorithm, the similar compact sensing matrix, is introduced in Section 2.1, and then the complete algorithm is described in Section 2.2. The complexity analysis of the proposed algorithm is presented in Section 2.3. 2.1. Construction of the similar compact sensing matrix The construction process of the similar compact sensing matrix is based on the similarity analysis of the original sensing matrix. In this paper, similarity is defined as the absolute and normalized inner product between any two different columns of the original sensing matrix Φ: λðφi ; φj Þ ¼
jφTi φj j ; J φi J J φj J
1 ri; j r N
and
i aj:
ð4Þ
J. Liu et al. / Signal Processing 95 (2014) 101–110
It can be seen that coherence is the largest similarity among the columns of a matrix. Any two columns with high similarity are coherent with each other, and vice versa. Therefore, the similarity can be used to distinguish the coherent columns from the incoherent columns of the original sensing matrix. A different way to understand similarity and coherence is by considering the Gram matrix G which is defined as ~ T Φ; ~ G¼Φ
ð5Þ
~ is the normalized sensing matrix obtained from where Φ the original sensing matrix with each column normalized. The off-diagonal entries in G are the similarity values defined in (4). The coherence is the off-diagonal entry with the largest magnitude. A threshold T1 is set properly to distinguish the highly coherent columns from the incoherent columns based on similarity. For each calculated similarity value, fλðφi ; φj Þ; i ¼ 1; …; N; j ¼ 1; …; N; ia jg, if λðφi ; φj Þ is greater than T1, the columns i and j are added to the set of highly coherent columns. The remaining columns that do not belong to the set of highly coherent columns form the set of incoherent columns. Fig. 1 shows the classification results in a similarity plane, where the incoherent columns and the highly coherent columns are indicated by black circles and black squares, respectively. In the similarity plane, the distance between any two columns ðφi ; φj Þ is defined as the similarity distance dsimilar, which is inversely proportional to the similarity value λðφi ; φj Þ, as follows: dsimilar ðφi ; φj Þ ¼ K s =λðφi ; φj Þ;
ð6Þ
where Ks is a constant. Therefore the closer the two columns located in the similarity plane, the more similar they are (with larger similarity). Fig. 1 also shows that any two highly coherent columns may not be close to each other (coherent with each other) in the similarity plane. The highly coherent columns are then further divided into a set of similar column groups (indicated by large circles in Fig. 1) to ensure that the highly coherent columns in a similar column group are similar with each other. Moreover, the similar column groups have the following property. Property 1. For any two similar column groups, which are defined as Γ1 ¼ fα1 ; …; αN g and Γ2 ¼ fβ1 ; …; βM g, the similarity between any two columns from these two different similar Highly coherent column
dsimilar Similar column group
dsimilar
dsimilar Incoherent column
Fig. 1. Classification results in a similarity plane.
103
column groups, e.g. αi A Γ1 and βj A Γ2 , is no larger than T1, as λðαi ; βj Þ rT 1 ;
αi A Γ1
and
βj A Γ2 ;
ð7Þ
while the similarity between any two columns in the same similar column group (e.g. Γ1 ) is larger than a threshold T2 as λðαi ; αj Þ 4T 2 ;
α i ; α j A Γ1
and
ð8Þ
ia j:
Property 1 stipulates that any column in one specific similar column group is highly coherent with other columns inside the same group, while incoherent with any column outside the group (including the highly coherent columns in other similar column groups, and incoherent columns). This is shown in Fig. 1 that any column in one large circle (similar column group) is very close to other columns inside the circle, while far from any column outside the circle. Therefore, any column in one specific similar column group represents the characteristics of other columns in the same group. In this paper, we consider reducing the original sensing matrix to a similar compact sensing matrix by drawing a column from each similar column group, while keeping the incoherent columns unchanged in the new sensing matrix. The newly built similar compact sensing matrix contains as much information as possible from the original sensing matrix. Next, the construction process of the similar compact sensing matrix is briefly introduced. Construction of the similar compact sensing matrix: The set of highly coherent columns and the set of incoherent columns are denoted by Shc ¼ fι1 ; ι2 ; ⋯; ιNhc g and Sic ¼ fζ 1 ; ζ 2 ; ⋯; ζ Nic g, respectively. Nhc is the number of highly coherent columns, and Nic is the number of incoherent columns with Nhc þN ic ¼ N. The set of highly coherent columns Shc is further divided into D similar column groups, fΓ1 ; Γ2 ; …; Γi ; …; ΓD g. We assume that M r D oN. Each similar column group contains more than one highly coherent columns, e.g. Γi ¼ fγi1 ; …; γiN i g, where N Γi indiΓ
cates the number of columns in Γi . Each similar column group is compacted to a condensed column. Here we just select a column from each similar column group randomly considering that the columns in one group are highly coherent and very similar to each other. Finally, we obtain a similar compact sensing matrix by combining the condensed columns and the incoherent columns, as Ψ ¼ 1 2 N ic ½γ1C ; γ 2C ; …; γ iC ; …; γD ,where γiC denotes a conC ; ζ ; ζ ; …; ζ
densed column from Γi . The similar compact sensing matrix Ψ is of size M ðD þ Nic Þ. We have the following propositions for the similar compact sensing matrix. Proposition 1. The coherence of the similar compact sensing matrix is less than or equal to T1. Proof. According to Property 1, the similarity between any two condensed columns is no larger than T1. Moreover, from the division process of the coherent columns and the incoherent columns, the similarity between any two incoherent columns, or between a condensed column and a incoherent column, is no larger than T1. As a result, the similarity between any two columns of the similar compact sensing matrix is no larger than T1. Therefore, the coherence
104
J. Liu et al. / Signal Processing 95 (2014) 101–110
of the similar compact sensing matrix, which is the largest similarity, is less than or equal to T1. □ Proposition 2 (Duarte and Eldar [28], Donoho et al. [29]). Let the signal x be a K-sparse vector and write y ¼ Φ x þe. Denote γ as the noise magnitude ðγ ¼ ‖e‖2 Þ, and mr as the measurement residual in the OMP algorithm. Suppose that K r ð1=μðΦÞ þ1Þ=4 and ε Z γ, where ε is an appropriately chosen bound on the noise magnitude. We have x^ ¼ arg min‖x‖1 x A RN
subject to ‖yΦx‖2 rε:
ð9Þ
Then the estimate x^ in (9) has error bounded by γ þε ^ 2 r pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; ‖xx‖ 1μðΦÞð4K1Þ
ð10Þ
where μðΦÞ denotes the coherence of the sensing matrix Φ. The estimate x^ of the OMP algorithm with halting criterion ‖mr‖2 rγ has the error bound γ ^ 2 r pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; ‖xx‖ 1μðΦÞð4K1Þ
ð11Þ
provided that γ rAð1μðΦÞð2K1ÞÞ=2 for OMP, with A being a positive lower bound on the magnitude of the nonzero entries of x. Proposition 3 (Setting of the threshold T1). In order to guarantee a perfect reconstruction of the sparse vector with large probability, the threshold T1 should satisfy (12) according to Propositions 1 and 2: T 1 r1=ð4K1Þ:
ð12Þ
In this paper the threshold T1 is set as T 1 ¼ 1=ð4K1Þ. However, a small T1 will result in a large number of highly coherent columns. While the restricted selection of T1 is true from a worst-case standpoint, it turns out that the coherence as defined previously does not do justice to the actual behavior of sparse representations and pursuit algorithms' performance. Thus, if we relax our expectations and allow a small fraction of failed reconstructions, then values which are substantially above the bound still lead to successful compressed sensing [30]. In this work, the threshold T1 is set to 0.4 in the simulation. The threshold T2 should be set to a large value to ensure that the columns in a similar column group are highly coherent (similar) with each other. The larger the T2 is, the more similar the columns in a similar column group will be, and vice versa. If T2 is set to a small value, the condensed column, which is selected from a similar column group randomly, cannot represent the characteristics of other columns in the same group. Here the threshold T2 is set to 0.9. Therefore, we can see that the proposed SSMP algorithm is not sensitive to the choice of the threshold values. We only need to set T1 to a small value to ensure the low coherence of the similar compact sensing matrix, and set T2 to a large value to ensure that the columns in a similar column group are highly coherent (similar) with each other.
2.2. The similar sensing matrix pursuit algorithm In this paper, a new algorithm called the similar sensing matrix pursuit, is used to cope with the problem induced by the highly coherent columns. Firstly, the original sensing matrix is transformed to a similar compact sensing matrix based on similarity analysis. The similar compact sensing matrix has low coherence, which guarantees a perfect reconstruction of the sparse vector with a large probability according to Propositions 1–3. An OMP-type algorithm (SP algorithm) is then used to find a rough estimate of the true support set, based on the similar compact sensing matrix and the measurement vector, which contains the indices of the columns contributing to the original sparse vector. The contributing columns may include the incoherent columns and the condensed columns. Since a condensed column is selected from a similar column group randomly, we can only obtain the indices of the contributing similar column groups. From the rough estimation process, we cannot decide the real contributing columns in the contributing similar column groups. Next, a refined estimation process is adopted to provide equal opportunities to all the columns in the contributing similar column groups. All the combinations of the columns in the contributing similar column groups and the contributing incoherent columns are listed, with each combination forming a candidate support set. Accordingly, we can obtain a candidate estimate of the original sparse vector using the least squares algorithm based on each candidate support set [16,19]. Finally, we can find the estimated sparse vector matching the residual best. The proposed algorithm consists of off-line processing and online processing. The off-line processing majors in the construction of the similar compact sensing matrix. The online processing consists of a rough estimation process and a refined estimation process. The main sparse reconstruction is carried out in the rough estimation process, where the SP algorithm is used to find a rough estimate of the true support set. The refined estimation process is a combinational searching process among the candidate support sets for the one matching the true support set best. The estimated support set may only include the indices of the contributing condensed columns or the contributing incoherent columns, besides containing both of them. Three individual refined estimation processes are carried out under the above three conditions. Under the first condition, the true nonzero entries of the original sparse vector fall in a number of similar column groups. All the combinations of the columns in the contributing similar column groups are listed with each combination forming a candidate support set. A candidate estimate of the original sparse vector is obtained based on each candidate support set. We can then find the estimated sparse vector matching the residual best. Under the second condition, when the contributing columns include the incoherent columns only, the nonzero entries are easily identified, since the compressed sensing algorithm can clearly identify the entries corresponding to the incoherent columns of the sensing matrix. Finally, the combination of the two algorithms from the first and the second conditions can be used to estimate the original sparse vector under the third condition. The detailed procedure of the proposed algorithm is presented next.
J. Liu et al. / Signal Processing 95 (2014) 101–110
Algorithm 1. Similar sensing matrix pursuit algorithm. Input: Sensing matrix Φ, measurement vector y Output: The estimated signal x^ O 1. Construction of the similar compact sensing matrix: The process is same with that in Section 2.1. 2. Rough estimation: The SP algorithm is used to find a rough estimate of the true support set based on the measurement vector y and the similar compact sensing matrix Ψ. The estimated support set is represented as 1 2 K a^ ¼ fa^ ; a^ ; …; a^ g. The estimated support set a^ contains the indices of the condensed columns or the indices of the incoherent columns, or both, which contribute to the original sparse vector. 3. Refined estimation: Three individual refined estimation processes are carried out under the following three conditions: (a) a^ contains the indices of the condensed columns only: (i) The estimated support set a^ is represented 1 2 K 1 as a^ ¼ fa^ ðγ iC Þ; a^ ðγ jC Þ; …; a^ ðγ lC Þg, where a^ ðγ iC Þ indicates that the first element of a^ corresponds to γ iC , the index of the ith condensed column. Accordingly, we can obtain a set b^ containing the indices of K similar column ^ as b^ ¼ fΓ i ; Γ j ⋯; Γ l g. groups corresponding to a, The indices of all the columns in the selected similar column groups in b^ are listed and form a final set f^ , where f^ ¼ fγ i1 ; …; γ iN i ; …; γ j1 ; …; Γ
γ jN j ; …; γ l1 ; …; γ lN l g. We assume that the total Γ
Γ
number of the indices in f^ is Hcc. combinations based on the column indices in f^ , and each combination forms a candidate support set, e.g. the pth candidate support set is represented as ϒp ¼ fΥ 1p ðγ i1 Þ; Υ 2p ðγ j1 Þ; …; Υ Kp ðγ l1 Þg; p ¼ 1; 2; …; Nco , where Υ 1p ðγ i1 Þ indicates the first element of ϒp corresponding to γ i1 , and Nco indicates the total number of combinations. The proposed algorithm then solves a least squares problem to approximate the nonzero entries of the original sparse vector on each candidate support set ðϒp ; p ¼ 1; 2; …; Nco Þ, and sets other entries to zero, resulting in an estimate of the original p sparse vector, x^ [16,19]. We can obtain the p estimate x^ as p x^ ϒp ¼ Φ†ϒp y;
ð13Þ
p x^ f1;2;…;Ngϒp ¼ 0;
ð14Þ
where † indicates the pseudo-inverse operation. The matrix Φϒp consists of the columns of p Φ with indices i A ϒp , x^ ϒp is composed of the p
entries of x^ indexed by i A ϒp , and x^ f1;2;…;Ngϒp is composed of the entries of x^ iA f1; 2; …; Ngϒp [19].
p
(iii) Final estimate: Among the obtained estimates 1 2 p N fx^ ; x^ ; …; x^ ; f⋯; x^ co g, find the estimate with the least residual, x^ min , and set x^ min as the output estimated signal: x^ O ¼ x^ min : ð15Þ (b) a^ contains the indices of the incoherent columns only: (i) The estimated support set a^ can be repreT
1 2 K sented as a^ ¼ fa^ ðζ i Þ; a^ ðζ j Þ; …; a^ ðζ l Þg , where 1 i a^ ðζ Þ indicates that the first element of a^
corresponds to ζ i , the ith incoherent column. The estimated support set a^ equals the true support set of the original K-sparse vector due to the incoherence between the contributing incoherent columns according to Propositions 1–3. (ii) Estimation: We can obtain the estimate of the original sparse vector using least squares algorithm [16,19] based on the true support set a^ as x^ a^ ¼ Φ†a^ y;
ð16Þ
x^ f1;2;…;Nga^ ¼ 0:
ð17Þ
(iii) Final estimate: x^ O ¼ x^
ð18Þ
(c) a^ contains the indices of both the condensed columns and the incoherent columns: (i) The estimated support set a^ is represented as 1 V V þ1 k K a^ ¼ fa^ ðγ iC Þ; …; a^ ðγ jC Þ; a^ ðζ Þ; …; a^ ðζ l Þg, 1
(ii) Refined estimation: List C KH cc
p
105
indexed by
V
where fa^ ðγ iC Þ; …; a^ ðγ jC Þg correspond selected condensed columns, and fa^
to
V þ1
V
ðζ k Þ;
K
⋯; a^ ðζ l Þg correspond to KV selected incoherent columns. Accordingly, we can obtain a set b^ containing the indices of V similar column 1 V groups corresponding to fa^ ðγ iC Þ; ⋯; a^ ðγ jC Þg, as b^ ¼ fΓ i ; ⋯; Γ j g. The indices of all the columns in
the selected similar column groups in b^ are listed and form a final set f^ , where f^ ¼ fγ i ; …; 1
γ iN i ; …; γ j1 ; …; γ jN j g. We assume that the total Γ
Γ
number of the indices in f^ is Hci. (ii) Refined estimation: List C VH ci
combinations based on the column indices in f^ , and the pth combination can be represented as rp ¼ fγ i1 ; …;
γ j1 g; p ¼ 1; 2; …; N cb , where Ncb is the total number of combinations. Each combination together with all the selected incoherent columns, fζ k ; …; ζ l g, forms a candidate support set. The pth candidate support set is represented as ϒp ¼ fΥ 1p ðγ i1 Þ; …; Υ Vp ðγ j1 Þ; Υ Vp þ 1 ðζ k Þ; …; Υ Kp ðζ l Þg; p ¼ 1; 2; …; N cb . The proposed algorithm then solves a least squares problem to approximate the nonzero entries of the original sparse vector on each candidate support set ðϒp ; p ¼ 1; 2; …; Ncb Þ, and sets other entries to
106
J. Liu et al. / Signal Processing 95 (2014) 101–110
zero, resulting in an estimate of the original p sparse vector, x^ [16,19]. We can obtain the p estimate x^ as p x^ ϒp ¼ Φ†ϒp y;
ð19Þ
p x^ f1;2;…;Ngϒp ¼ 0:
ð20Þ
1 2 (iii) Final estimate: Among the estimates fx^ ; x^ ; …; p N x^ ; …; x^ cb g, find the estimate with the least residual, x^ min , and set x^ min as the output estimated signal: x^ O ¼ x^ min : ð21Þ
2.3. Complexity analysis In this section we will analyze the computational complexity of the proposed similar sensing matrix pursuit algorithm, and further make a comparison with the SP and BP algorithms based on the computational complexity. The proposed approach consists of off-line processing and online processing. The off-line processing transforms the M N original sensing matrix to a M N′ðN′ oNÞ similar compact sensing matrix based on similarity analysis. The computational complexity of the off-line processing mainly focuses on the computation of the similarity between any two columns of the original sensing matrix, which is of the order of OðððN1ÞN=2ÞM). The online processing consists of the rough and refined estimation processes. In the rough estimation process, an SP algorithm is used to find a rough estimate of the true support set. The computational complexity of the rough estimation process is same as that of the SP algorithm and is of the order of OðMN′Þ according to [20]. Three individual refined estimation processes are then carried out under three different conditions mentioned in Section 2.2. The complexity analysis will be based on these three conditions: 1. a^ contains the indices of the condensed columns only. ^ gives The rough estimate of the true support set, a, correct positions of the nonzero entries corresponding to the contributing condensed columns according to Propositions 1–3 They are the indices of K similar column groups contributing to the original sparse vector. In the refined estimation process, C KHcc combinations are listed out as the candidate support sets, where Hcc indicates the total number of columns in the K contributing similar column groups. Assuming Cmax as the size of the largest similar column group, i.e. the number of columns contained in the largest similar column group, we have H cc r KC max . Based on each candidate support set, the nonzero entries of the estimated sparse vector are calculated using least squares algorithm. The computational cost of this refined estimation process is of the order of C KHcc OðLSÞ, where O(LS) indicates the order of the computational complexity for the least squares algorithm and OðLSÞ ¼ OðKMÞ according to [20]. 2. a^ contains the indices of the incoherent columns only. In this condition, the estimated support set a^ equals the true support set of the original K-sparse vector due to
the incoherence between the contributing incoherent columns according to Propositions 1–3. We can reconstruct the original sparse vector using the least squares ^ and the computational cost is of algorithm based on a, the order of O(KM). 3. a^ contains the indices of both the condensed columns and the incoherent columns. ^ gives The rough estimate of the true support set, a, correct positions of the nonzero entries corresponding to the contributing condensed columns and the contributing incoherent columns according to Propositions 1–3. We assume that a^ contains the indices of V ðV o KÞ contributing similar column groups and the indices of KV contributing incoherent columns. The combinations from the columns of the contributing similar column groups, together with the contributing incoherent columns form a number of candidate support sets, among which we can find the true one that matches the residual best. So the computational cost is of the order of C VHci OðKMÞ, where Hci denotes the total number of columns from the V contributing similar column groups. We have H ci rVC max .
The computational complexity of the proposed algorithm is summarized in Table 1. The computational costs of the SP and BP algorithms are of the order of O(MN) and LPðM; NÞ, respectively, according to [20]. The symbol LPðM; NÞ indicates the cost of solving a linear program with N variables and M constraints, which is OðM2 N1:5 Þ for an interior-point algorithm. The computational complexity of the SP and BP algorithms is shown in Table 2. The proposed SSMP algorithm is compared with the SP and BP algorithms based on the computational complexity of the online processing considering that the off-line processing is carried out before the measurements arrive and performed only once during the whole process with fixed original sensing matrix. From Tables 1 and 2, it can be seen that under the second condition, the computational complexity of the proposed algorithm is the lowest and comparable to that of the SP algorithm, considering that N′o N and K is far less than M and N. It can also be Table 1 Computational complexity of the proposed algorithm. Condition
Off-line processing
Condition 1
O
Condition 2 Condition 3
O O
ðN1ÞN M 2 ðN1ÞN M 2 ðN1ÞN M 2
Rough estimation
Refined estimation
OðMN′Þ
C KHcc OðKMÞ
OðMN′Þ
O(KM)
OðMN′Þ
C VHci OðKMÞ
Table 2 Computational complexity of the SP and BP algorithms.
Computational complexity
SP
BP
O(MN)
OðM 2 N 1:5 Þ
J. Liu et al. / Signal Processing 95 (2014) 101–110
seen that under the first condition, the proposed algorithm has the highest computational complexity among the three conditions, which relies mainly on Hcc. If the sensing matrix is highly coherent with a large Cmax, Hcc may be very large, resulting in high computational complexity. On the contrary, if the sensing matrix is less coherent with a small Cmax, Hcc will be small and the computational complexity of the proposed algorithm will be comparable to that of the SP algorithm and lower than that of the BP algorithm. In Section 3, three specific examples are presented to show the comparison results of the three algorithms based on computational complexity. 3. Simulation results and analysis The similar sensing matrix pursuit algorithm is proposed to cope with deterministic sensing matrix with high coherence. The deterministic sensing matrix is always related to specific applications. In this simulation, we consider the direction-of-arrival (DOA) estimation problem using compressed sensing in sensor array processing. A sensor array system with half-wavelength spaced uniform linear arrays (ULA) is considered in this example. The ULA system consists of Mantenna antennas, which are used to emit the transmitted signal. Targets may appear at directions represented by DOA angles. The task of signal processing is to estimate the directions to the targets and the corresponding complex amplitudes (DOA estimation, see [31]). We assume that the other parameters such as range and Doppler frequency have been isolated before by appropriate processing. Since the radar scene is generally sparse in practice, compressed sensing is a valid candidate for estimating the DOA angles for multiple targets. To do so, the DOA angle plane is divided into Ngrid fine grids, each cell generally with the same size Δθ. The ith grid represents a DOA angle of θ0 þ ði1ÞΔθ, where θ0 is the initial angle of the DOA plane. As the system has no knowledge of the number and locations of the targets, the information of all the grids in the DOA plane should be considered. Therefore, the steering matrix Φsteering is defined as Φsteering ¼ ½uðθ0 Þu ðθ0 þ ΔθÞ⋯uðθ0 þ ði1ÞΔθÞ⋯uðθ0 þðNgrid 1ÞΔθÞ, where uðÞ is a M antenna 1 steering vector, e.g. h iT ð22Þ uðθ0 Þ ¼ 1 e j2πd sin θ0 =λ ⋯e jðMantenna 1Þ2πd sin θ0 =λ : In (22), d is the distance between the elements of the arrays and λ denotes wavelength. Since a small number of grids are occupied by the targets, the reflection vector x is an N grid 1 sparse vector with the ith element defined as xðiÞ ¼ rðθ0 þ ði1ÞΔθÞ if the ith grid is occupied by a target; otherwise, xðiÞ ¼ 0. rðθ0 þ ði1ÞΔθÞ is the reflection coefficient of the hypothetical target at a DOA angle θ0 þ ði1ÞΔθ. As a result, the Mantenna 1 received complex vector of array observations y could be written as y ¼ Φsteering x þ e;
ð23Þ
where e is a M antenna 1 complex Gaussian noise vector. Though the radar vectors and matrices in (23) are complex valued in contrast to the original compressed sensing environment, it is easy to transfer them to real variables
107
according to [32,8], where the radar case is mapped into the standard compressed sensing framework by breaking the problem into its real and imaginary parts. In the simulation setup, the number of transmit/receive antennas (Mantenna) is set to 10. The antennas transmit independent orthogonal quadrature phase shift keyed (QPSK) waveforms and the carrier frequency is 8.62 GHz. The range of the DOA plane is [01,1801], which is divided into 20 (Ngrid) cells with the initial angle ðθ0 Þ) and angle interval ðΔθÞ equaling 01 and 91, respectively. A signal sparsity level K is chosen such that K r Mantenna =2. The deterministic steering matrix Φsteering is chosen as the original sensing matrix. A support set S of size K is selected uniformly at random, and the original sparse signal vector x is chosen as either Gaussian signal or zero-one signal [19]. In the proposed algorithm, the Gram matrix G of the original sensing matrix Φ (comprised of 20 columns) is built via (5), and the histogram of the similarity values of the original sensing matrix is shown in Fig. 2. Next, we will set the threshold T1 to distinguish the highly coherent columns from the incoherent columns. In this simulation example, the threshold T1 is set to 0.4, resulting in 19 highly coherent columns and 1 incoherent column. The 19 coherent columns are further divided into 7 similar column groups according to the threshold T2 (0.9). Each similar column group is compacted to a condensed column. The obtained condensed columns together with the incoherent columns form a similar compact sensing matrix (comprised of 8 columns). The histogram of the similarity values of the similar compact sensing matrix is shown in Fig. 3. As can be seen, the size of the similar compact sensing matrix is greatly reduced, which reduces the computational complexity efficiently. Furthermore, there is a shift toward the origin of the histogram from Fig. 2 (the original sensing matrix) to Fig. 3 (the similar compact sensing matrix). The tail representing the higher values in Fig. 2 disappears in Fig. 3. Therefore the coherence of the similar compact sensing matrix is far less than that of the original sensing matrix. The simulations are carried out to compare the accuracy of different reconstruction algorithms empirically. The proposed algorithm is compared with the SP and BP
350 300 250 200 150 100 50 0
0
0.2
0.4
0.6
0.8
1
Fig. 2. Histogram of the similarity values (absolute off-diagonal entries of G) of the original sensing matrix.
108
J. Liu et al. / Signal Processing 95 (2014) 101–110
1
40
0.9
Frequency of Exact Reconstruction
45
35 30 25 20 15 10 5
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1
0
0
0.2
0.4
0.6
0.8
Fig. 3. Histogram of the similarity values (absolute off-diagonal entries of G) of the similar compact sensing matrix.
2
2.5
3
3.5
4
4.5
5
Sparsity Level: K Fig. 5. Frequency of exact reconstruction for the Gaussian signal.
Table 3 Parameter setting of the three examples.
1
Frequency of Exact Reconstruction
1.5
1
0.9 0.8
Example number
N
M
K N’
D
Cmax Hcc
Example 1
1000 500 5
Example 2
1000 500 5 100 100
10
Example 3
1000 500 5 500 500
2
C KHcc
0.7 0.6 0.5
10
10 100
500 2:5524 1011 50 2:1188 106 10 252
0.4
Table 4 Computational complexity of the three algorithms for the three specific examples.
0.3 0.2 0.1
Example number
Proposed algorithm
SP
BP
Example 1
Oð1014 Þ
Oð105 Þ
0 1
1.5
2
2.5
3
3.5
4
4.5
5
Sparsity Level: K Fig. 4. Frequency of exact reconstruction for the zero-one signal.
algorithms. In the experiments, the BP algorithm uses the default settings [33] (BP tools are from SparseLab [34]), and the SP algorithm uses the parameters given in [19]. Five hundred Monte Carlo simulations are performed for each fixed value of K (size of the support set). The reconstruction is considered to be exact when the l2 norm of the difference between the original signal x and the ^ 2 reconstructed one x^ is smaller than 10 5, that is ‖xx‖ o 105 . Figs. 4 and 5 present the reconstruction results for binary zero-one and Gaussian sparse signals, respectively. This shows that the reconstruction performance of the proposed algorithm is much better than that of the SP and BP algorithms, with a successful reconstruction probability greater than 0.9. Finally, the proposed SSMP algorithm is compared with the SP and BP algorithms for computational complexity. We consider the worst case where the estimated support set a^ contains the indices of the condensed columns only. Three examples with different Cmax are provided. The first example considers a highly coherent sensing matrix with a large Cmax (100). The original sensing matrix is comprised of 1000 highly coherent columns, which are divided into 10 similar column groups with equal length of 100. While in the second and third examples, less coherent sensing matrices are adopted with Cmax as 10 and 2, respectively.
5
Oð108 Þ
Example 2
Oð10 Þ
Oð10 Þ
Oð108 Þ
Example 3
Oð105 Þ
Oð105 Þ
Oð108 Þ
8
The parameters of the three examples are presented in Table 3. The last column of Table 3 lists the values of C KHcc , which is proportional to the computational complexity of the proposed algorithm. It can be seen that the more coherent the sensing matrix is, the larger the C KHcc will be. The computational complexity of the three algorithms for the three specific examples is listed in Table 4. When the sensing matrix is highly coherent with a large Cmax (100, in the first example), the computational complexity of the proposed algorithm is higher than that of the SP and BP algorithms. On the contrary, if the sensing matrix is less coherent with a small Cmax (2, in the third example), the computational complexity of the proposed algorithm is comparable to that of the SP algorithm and lower than BP algorithm. 4. Conclusions In this paper, a novel algorithm called the similar sensing matrix pursuit (SSMP) is proposed to reconstruct a K-sparse vector based on a deterministic sensing matrix with high coherence. The proposed algorithm consists of two parts: the off-line processing and online processing. The goal of the off-line processing is to find a similar
J. Liu et al. / Signal Processing 95 (2014) 101–110
compact sensing matrix, which contains as much information as possible from the original sensing matrix. The online processing consists of a rough estimation process and a refined estimation process. In the rough estimation process, a subspace pursuit algorithm is used to find a rough estimate of the true support set, which contains the indices of the columns that contribute to the original sparse vector. Three kinds of structures of the estimated support set are considered, and three individual refined estimation processes are carried out under these three conditions. From the simulation results, it can be seen that the proposed algorithm obtains much better performance when coping with the deterministic sensing matrix with high coherence, compared with the SP and BP algorithms.
Acknowledgments This work was supported in part by National 973 project of China through Grant 2013CB329405, the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (61221063), and the Natural Science Table A1 Notations. Variable
Definition
K x N y M Φ δK μðΦÞ φi H λðφi ; φj Þ G T1
Sparsity level K-sparse vector Length of K-sparse vector x Measurement vector Length of measurement vector y Sensing matrix Restricted isometry constant (RIC) Coherence of Φ ith column of Φ Additional sensing matrix Similarity between φi ; φj Gram matrix Threshold set to distinguish the highly coherent columns from the incoherent columns Threshold set to distinguish the highly coherent columns inside a similar column group Similarity distance between φi ; φj
T2 dsimilar ðφi ; φj Þ Ks Γi Shc Sic ιi ζi Nhc Nic D γ ij N Γi Ψ γ iC
A constant in calculating dsimilar The ith similar column group Set of highly coherent columns Set of incoherent columns The ith highly coherent column The ith incoherent column Number of highly coherent columns Number of incoherent columns Number of similar column groups The jth column in the ith similar column group Γi Number of columns in Γi Similar compact sensing matrix
The condensed column from Γi in the similar compact sensing matrix Ψ γ The noise magnitude mr The measurement residual in the OMP algorithm ε An appropriately chosen bound on the noise magnitude The estimated signal x^ O 1 2 K a^ ¼ fa^ ; a^ ; …; a^ g Estimated support set The first element of a^ corresponds to γ iC , the ^a 1 ðγ iC Þ index of the ith condensed column.
109
Table A2 Notations. b^ ¼ fΓ i ; Γ j ⋯; Γ l g A set containing the indices of K similar column groups corresponding to a^ The indices of all the columns in the selected f^ similar column groups in b^ Hcc Hci ϒp
Total number of columns in the K contributing similar column groups Total number of columns from the V contributing similar column groups The pth candidate support set
Υ 1p ðγ i1 Þ
The first element of ϒp corresponding to γ i1
Nco p x^
Total number of combinations The estimate of the original sparse vector based on ϒp
p x^ ϒp p x^ f1;2;…;NgΥ p
p The entries of x^ indexed by iA ϒp p The entries of x^ indexed by iA f1; 2; …; Ngϒp
Φpϒ † x^ min
The columns of Φ indexed by iA ϒp Pseudo-inverse operation The estimate with the least residual The first element of a^ corresponds to ζ i , the ith incoherent column. Number of selected condensed columns The pth combination
1 a^ ðζ i Þ
V rp ¼ fγ i1 ; …; γ j1 g Ncb O(LS) Cmax Mantenna Ngrid Δθ θ0 Φsteering u d λ r S
Number of combinations Order of computational cost for least squares algorithm Size of the largest similar column group Number of antennas Number of fine grids in the DOA angle plane Size of each cell Initial angle of the DOA plane Steering matrix M antenna 1 steering vector Distance between the elements of the arrays Wavelength Reflection coefficient Support set
Foundations of China (nos. 61104051, 61074176 and 61004087), the Fundamental Research Funds for the Central Universities, and the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry. Appendix A. Notations The definitions of variables used in the paper are listed in Tables A1 and A2. References [1] J. Romberg, Imaging via compressive sampling [introduction to compressive sampling and recovery via convex programming], IEEE Signal Processing Magazine 25 (2) (2008) 14–20. [2] Xue Bi, Xiang-dong Chen, Yu Zhang, Bin Liu, Image compressed sensing based on wavelet transform in contourlet domain, Signal Processing 91 (5) (2011) 1085–1092. [3] C. Berger, S. Zhou, J. Preisig, P. Willett, Sparse channel estimation for multicarrier underwater acoustic communication: from subspace methods to compressed sensing, IEEE Transactions on Signal Processing 58 (3) (2010) 1708–1721. [4] W.U. Bajwa, J.D. Haupt, A.M. Sayeed, R.D. Nowak, Compressed channel sensing: a new approach to estimating sparse multipath channels, Proceedings of the IEEE 98 (6) (2010) 1058–1076. [5] M. Herman, T. Strohmer, Compressed sensing radar, IEEE Transactions on Signal Processing 13 (10) (2006) 589–592. [6] Guanghui Zhao, Zhengyang Wang, Qi Wang, Guangming Shi, Fangfang Shen, Robust ISAR imaging based on compressive sensing from noisy measurements, Signal Processing 92 (1) (2012) 120–129.
110
J. Liu et al. / Signal Processing 95 (2014) 101–110
[7] C. Berger, B. Demissie, J. Heckenbach, P. Willett, S. Zhou, Signal processing for passive radar using OFDM waveforms, IEEE Journal of Selected Topics in Signal Processing 4 (1) (2010) 226–238. [8] Joachim H.G. Ender, On compressive sensing applied to radar, Signal Processing 90 (5) (2010) 1402–1414. [9] Y. Yu, A. Petropulu, H. Poor, Measurement matrix design for compressive sensingcbased MIMO radar, IEEE Transactions on Signal Processing 59 (11) (2011) 5338–5352. [10] E. Candès, J. Romberg, T. Tao, Stable signal recovery from incomplete and inaccurate measurements, Communications on Pure and Applied Mathematics 59 (8) (2006) 1207–1223. [11] D. Donoho, Compressed sensing, IEEE Transactions on Information Theory 52 (4) (2006) 1289–1306. [12] E. Candès, Compressive sampling, Proceedings of the International Congress of Mathematicians 3 (2006) 1433–1452. [13] E. Candès, J. Romberg, T. Tao, Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information, IEEE Transactions on Information Theory 52 (2) (2006) 489–509. [14] E. Candès, T. Tao, Decoding by linear programming, IEEE Transactions on Information Theory 51 (12) (2005) 4203–4215. [15] S. Mallat, Z. Zhang, Matching pursuit with time-frequency dictionaries, IEEE Transactions on Signal Processing 41 (12) (1993) 3397–3415. [16] J.A. Tropp, Greed is good: algorithmic results for sparse approximation, IEEE Transactions on Information Theory 50 (10) (2004) 2231–2242. [17] D. Needell, R. Vershynin, Uniform uncertainty principle and signal recovery via regularized orthogonal matching pursuit, Foundations of Computational Mathematics 9 (2009) 317–334. [18] D.L. Donoho, Y. Tsaig, I. Drori, J.-L. Starck, Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit, IEEE Transactions on Information Theory 58 (2) (2012) 1094–1121. [19] W. Dai, O. Milenkovic, Subspace pursuit for compressive sensing signal reconstruction, IEEE Transactions on Information Theory 55 (5) (2009) 2230–2249. [20] D. Needell, J.A. Tropp, COSAMP: iterative signal recovery from incomplete and inaccurate samples, Applied and Computational Harmonic Analysis 26 (3) (2009) 301–321.
[21] H. Huang, A. Makur, Backtracking-based matching pursuit method for sparse signal reconstruction, IEEE Signal Processing Letters 18 (7) (2011) 391–394. [22] J. Bourgain, S. Dilworth, K. Ford, S. Konyagin, D. Kutzarova, Explicit constructions of RIP matrices and related problems, Duke Mathematics Journal 159 (1) (2011) 145–185. [23] A. Cohen, W. Dahmen, R. DeVore, Compressed sensing and best k-term approximation, Journal of American Mathematical Society 22 (2009) 211–231. [24] S.D. Howard, A.R. Calderbank, S.J. Searle, A fast reconstruction algorithm for deterministic compressive sensing using second order reed-muller codes, in: IEEE Conference on Information Sciences and Systems (CISS2008), 2008. [25] A. Amini, F. Marvasti, Deterministic construction of binary, bipolar, and ternary compressed sensing matrices, IEEE Transactions on Information Theory 57 (4) (2011) 2360–2370. [26] E. Candès, J. Romberg, Sparsity and incoherence in compressive sampling, Inverse Problems 23 (3) (2007) 969–985. [27] J. Zhang, D. Zhu, G. Zhang, Adaptive compressed sensing radar oriented toward cognitive detection in dynamic sparse target scene, IEEE Transactions on Signal Processing 60 (4) (2012) 1718–1729. [28] M.F. Duarte, Y.C. Eldar, Structured compressed sensing: from theory to applications, IEEE Transactions on Signal Processing 59 (9) (2011) 4053–4085. [29] D.L. Donoho, M. Elad, V.N. Temlyakov, Stable recovery of sparse overcomplete representations in the presence of noise, IEEE Transactions on Information Theory 52 (1) (2006) 6–18. [30] M. Elad, Optimized projections for compressed sensing, IEEE Transactions on Signal Processing 55 (12) (2007) 5695–5702. [31] A. Gurbuz, J. McClellan, V. Cevher, A compressive beamforming method, in: Proceedings of IEEE International Conference on Acoustics and Speech Signal Processing, 2008, pp. 2617–2620. [32] I. Stojanovic, W. Karl, M. Cetin, Compressed sensing of mono-static and multi-static SAR, in: SPIE Defense and Security Symposium, Algorithms for Synthetic Aperture Radar Imagery XVI. [33] S. Chen, D.L. Donoho, M. Saunders, Atomic decomposition by basis pursuit, SIAM Journal of Science Computing 20 (1998) 33–61. [34] SparseLab [Online]. Available: 〈http://dsp.rice.edu/cs〉.
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Similar sensing matrix pursuit: An efficient reconstruction algorithm to cope with deterministic sensing matrix Jing Liu a,n, Mahendra Mallick b, ChongZhao Han a, XiangHua Yao a, Feng Lian a a b
School of Electronics and Information Engineering, Xi'an Jiaotong University, Xi'an, PR China Anacortes, WA, USA
a r t i c l e i n f o
abstract
Article history: Received 24 September 2012 Received in revised form 1 August 2013 Accepted 19 August 2013 Available online 8 September 2013
Deterministic sensing matrices are useful, because in practice, the sampler has to be a deterministic matrix. It is quite challenging to design a deterministic sensing matrix with low coherence. In this paper, we consider a more general condition, when the deterministic sensing matrix has high coherence and does not satisfy the restricted isometry property (RIP). A novel algorithm, called the similar sensing matrix pursuit (SSMP), is proposed to reconstruct a K-sparse signal, based on the original deterministic sensing matrix. The proposed algorithm consists of off-line and online processing. The goal of the off-line processing is to construct a similar compact sensing matrix containing as much information as possible from the original sensing matrix. The similar compact sensing matrix has low coherence, which guarantees a perfect reconstruction of the sparse vector with high probability. The online processing begins when measurements arrive, and consists of rough and refined estimation processes. Results from our simulation show that the proposed algorithm obtains much better performance while coping with a deterministic sensing matrix with high coherence compared with the subspace pursuit (SP) and basis pursuit (BP) algorithms. & 2013 Elsevier B.V. All rights reserved.
Keywords: Compressed sensing Deterministic sensing matrix Reconstruction algorithm Similar sensing matrix pursuit algorithm
1. Introduction Compressed sensing has received considerable attention recently, and has been applied successfully in diverse fields, e.g. image processing [1,2], underwater acoustic communication [3], wireless communication [4] and radar systems [5–9]. The central goal of compressed sensing is to capture attributes of a signal using very few measurements. In most work to date, this broader objective is exemplified by the important special case in which a K-sparse vector x A RN (with N large) is to be reconstructed from a small number M of linear measurements with K oM oN. K-sparse signals are the signals that can be represented by K significant coefficients over an N-dimensional basis. This can be compactly described via y ¼ Φx þe: n
ð1Þ
Corresponding author. E-mail address: [email protected] (J. Liu).
0165-1684/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sigpro.2013.08.009
Here, y A RM denotes a measurement vector, Φ represents an M N sensing matrix, and e is an M 1 noise vector. The two fundamental questions in compressed sensing are as follows: how to construct suitable sensing matrix Φ, and how to recover the K-sparse vector x from the measurement vector y efficiently. Tables A1 and A2 in Appendix A list notations for variables used in the paper. In early work of compressed sensing, the entries of the sensing matrix Φ are generated by an independent and identically distributed (i.i.d) Gaussian or Bernoulli process, or from random Fourier ensembles [10–12]. In general, the exact solution to the above second question is shown to be an NP-hard problem [13,14]. If the number of samples (M) exceeds the lower bound of M 4OðK log ðN=KÞÞ, l1 minimization (e.g. BP algorithm) can be performed instead of the exact l0 minimization with the same solution for almost all the possible inputs [14]. An alternative approach to sparse signal recovery is based on the idea of iterative greedy pursuit, and tries to approximate the solution to
102
J. Liu et al. / Signal Processing 95 (2014) 101–110
l0 minimization directly. The greedy algorithms include matching pursuit (MP) [15], orthogonal matching pursuit (OMP) [16], regularized OMP (ROMP) [17], stagewise OMP (StOMP) [18], SP [19], compressive sampling matching pursuit (CoSaMP) [20] and backtracking-based matching pursuit (BAOMP) [21], etc. The reconstruction complexity of these approximate algorithms is significantly lower than that of BP algorithm. In compressed sensing, one of the well-studied conditions on the sensing matrix, which guarantees stable recovery for a number of reconstruction algorithms, is the RIP [13,14]. If a sensing matrix Φ whose column vectors have unit norm and satisfies ð1δK Þ‖x‖22 r‖Φx‖22 rð1 þ δK Þ‖x‖22
ð2Þ
for all possible K-sparse vectors with restricted isometry constant (RIC) δK , then Φ is said to obey K-RIP with δK . The RIP with suitable constant δK guarantees perfect reconstruction [13,14], but it is very hard to check whether a sensing matrix satisfies RIP or not. Coherence, the maximal correlation between two columns in a sensing matrix, is also a well-known performance measure for sensing matrices. For a matrix Φ with columns φ1 ; φ2 ; …; φN , the coherence of Φ is defined as μðΦÞ ¼
max
1 r i;j r N and
jφTi φj j : i a j J φ i J J φj J
ð3Þ
Coherence plays a central role in the sensing matrix construction, because small coherence implies the RIP [22]. In this paper, we are interested in deterministic sensing matrices. Deterministic sensing matrices are useful because in practice, the sampler has to be a deterministic matrix. Although random matrices perform quite well on the average, there is no guarantee that a specific realization works. For the deterministic approaches, the Vandermond matrices seem to be good options, since any K columns of a K N Vandermond matrix are linearly independent. However, when N increases, the constant δK rapidly approaches 1 and some of the K K submatrices become ill-conditioned [23]. A connection between the coding theory and sensing matrices is established in [24] where second order Reed– Muller codes are used to construct bipolar matrices. However, they lack a guarantee on the RIP order. In [25], the authors propose a series of deterministic sensing matrices, the binary, bipolar, and ternary compressed sensing matrices which satisfy the RIP condition. The key concept of coherence is extended to pairs of orthonormal bases. This enables a new choice of the sensing matrices: one simply selects an orthonormal basis that is incoherent with the sparsity basis, and obtains measurements by selecting a subset of the coefficients of the signal in the chosen basis [26]. This approach has successful applications in radar systems [9,27], where an additional sensing matrix H is introduced and the received signal is compressed further by making nonadaptive, linear projections of the direct data sampled at the Nyquist frequency. However, neither of these algorithms mention the hardware implementation of the additional sensing matrix, which is very complex and expensive. In practice, it is challenging to design a deterministic sensing matrix having low coherence. In this paper, we
consider a more general condition when the deterministic sensing matrix has high coherence and does not satisfy the RIP condition. A novel algorithm, called the SSMP algorithm, is proposed to reconstruct the K-sparse signal based on the original deterministic sensing matrix. The proposed algorithm consists of two parts: the off-line processing and the online processing. The goal of the off-line processing is to construct a similar compact sensing matrix with low coherence, which contains as much information as possible from the original sensing matrix. The online processing begins when the measurements arrive, which consists of a rough estimation process and a refined estimation process. In the rough estimation process, an SP algorithm is used to find a rough estimate of the true support set, which contains the indices of the columns that contribute to the original sparse vector. Three kinds of structures of the estimated support set are considered, and three individual refined estimation processes are carried out under these three conditions. We observe from simulation results that the proposed algorithm obtains much better performance when coping with the deterministic sensing matrix with high coherence compared with the SP and BP algorithms. The paper is organized as follows. Section 2 introduces the proposed similar sensing matrix pursuit algorithm. Section 3 presents simulation results, and Section 4 summarizes conclusions. 2. Similar sensing matrix pursuit (SSMP) algorithm Recently, algorithms used to cope with deterministic sensing matrices focus on designing a sensing matrix which satisfies the RIP condition (or with low coherence). However, in practice it is challenging to design a deterministic sensing matrix having a very small restricted isometry constant (coherence). In this paper, we consider a more general condition when a deterministic sensing matrix has high coherence and does not satisfy the RIP condition. We concentrate in developing a novel reconstruction algorithm rather than building a deterministic sensing matrix with low coherence. A novel algorithm called the SSMP is proposed to reconstruct the K-sparse signal, based on the original deterministic sensing matrix. This section introduces the proposed SSMP algorithm. First, the key component of the proposed algorithm, the similar compact sensing matrix, is introduced in Section 2.1, and then the complete algorithm is described in Section 2.2. The complexity analysis of the proposed algorithm is presented in Section 2.3. 2.1. Construction of the similar compact sensing matrix The construction process of the similar compact sensing matrix is based on the similarity analysis of the original sensing matrix. In this paper, similarity is defined as the absolute and normalized inner product between any two different columns of the original sensing matrix Φ: λðφi ; φj Þ ¼
jφTi φj j ; J φi J J φj J
1 ri; j r N
and
i aj:
ð4Þ
J. Liu et al. / Signal Processing 95 (2014) 101–110
It can be seen that coherence is the largest similarity among the columns of a matrix. Any two columns with high similarity are coherent with each other, and vice versa. Therefore, the similarity can be used to distinguish the coherent columns from the incoherent columns of the original sensing matrix. A different way to understand similarity and coherence is by considering the Gram matrix G which is defined as ~ T Φ; ~ G¼Φ
ð5Þ
~ is the normalized sensing matrix obtained from where Φ the original sensing matrix with each column normalized. The off-diagonal entries in G are the similarity values defined in (4). The coherence is the off-diagonal entry with the largest magnitude. A threshold T1 is set properly to distinguish the highly coherent columns from the incoherent columns based on similarity. For each calculated similarity value, fλðφi ; φj Þ; i ¼ 1; …; N; j ¼ 1; …; N; ia jg, if λðφi ; φj Þ is greater than T1, the columns i and j are added to the set of highly coherent columns. The remaining columns that do not belong to the set of highly coherent columns form the set of incoherent columns. Fig. 1 shows the classification results in a similarity plane, where the incoherent columns and the highly coherent columns are indicated by black circles and black squares, respectively. In the similarity plane, the distance between any two columns ðφi ; φj Þ is defined as the similarity distance dsimilar, which is inversely proportional to the similarity value λðφi ; φj Þ, as follows: dsimilar ðφi ; φj Þ ¼ K s =λðφi ; φj Þ;
ð6Þ
where Ks is a constant. Therefore the closer the two columns located in the similarity plane, the more similar they are (with larger similarity). Fig. 1 also shows that any two highly coherent columns may not be close to each other (coherent with each other) in the similarity plane. The highly coherent columns are then further divided into a set of similar column groups (indicated by large circles in Fig. 1) to ensure that the highly coherent columns in a similar column group are similar with each other. Moreover, the similar column groups have the following property. Property 1. For any two similar column groups, which are defined as Γ1 ¼ fα1 ; …; αN g and Γ2 ¼ fβ1 ; …; βM g, the similarity between any two columns from these two different similar Highly coherent column
dsimilar Similar column group
dsimilar
dsimilar Incoherent column
Fig. 1. Classification results in a similarity plane.
103
column groups, e.g. αi A Γ1 and βj A Γ2 , is no larger than T1, as λðαi ; βj Þ rT 1 ;
αi A Γ1
and
βj A Γ2 ;
ð7Þ
while the similarity between any two columns in the same similar column group (e.g. Γ1 ) is larger than a threshold T2 as λðαi ; αj Þ 4T 2 ;
α i ; α j A Γ1
and
ð8Þ
ia j:
Property 1 stipulates that any column in one specific similar column group is highly coherent with other columns inside the same group, while incoherent with any column outside the group (including the highly coherent columns in other similar column groups, and incoherent columns). This is shown in Fig. 1 that any column in one large circle (similar column group) is very close to other columns inside the circle, while far from any column outside the circle. Therefore, any column in one specific similar column group represents the characteristics of other columns in the same group. In this paper, we consider reducing the original sensing matrix to a similar compact sensing matrix by drawing a column from each similar column group, while keeping the incoherent columns unchanged in the new sensing matrix. The newly built similar compact sensing matrix contains as much information as possible from the original sensing matrix. Next, the construction process of the similar compact sensing matrix is briefly introduced. Construction of the similar compact sensing matrix: The set of highly coherent columns and the set of incoherent columns are denoted by Shc ¼ fι1 ; ι2 ; ⋯; ιNhc g and Sic ¼ fζ 1 ; ζ 2 ; ⋯; ζ Nic g, respectively. Nhc is the number of highly coherent columns, and Nic is the number of incoherent columns with Nhc þN ic ¼ N. The set of highly coherent columns Shc is further divided into D similar column groups, fΓ1 ; Γ2 ; …; Γi ; …; ΓD g. We assume that M r D oN. Each similar column group contains more than one highly coherent columns, e.g. Γi ¼ fγi1 ; …; γiN i g, where N Γi indiΓ
cates the number of columns in Γi . Each similar column group is compacted to a condensed column. Here we just select a column from each similar column group randomly considering that the columns in one group are highly coherent and very similar to each other. Finally, we obtain a similar compact sensing matrix by combining the condensed columns and the incoherent columns, as Ψ ¼ 1 2 N ic ½γ1C ; γ 2C ; …; γ iC ; …; γD ,where γiC denotes a conC ; ζ ; ζ ; …; ζ
densed column from Γi . The similar compact sensing matrix Ψ is of size M ðD þ Nic Þ. We have the following propositions for the similar compact sensing matrix. Proposition 1. The coherence of the similar compact sensing matrix is less than or equal to T1. Proof. According to Property 1, the similarity between any two condensed columns is no larger than T1. Moreover, from the division process of the coherent columns and the incoherent columns, the similarity between any two incoherent columns, or between a condensed column and a incoherent column, is no larger than T1. As a result, the similarity between any two columns of the similar compact sensing matrix is no larger than T1. Therefore, the coherence
104
J. Liu et al. / Signal Processing 95 (2014) 101–110
of the similar compact sensing matrix, which is the largest similarity, is less than or equal to T1. □ Proposition 2 (Duarte and Eldar [28], Donoho et al. [29]). Let the signal x be a K-sparse vector and write y ¼ Φ x þe. Denote γ as the noise magnitude ðγ ¼ ‖e‖2 Þ, and mr as the measurement residual in the OMP algorithm. Suppose that K r ð1=μðΦÞ þ1Þ=4 and ε Z γ, where ε is an appropriately chosen bound on the noise magnitude. We have x^ ¼ arg min‖x‖1 x A RN
subject to ‖yΦx‖2 rε:
ð9Þ
Then the estimate x^ in (9) has error bounded by γ þε ^ 2 r pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; ‖xx‖ 1μðΦÞð4K1Þ
ð10Þ
where μðΦÞ denotes the coherence of the sensing matrix Φ. The estimate x^ of the OMP algorithm with halting criterion ‖mr‖2 rγ has the error bound γ ^ 2 r pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; ‖xx‖ 1μðΦÞð4K1Þ
ð11Þ
provided that γ rAð1μðΦÞð2K1ÞÞ=2 for OMP, with A being a positive lower bound on the magnitude of the nonzero entries of x. Proposition 3 (Setting of the threshold T1). In order to guarantee a perfect reconstruction of the sparse vector with large probability, the threshold T1 should satisfy (12) according to Propositions 1 and 2: T 1 r1=ð4K1Þ:
ð12Þ
In this paper the threshold T1 is set as T 1 ¼ 1=ð4K1Þ. However, a small T1 will result in a large number of highly coherent columns. While the restricted selection of T1 is true from a worst-case standpoint, it turns out that the coherence as defined previously does not do justice to the actual behavior of sparse representations and pursuit algorithms' performance. Thus, if we relax our expectations and allow a small fraction of failed reconstructions, then values which are substantially above the bound still lead to successful compressed sensing [30]. In this work, the threshold T1 is set to 0.4 in the simulation. The threshold T2 should be set to a large value to ensure that the columns in a similar column group are highly coherent (similar) with each other. The larger the T2 is, the more similar the columns in a similar column group will be, and vice versa. If T2 is set to a small value, the condensed column, which is selected from a similar column group randomly, cannot represent the characteristics of other columns in the same group. Here the threshold T2 is set to 0.9. Therefore, we can see that the proposed SSMP algorithm is not sensitive to the choice of the threshold values. We only need to set T1 to a small value to ensure the low coherence of the similar compact sensing matrix, and set T2 to a large value to ensure that the columns in a similar column group are highly coherent (similar) with each other.
2.2. The similar sensing matrix pursuit algorithm In this paper, a new algorithm called the similar sensing matrix pursuit, is used to cope with the problem induced by the highly coherent columns. Firstly, the original sensing matrix is transformed to a similar compact sensing matrix based on similarity analysis. The similar compact sensing matrix has low coherence, which guarantees a perfect reconstruction of the sparse vector with a large probability according to Propositions 1–3. An OMP-type algorithm (SP algorithm) is then used to find a rough estimate of the true support set, based on the similar compact sensing matrix and the measurement vector, which contains the indices of the columns contributing to the original sparse vector. The contributing columns may include the incoherent columns and the condensed columns. Since a condensed column is selected from a similar column group randomly, we can only obtain the indices of the contributing similar column groups. From the rough estimation process, we cannot decide the real contributing columns in the contributing similar column groups. Next, a refined estimation process is adopted to provide equal opportunities to all the columns in the contributing similar column groups. All the combinations of the columns in the contributing similar column groups and the contributing incoherent columns are listed, with each combination forming a candidate support set. Accordingly, we can obtain a candidate estimate of the original sparse vector using the least squares algorithm based on each candidate support set [16,19]. Finally, we can find the estimated sparse vector matching the residual best. The proposed algorithm consists of off-line processing and online processing. The off-line processing majors in the construction of the similar compact sensing matrix. The online processing consists of a rough estimation process and a refined estimation process. The main sparse reconstruction is carried out in the rough estimation process, where the SP algorithm is used to find a rough estimate of the true support set. The refined estimation process is a combinational searching process among the candidate support sets for the one matching the true support set best. The estimated support set may only include the indices of the contributing condensed columns or the contributing incoherent columns, besides containing both of them. Three individual refined estimation processes are carried out under the above three conditions. Under the first condition, the true nonzero entries of the original sparse vector fall in a number of similar column groups. All the combinations of the columns in the contributing similar column groups are listed with each combination forming a candidate support set. A candidate estimate of the original sparse vector is obtained based on each candidate support set. We can then find the estimated sparse vector matching the residual best. Under the second condition, when the contributing columns include the incoherent columns only, the nonzero entries are easily identified, since the compressed sensing algorithm can clearly identify the entries corresponding to the incoherent columns of the sensing matrix. Finally, the combination of the two algorithms from the first and the second conditions can be used to estimate the original sparse vector under the third condition. The detailed procedure of the proposed algorithm is presented next.
J. Liu et al. / Signal Processing 95 (2014) 101–110
Algorithm 1. Similar sensing matrix pursuit algorithm. Input: Sensing matrix Φ, measurement vector y Output: The estimated signal x^ O 1. Construction of the similar compact sensing matrix: The process is same with that in Section 2.1. 2. Rough estimation: The SP algorithm is used to find a rough estimate of the true support set based on the measurement vector y and the similar compact sensing matrix Ψ. The estimated support set is represented as 1 2 K a^ ¼ fa^ ; a^ ; …; a^ g. The estimated support set a^ contains the indices of the condensed columns or the indices of the incoherent columns, or both, which contribute to the original sparse vector. 3. Refined estimation: Three individual refined estimation processes are carried out under the following three conditions: (a) a^ contains the indices of the condensed columns only: (i) The estimated support set a^ is represented 1 2 K 1 as a^ ¼ fa^ ðγ iC Þ; a^ ðγ jC Þ; …; a^ ðγ lC Þg, where a^ ðγ iC Þ indicates that the first element of a^ corresponds to γ iC , the index of the ith condensed column. Accordingly, we can obtain a set b^ containing the indices of K similar column ^ as b^ ¼ fΓ i ; Γ j ⋯; Γ l g. groups corresponding to a, The indices of all the columns in the selected similar column groups in b^ are listed and form a final set f^ , where f^ ¼ fγ i1 ; …; γ iN i ; …; γ j1 ; …; Γ
γ jN j ; …; γ l1 ; …; γ lN l g. We assume that the total Γ
Γ
number of the indices in f^ is Hcc. combinations based on the column indices in f^ , and each combination forms a candidate support set, e.g. the pth candidate support set is represented as ϒp ¼ fΥ 1p ðγ i1 Þ; Υ 2p ðγ j1 Þ; …; Υ Kp ðγ l1 Þg; p ¼ 1; 2; …; Nco , where Υ 1p ðγ i1 Þ indicates the first element of ϒp corresponding to γ i1 , and Nco indicates the total number of combinations. The proposed algorithm then solves a least squares problem to approximate the nonzero entries of the original sparse vector on each candidate support set ðϒp ; p ¼ 1; 2; …; Nco Þ, and sets other entries to zero, resulting in an estimate of the original p sparse vector, x^ [16,19]. We can obtain the p estimate x^ as p x^ ϒp ¼ Φ†ϒp y;
ð13Þ
p x^ f1;2;…;Ngϒp ¼ 0;
ð14Þ
where † indicates the pseudo-inverse operation. The matrix Φϒp consists of the columns of p Φ with indices i A ϒp , x^ ϒp is composed of the p
entries of x^ indexed by i A ϒp , and x^ f1;2;…;Ngϒp is composed of the entries of x^ iA f1; 2; …; Ngϒp [19].
p
(iii) Final estimate: Among the obtained estimates 1 2 p N fx^ ; x^ ; …; x^ ; f⋯; x^ co g, find the estimate with the least residual, x^ min , and set x^ min as the output estimated signal: x^ O ¼ x^ min : ð15Þ (b) a^ contains the indices of the incoherent columns only: (i) The estimated support set a^ can be repreT
1 2 K sented as a^ ¼ fa^ ðζ i Þ; a^ ðζ j Þ; …; a^ ðζ l Þg , where 1 i a^ ðζ Þ indicates that the first element of a^
corresponds to ζ i , the ith incoherent column. The estimated support set a^ equals the true support set of the original K-sparse vector due to the incoherence between the contributing incoherent columns according to Propositions 1–3. (ii) Estimation: We can obtain the estimate of the original sparse vector using least squares algorithm [16,19] based on the true support set a^ as x^ a^ ¼ Φ†a^ y;
ð16Þ
x^ f1;2;…;Nga^ ¼ 0:
ð17Þ
(iii) Final estimate: x^ O ¼ x^
ð18Þ
(c) a^ contains the indices of both the condensed columns and the incoherent columns: (i) The estimated support set a^ is represented as 1 V V þ1 k K a^ ¼ fa^ ðγ iC Þ; …; a^ ðγ jC Þ; a^ ðζ Þ; …; a^ ðζ l Þg, 1
(ii) Refined estimation: List C KH cc
p
105
indexed by
V
where fa^ ðγ iC Þ; …; a^ ðγ jC Þg correspond selected condensed columns, and fa^
to
V þ1
V
ðζ k Þ;
K
⋯; a^ ðζ l Þg correspond to KV selected incoherent columns. Accordingly, we can obtain a set b^ containing the indices of V similar column 1 V groups corresponding to fa^ ðγ iC Þ; ⋯; a^ ðγ jC Þg, as b^ ¼ fΓ i ; ⋯; Γ j g. The indices of all the columns in
the selected similar column groups in b^ are listed and form a final set f^ , where f^ ¼ fγ i ; …; 1
γ iN i ; …; γ j1 ; …; γ jN j g. We assume that the total Γ
Γ
number of the indices in f^ is Hci. (ii) Refined estimation: List C VH ci
combinations based on the column indices in f^ , and the pth combination can be represented as rp ¼ fγ i1 ; …;
γ j1 g; p ¼ 1; 2; …; N cb , where Ncb is the total number of combinations. Each combination together with all the selected incoherent columns, fζ k ; …; ζ l g, forms a candidate support set. The pth candidate support set is represented as ϒp ¼ fΥ 1p ðγ i1 Þ; …; Υ Vp ðγ j1 Þ; Υ Vp þ 1 ðζ k Þ; …; Υ Kp ðζ l Þg; p ¼ 1; 2; …; N cb . The proposed algorithm then solves a least squares problem to approximate the nonzero entries of the original sparse vector on each candidate support set ðϒp ; p ¼ 1; 2; …; Ncb Þ, and sets other entries to
106
J. Liu et al. / Signal Processing 95 (2014) 101–110
zero, resulting in an estimate of the original p sparse vector, x^ [16,19]. We can obtain the p estimate x^ as p x^ ϒp ¼ Φ†ϒp y;
ð19Þ
p x^ f1;2;…;Ngϒp ¼ 0:
ð20Þ
1 2 (iii) Final estimate: Among the estimates fx^ ; x^ ; …; p N x^ ; …; x^ cb g, find the estimate with the least residual, x^ min , and set x^ min as the output estimated signal: x^ O ¼ x^ min : ð21Þ
2.3. Complexity analysis In this section we will analyze the computational complexity of the proposed similar sensing matrix pursuit algorithm, and further make a comparison with the SP and BP algorithms based on the computational complexity. The proposed approach consists of off-line processing and online processing. The off-line processing transforms the M N original sensing matrix to a M N′ðN′ oNÞ similar compact sensing matrix based on similarity analysis. The computational complexity of the off-line processing mainly focuses on the computation of the similarity between any two columns of the original sensing matrix, which is of the order of OðððN1ÞN=2ÞM). The online processing consists of the rough and refined estimation processes. In the rough estimation process, an SP algorithm is used to find a rough estimate of the true support set. The computational complexity of the rough estimation process is same as that of the SP algorithm and is of the order of OðMN′Þ according to [20]. Three individual refined estimation processes are then carried out under three different conditions mentioned in Section 2.2. The complexity analysis will be based on these three conditions: 1. a^ contains the indices of the condensed columns only. ^ gives The rough estimate of the true support set, a, correct positions of the nonzero entries corresponding to the contributing condensed columns according to Propositions 1–3 They are the indices of K similar column groups contributing to the original sparse vector. In the refined estimation process, C KHcc combinations are listed out as the candidate support sets, where Hcc indicates the total number of columns in the K contributing similar column groups. Assuming Cmax as the size of the largest similar column group, i.e. the number of columns contained in the largest similar column group, we have H cc r KC max . Based on each candidate support set, the nonzero entries of the estimated sparse vector are calculated using least squares algorithm. The computational cost of this refined estimation process is of the order of C KHcc OðLSÞ, where O(LS) indicates the order of the computational complexity for the least squares algorithm and OðLSÞ ¼ OðKMÞ according to [20]. 2. a^ contains the indices of the incoherent columns only. In this condition, the estimated support set a^ equals the true support set of the original K-sparse vector due to
the incoherence between the contributing incoherent columns according to Propositions 1–3. We can reconstruct the original sparse vector using the least squares ^ and the computational cost is of algorithm based on a, the order of O(KM). 3. a^ contains the indices of both the condensed columns and the incoherent columns. ^ gives The rough estimate of the true support set, a, correct positions of the nonzero entries corresponding to the contributing condensed columns and the contributing incoherent columns according to Propositions 1–3. We assume that a^ contains the indices of V ðV o KÞ contributing similar column groups and the indices of KV contributing incoherent columns. The combinations from the columns of the contributing similar column groups, together with the contributing incoherent columns form a number of candidate support sets, among which we can find the true one that matches the residual best. So the computational cost is of the order of C VHci OðKMÞ, where Hci denotes the total number of columns from the V contributing similar column groups. We have H ci rVC max .
The computational complexity of the proposed algorithm is summarized in Table 1. The computational costs of the SP and BP algorithms are of the order of O(MN) and LPðM; NÞ, respectively, according to [20]. The symbol LPðM; NÞ indicates the cost of solving a linear program with N variables and M constraints, which is OðM2 N1:5 Þ for an interior-point algorithm. The computational complexity of the SP and BP algorithms is shown in Table 2. The proposed SSMP algorithm is compared with the SP and BP algorithms based on the computational complexity of the online processing considering that the off-line processing is carried out before the measurements arrive and performed only once during the whole process with fixed original sensing matrix. From Tables 1 and 2, it can be seen that under the second condition, the computational complexity of the proposed algorithm is the lowest and comparable to that of the SP algorithm, considering that N′o N and K is far less than M and N. It can also be Table 1 Computational complexity of the proposed algorithm. Condition
Off-line processing
Condition 1
O
Condition 2 Condition 3
O O
ðN1ÞN M 2 ðN1ÞN M 2 ðN1ÞN M 2
Rough estimation
Refined estimation
OðMN′Þ
C KHcc OðKMÞ
OðMN′Þ
O(KM)
OðMN′Þ
C VHci OðKMÞ
Table 2 Computational complexity of the SP and BP algorithms.
Computational complexity
SP
BP
O(MN)
OðM 2 N 1:5 Þ
J. Liu et al. / Signal Processing 95 (2014) 101–110
seen that under the first condition, the proposed algorithm has the highest computational complexity among the three conditions, which relies mainly on Hcc. If the sensing matrix is highly coherent with a large Cmax, Hcc may be very large, resulting in high computational complexity. On the contrary, if the sensing matrix is less coherent with a small Cmax, Hcc will be small and the computational complexity of the proposed algorithm will be comparable to that of the SP algorithm and lower than that of the BP algorithm. In Section 3, three specific examples are presented to show the comparison results of the three algorithms based on computational complexity. 3. Simulation results and analysis The similar sensing matrix pursuit algorithm is proposed to cope with deterministic sensing matrix with high coherence. The deterministic sensing matrix is always related to specific applications. In this simulation, we consider the direction-of-arrival (DOA) estimation problem using compressed sensing in sensor array processing. A sensor array system with half-wavelength spaced uniform linear arrays (ULA) is considered in this example. The ULA system consists of Mantenna antennas, which are used to emit the transmitted signal. Targets may appear at directions represented by DOA angles. The task of signal processing is to estimate the directions to the targets and the corresponding complex amplitudes (DOA estimation, see [31]). We assume that the other parameters such as range and Doppler frequency have been isolated before by appropriate processing. Since the radar scene is generally sparse in practice, compressed sensing is a valid candidate for estimating the DOA angles for multiple targets. To do so, the DOA angle plane is divided into Ngrid fine grids, each cell generally with the same size Δθ. The ith grid represents a DOA angle of θ0 þ ði1ÞΔθ, where θ0 is the initial angle of the DOA plane. As the system has no knowledge of the number and locations of the targets, the information of all the grids in the DOA plane should be considered. Therefore, the steering matrix Φsteering is defined as Φsteering ¼ ½uðθ0 Þu ðθ0 þ ΔθÞ⋯uðθ0 þ ði1ÞΔθÞ⋯uðθ0 þðNgrid 1ÞΔθÞ, where uðÞ is a M antenna 1 steering vector, e.g. h iT ð22Þ uðθ0 Þ ¼ 1 e j2πd sin θ0 =λ ⋯e jðMantenna 1Þ2πd sin θ0 =λ : In (22), d is the distance between the elements of the arrays and λ denotes wavelength. Since a small number of grids are occupied by the targets, the reflection vector x is an N grid 1 sparse vector with the ith element defined as xðiÞ ¼ rðθ0 þ ði1ÞΔθÞ if the ith grid is occupied by a target; otherwise, xðiÞ ¼ 0. rðθ0 þ ði1ÞΔθÞ is the reflection coefficient of the hypothetical target at a DOA angle θ0 þ ði1ÞΔθ. As a result, the Mantenna 1 received complex vector of array observations y could be written as y ¼ Φsteering x þ e;
ð23Þ
where e is a M antenna 1 complex Gaussian noise vector. Though the radar vectors and matrices in (23) are complex valued in contrast to the original compressed sensing environment, it is easy to transfer them to real variables
107
according to [32,8], where the radar case is mapped into the standard compressed sensing framework by breaking the problem into its real and imaginary parts. In the simulation setup, the number of transmit/receive antennas (Mantenna) is set to 10. The antennas transmit independent orthogonal quadrature phase shift keyed (QPSK) waveforms and the carrier frequency is 8.62 GHz. The range of the DOA plane is [01,1801], which is divided into 20 (Ngrid) cells with the initial angle ðθ0 Þ) and angle interval ðΔθÞ equaling 01 and 91, respectively. A signal sparsity level K is chosen such that K r Mantenna =2. The deterministic steering matrix Φsteering is chosen as the original sensing matrix. A support set S of size K is selected uniformly at random, and the original sparse signal vector x is chosen as either Gaussian signal or zero-one signal [19]. In the proposed algorithm, the Gram matrix G of the original sensing matrix Φ (comprised of 20 columns) is built via (5), and the histogram of the similarity values of the original sensing matrix is shown in Fig. 2. Next, we will set the threshold T1 to distinguish the highly coherent columns from the incoherent columns. In this simulation example, the threshold T1 is set to 0.4, resulting in 19 highly coherent columns and 1 incoherent column. The 19 coherent columns are further divided into 7 similar column groups according to the threshold T2 (0.9). Each similar column group is compacted to a condensed column. The obtained condensed columns together with the incoherent columns form a similar compact sensing matrix (comprised of 8 columns). The histogram of the similarity values of the similar compact sensing matrix is shown in Fig. 3. As can be seen, the size of the similar compact sensing matrix is greatly reduced, which reduces the computational complexity efficiently. Furthermore, there is a shift toward the origin of the histogram from Fig. 2 (the original sensing matrix) to Fig. 3 (the similar compact sensing matrix). The tail representing the higher values in Fig. 2 disappears in Fig. 3. Therefore the coherence of the similar compact sensing matrix is far less than that of the original sensing matrix. The simulations are carried out to compare the accuracy of different reconstruction algorithms empirically. The proposed algorithm is compared with the SP and BP
350 300 250 200 150 100 50 0
0
0.2
0.4
0.6
0.8
1
Fig. 2. Histogram of the similarity values (absolute off-diagonal entries of G) of the original sensing matrix.
108
J. Liu et al. / Signal Processing 95 (2014) 101–110
1
40
0.9
Frequency of Exact Reconstruction
45
35 30 25 20 15 10 5
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1
0
0
0.2
0.4
0.6
0.8
Fig. 3. Histogram of the similarity values (absolute off-diagonal entries of G) of the similar compact sensing matrix.
2
2.5
3
3.5
4
4.5
5
Sparsity Level: K Fig. 5. Frequency of exact reconstruction for the Gaussian signal.
Table 3 Parameter setting of the three examples.
1
Frequency of Exact Reconstruction
1.5
1
0.9 0.8
Example number
N
M
K N’
D
Cmax Hcc
Example 1
1000 500 5
Example 2
1000 500 5 100 100
10
Example 3
1000 500 5 500 500
2
C KHcc
0.7 0.6 0.5
10
10 100
500 2:5524 1011 50 2:1188 106 10 252
0.4
Table 4 Computational complexity of the three algorithms for the three specific examples.
0.3 0.2 0.1
Example number
Proposed algorithm
SP
BP
Example 1
Oð1014 Þ
Oð105 Þ
0 1
1.5
2
2.5
3
3.5
4
4.5
5
Sparsity Level: K Fig. 4. Frequency of exact reconstruction for the zero-one signal.
algorithms. In the experiments, the BP algorithm uses the default settings [33] (BP tools are from SparseLab [34]), and the SP algorithm uses the parameters given in [19]. Five hundred Monte Carlo simulations are performed for each fixed value of K (size of the support set). The reconstruction is considered to be exact when the l2 norm of the difference between the original signal x and the ^ 2 reconstructed one x^ is smaller than 10 5, that is ‖xx‖ o 105 . Figs. 4 and 5 present the reconstruction results for binary zero-one and Gaussian sparse signals, respectively. This shows that the reconstruction performance of the proposed algorithm is much better than that of the SP and BP algorithms, with a successful reconstruction probability greater than 0.9. Finally, the proposed SSMP algorithm is compared with the SP and BP algorithms for computational complexity. We consider the worst case where the estimated support set a^ contains the indices of the condensed columns only. Three examples with different Cmax are provided. The first example considers a highly coherent sensing matrix with a large Cmax (100). The original sensing matrix is comprised of 1000 highly coherent columns, which are divided into 10 similar column groups with equal length of 100. While in the second and third examples, less coherent sensing matrices are adopted with Cmax as 10 and 2, respectively.
5
Oð108 Þ
Example 2
Oð10 Þ
Oð10 Þ
Oð108 Þ
Example 3
Oð105 Þ
Oð105 Þ
Oð108 Þ
8
The parameters of the three examples are presented in Table 3. The last column of Table 3 lists the values of C KHcc , which is proportional to the computational complexity of the proposed algorithm. It can be seen that the more coherent the sensing matrix is, the larger the C KHcc will be. The computational complexity of the three algorithms for the three specific examples is listed in Table 4. When the sensing matrix is highly coherent with a large Cmax (100, in the first example), the computational complexity of the proposed algorithm is higher than that of the SP and BP algorithms. On the contrary, if the sensing matrix is less coherent with a small Cmax (2, in the third example), the computational complexity of the proposed algorithm is comparable to that of the SP algorithm and lower than BP algorithm. 4. Conclusions In this paper, a novel algorithm called the similar sensing matrix pursuit (SSMP) is proposed to reconstruct a K-sparse vector based on a deterministic sensing matrix with high coherence. The proposed algorithm consists of two parts: the off-line processing and online processing. The goal of the off-line processing is to find a similar
J. Liu et al. / Signal Processing 95 (2014) 101–110
compact sensing matrix, which contains as much information as possible from the original sensing matrix. The online processing consists of a rough estimation process and a refined estimation process. In the rough estimation process, a subspace pursuit algorithm is used to find a rough estimate of the true support set, which contains the indices of the columns that contribute to the original sparse vector. Three kinds of structures of the estimated support set are considered, and three individual refined estimation processes are carried out under these three conditions. From the simulation results, it can be seen that the proposed algorithm obtains much better performance when coping with the deterministic sensing matrix with high coherence, compared with the SP and BP algorithms.
Acknowledgments This work was supported in part by National 973 project of China through Grant 2013CB329405, the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (61221063), and the Natural Science Table A1 Notations. Variable
Definition
K x N y M Φ δK μðΦÞ φi H λðφi ; φj Þ G T1
Sparsity level K-sparse vector Length of K-sparse vector x Measurement vector Length of measurement vector y Sensing matrix Restricted isometry constant (RIC) Coherence of Φ ith column of Φ Additional sensing matrix Similarity between φi ; φj Gram matrix Threshold set to distinguish the highly coherent columns from the incoherent columns Threshold set to distinguish the highly coherent columns inside a similar column group Similarity distance between φi ; φj
T2 dsimilar ðφi ; φj Þ Ks Γi Shc Sic ιi ζi Nhc Nic D γ ij N Γi Ψ γ iC
A constant in calculating dsimilar The ith similar column group Set of highly coherent columns Set of incoherent columns The ith highly coherent column The ith incoherent column Number of highly coherent columns Number of incoherent columns Number of similar column groups The jth column in the ith similar column group Γi Number of columns in Γi Similar compact sensing matrix
The condensed column from Γi in the similar compact sensing matrix Ψ γ The noise magnitude mr The measurement residual in the OMP algorithm ε An appropriately chosen bound on the noise magnitude The estimated signal x^ O 1 2 K a^ ¼ fa^ ; a^ ; …; a^ g Estimated support set The first element of a^ corresponds to γ iC , the ^a 1 ðγ iC Þ index of the ith condensed column.
109
Table A2 Notations. b^ ¼ fΓ i ; Γ j ⋯; Γ l g A set containing the indices of K similar column groups corresponding to a^ The indices of all the columns in the selected f^ similar column groups in b^ Hcc Hci ϒp
Total number of columns in the K contributing similar column groups Total number of columns from the V contributing similar column groups The pth candidate support set
Υ 1p ðγ i1 Þ
The first element of ϒp corresponding to γ i1
Nco p x^
Total number of combinations The estimate of the original sparse vector based on ϒp
p x^ ϒp p x^ f1;2;…;NgΥ p
p The entries of x^ indexed by iA ϒp p The entries of x^ indexed by iA f1; 2; …; Ngϒp
Φpϒ † x^ min
The columns of Φ indexed by iA ϒp Pseudo-inverse operation The estimate with the least residual The first element of a^ corresponds to ζ i , the ith incoherent column. Number of selected condensed columns The pth combination
1 a^ ðζ i Þ
V rp ¼ fγ i1 ; …; γ j1 g Ncb O(LS) Cmax Mantenna Ngrid Δθ θ0 Φsteering u d λ r S
Number of combinations Order of computational cost for least squares algorithm Size of the largest similar column group Number of antennas Number of fine grids in the DOA angle plane Size of each cell Initial angle of the DOA plane Steering matrix M antenna 1 steering vector Distance between the elements of the arrays Wavelength Reflection coefficient Support set
Foundations of China (nos. 61104051, 61074176 and 61004087), the Fundamental Research Funds for the Central Universities, and the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry. Appendix A. Notations The definitions of variables used in the paper are listed in Tables A1 and A2. References [1] J. Romberg, Imaging via compressive sampling [introduction to compressive sampling and recovery via convex programming], IEEE Signal Processing Magazine 25 (2) (2008) 14–20. [2] Xue Bi, Xiang-dong Chen, Yu Zhang, Bin Liu, Image compressed sensing based on wavelet transform in contourlet domain, Signal Processing 91 (5) (2011) 1085–1092. [3] C. Berger, S. Zhou, J. Preisig, P. Willett, Sparse channel estimation for multicarrier underwater acoustic communication: from subspace methods to compressed sensing, IEEE Transactions on Signal Processing 58 (3) (2010) 1708–1721. [4] W.U. Bajwa, J.D. Haupt, A.M. Sayeed, R.D. Nowak, Compressed channel sensing: a new approach to estimating sparse multipath channels, Proceedings of the IEEE 98 (6) (2010) 1058–1076. [5] M. Herman, T. Strohmer, Compressed sensing radar, IEEE Transactions on Signal Processing 13 (10) (2006) 589–592. [6] Guanghui Zhao, Zhengyang Wang, Qi Wang, Guangming Shi, Fangfang Shen, Robust ISAR imaging based on compressive sensing from noisy measurements, Signal Processing 92 (1) (2012) 120–129.
110
J. Liu et al. / Signal Processing 95 (2014) 101–110
[7] C. Berger, B. Demissie, J. Heckenbach, P. Willett, S. Zhou, Signal processing for passive radar using OFDM waveforms, IEEE Journal of Selected Topics in Signal Processing 4 (1) (2010) 226–238. [8] Joachim H.G. Ender, On compressive sensing applied to radar, Signal Processing 90 (5) (2010) 1402–1414. [9] Y. Yu, A. Petropulu, H. Poor, Measurement matrix design for compressive sensingcbased MIMO radar, IEEE Transactions on Signal Processing 59 (11) (2011) 5338–5352. [10] E. Candès, J. Romberg, T. Tao, Stable signal recovery from incomplete and inaccurate measurements, Communications on Pure and Applied Mathematics 59 (8) (2006) 1207–1223. [11] D. Donoho, Compressed sensing, IEEE Transactions on Information Theory 52 (4) (2006) 1289–1306. [12] E. Candès, Compressive sampling, Proceedings of the International Congress of Mathematicians 3 (2006) 1433–1452. [13] E. Candès, J. Romberg, T. Tao, Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information, IEEE Transactions on Information Theory 52 (2) (2006) 489–509. [14] E. Candès, T. Tao, Decoding by linear programming, IEEE Transactions on Information Theory 51 (12) (2005) 4203–4215. [15] S. Mallat, Z. Zhang, Matching pursuit with time-frequency dictionaries, IEEE Transactions on Signal Processing 41 (12) (1993) 3397–3415. [16] J.A. Tropp, Greed is good: algorithmic results for sparse approximation, IEEE Transactions on Information Theory 50 (10) (2004) 2231–2242. [17] D. Needell, R. Vershynin, Uniform uncertainty principle and signal recovery via regularized orthogonal matching pursuit, Foundations of Computational Mathematics 9 (2009) 317–334. [18] D.L. Donoho, Y. Tsaig, I. Drori, J.-L. Starck, Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit, IEEE Transactions on Information Theory 58 (2) (2012) 1094–1121. [19] W. Dai, O. Milenkovic, Subspace pursuit for compressive sensing signal reconstruction, IEEE Transactions on Information Theory 55 (5) (2009) 2230–2249. [20] D. Needell, J.A. Tropp, COSAMP: iterative signal recovery from incomplete and inaccurate samples, Applied and Computational Harmonic Analysis 26 (3) (2009) 301–321.
[21] H. Huang, A. Makur, Backtracking-based matching pursuit method for sparse signal reconstruction, IEEE Signal Processing Letters 18 (7) (2011) 391–394. [22] J. Bourgain, S. Dilworth, K. Ford, S. Konyagin, D. Kutzarova, Explicit constructions of RIP matrices and related problems, Duke Mathematics Journal 159 (1) (2011) 145–185. [23] A. Cohen, W. Dahmen, R. DeVore, Compressed sensing and best k-term approximation, Journal of American Mathematical Society 22 (2009) 211–231. [24] S.D. Howard, A.R. Calderbank, S.J. Searle, A fast reconstruction algorithm for deterministic compressive sensing using second order reed-muller codes, in: IEEE Conference on Information Sciences and Systems (CISS2008), 2008. [25] A. Amini, F. Marvasti, Deterministic construction of binary, bipolar, and ternary compressed sensing matrices, IEEE Transactions on Information Theory 57 (4) (2011) 2360–2370. [26] E. Candès, J. Romberg, Sparsity and incoherence in compressive sampling, Inverse Problems 23 (3) (2007) 969–985. [27] J. Zhang, D. Zhu, G. Zhang, Adaptive compressed sensing radar oriented toward cognitive detection in dynamic sparse target scene, IEEE Transactions on Signal Processing 60 (4) (2012) 1718–1729. [28] M.F. Duarte, Y.C. Eldar, Structured compressed sensing: from theory to applications, IEEE Transactions on Signal Processing 59 (9) (2011) 4053–4085. [29] D.L. Donoho, M. Elad, V.N. Temlyakov, Stable recovery of sparse overcomplete representations in the presence of noise, IEEE Transactions on Information Theory 52 (1) (2006) 6–18. [30] M. Elad, Optimized projections for compressed sensing, IEEE Transactions on Signal Processing 55 (12) (2007) 5695–5702. [31] A. Gurbuz, J. McClellan, V. Cevher, A compressive beamforming method, in: Proceedings of IEEE International Conference on Acoustics and Speech Signal Processing, 2008, pp. 2617–2620. [32] I. Stojanovic, W. Karl, M. Cetin, Compressed sensing of mono-static and multi-static SAR, in: SPIE Defense and Security Symposium, Algorithms for Synthetic Aperture Radar Imagery XVI. [33] S. Chen, D.L. Donoho, M. Saunders, Atomic decomposition by basis pursuit, SIAM Journal of Science Computing 20 (1998) 33–61. [34] SparseLab [Online]. Available: 〈http://dsp.rice.edu/cs〉.