Splitting Matching Pursuit Method for Reconstructing Sparse Signal in Compressed Sensing

In this paper, a novel method named as splitting matching pursuit (SMP) is proposed to reconstructK-sparse signal in compressed sensing.The proposedmethod selectsFl (Fl > 2K) largest components of the correlation vector c, which are divided intoF split sets with equal length l.The searching area is thus expanded to incorporatemore candidate components, which increases the probability of finding the true components at one iteration. The proposed method does not require the sparsity level K to be known in prior. The Merging, Estimation and Pruning steps are carried out for each split set independently, which makes it especially suitable for parallel computation. The proposed SMP method is then extended to more practical condition, e.g. the direction of arrival (DOA) estimation problem in phased array radar systemusing compressed sensing. Numerical simulations show that the proposedmethod succeeds in identifying multiple targets in a sparse radar scene, outperforming other OMP-type methods. The proposed method also obtains more precise estimation of DOA angle using one snapshot compared with the traditional estimation methods such as Capon, APES (amplitude and phase estimation) and GLRT (generalized likelihood ratio test) based on hundreds of snapshots.


Introduction
The standard noiseless model in compressed sensing is where  ∈ R  is a -sparse signal ( ≪ ),  ∈ R  is a measurement of , and Φ is an  ×  sensing matrix.The compressed sensing recovery problem is defined as follows: given  and Φ, find a signal  within the class of interest satisfies (1) exactly.The compressed sensing recovery process consists of a search for the sparsest signal  that yields the measurement .By defining the  0 "norm" of a vector ‖‖ 0 as the number of nonzero entries in , the simplest way to pose a recovery algorithm is using the optimization Solutions to (2) therefore lead to algorithms for recovering -sparse signals from  linear measurements.In general, the minimization of ( 2) is NP-hard.An alternative to the  0 "norm" used in (2) is to use the  1 "norm", defined as ‖‖ 1 = ∑  =1 |()|.The resulting adaptation of (2), known as basis pursuit (BP) [1], is formally defined as Since the  1 "norm" is convex, (3) can be seen as a convex relaxation of (2).The optimization (3) can be modified to allow for noise in the measurements  = Φ + , where  denotes an  × 1 measurement noise vector.We simply change the constraint on the solution to where  > ‖‖ 2 is an appropriately chosen bound on the noise magnitude.This modified optimization is known as basis pursuit with inequality constraints (BPIC) and is a quadratic program with polynomial complexity solvers [2].The Lagrangian relaxation of this quadratic program is written as  = arg min ‖‖ 1 +       − Φ    2 (5) and is known as basis pursuit denoising (BPDN).There exist many efficient solvers to find BP, BPIC, and BPDN solutions; for an overview, see [3].Unfortunately, the complexity of the linear programming algorithms for solving (3)∼( 5) is highly impractical for large-scale applications.
An alternative approach to sparse signal recovery is based on the idea of iterative greedy pursuit.The basic greedy algorithm is the matching pursuit (MP) [4].OMP [5] is a variation of MP method, which adds a least-squares minimization step to MP method to obtain the best approximation over the chosen atoms.Unlike MP and OMP choosing just one atom at each time, ROMP, StOMP, SP, CoSaMP, and BAOMP methods choose several atoms at each iteration.Furthermore, regularized OMP (ROMP) [6], subspace pursuit (SP) [7], compressive sampling matching pursuit (CoSaMP) [8], and backtracking-based matching pursuit (BAOMP) [9] methods also use a two-step selection technique to carefully choose the atoms.While SP and CoSaMP have offered comparable theoretical reconstruction quality to the linear programming methods along with low reconstruction complexity, they require the sparsity level to be known for exact recovery.As an improvement, the BAOMP algorithm achieves the blind sparse signal reconstruction without requiring the sparsity level .However, the BAOMP algorithm adopts two parameters (atom-adding constant  1 and atom-deleting constant  2 ) which are related with sparsity level and required to be tuned online.
In this paper, a novel method named as splitting matching pursuit (SMP) is proposed to reconstruct sparse signal in compressed sensing.The SP/CoSaMP algorithm assumes that the true indices in the support set of the original signal correspond to the large components of the correlation vector and choose /2 largest components at each iteration.However, in practice some true indices may correspond to a set of small components.The proposed method selects  ( > 2) largest components of the correlation vector, which expands the searching area and increases the probability of finding the true components at one iteration.The  candidate components are divided into  split sets with equal length .The proposed method does not require the sparsity level  to be known in prior.Different from BAOMP algorithm which adopts two tuning parameters related with sparsity level, the proposed algorithm uses two parameters  and  which are preset and determined by  and .The candidate components are divided into  split sets, which could be processed simultaneously and suitable for parallel computation.
The proposed SMP method is then extended to more practical condition, for example, the direction of arrival (DOA) estimation problem in phased-array radar system using compressed sensing.Numerical simulations show that the proposed method succeeds in identifying multiple targets in a sparse radar scene, outperforming other OMP-type methods.The proposed method also obtains more precise estimation of DOA angle using one snapshot compared with the traditional estimation methods such as Capon [10], APES [11], and GLRT [12] based on hundreds of snapshots.
The rest of the sections are organized as follows.The proposed SMP method is introduced in Section 2. The simulation results are listed in Section 3, and the paper is summarized in Section 4.

The Splitting Matching Pursuit Method
The most difficult part of signal reconstruction is to identify the locations of  largest components in the target signal.The SP/CoSaMP algorithm assumes that the true indices in the support set of the original signal correspond to the large components of the correlation vector and choose /2 largest components at each iteration.However, in practice some true indices may correspond to a set of small components.The proposed SMP selects  ( > 2) largest components of the correlation vector, which expands the searching area and increases the probability of finding the true components at one iteration.In the following, a short introduction of the proposed method is given in Section 2.1, and the detailed algorithm is introduced in Section 2.2.The convergency analysis, parameters setting, complexity analysis, and convergency speed analysis of the proposed method are provided in Sections 2.3 to 2.6, respectively.

Brief Introduction to Splitting Matching Pursuit Method.
A schematic diagram of the proposed algorithm is depicted in Figure 1.At the beginning of each iteration, the proposed method selects  ( > 2) largest components of the correlation vector, which are divided into  split sets with equal length .Each split set is merged with the estimated support set from the previous iteration, resulting in a sequence of merged split sets.The proposed algorithm then approximates the target signal on each merged split set to obtain a sequence of split estimates using least-square algorithm.A set of pruned split sets are then built by retaining only the  largest magnitude entries in the split estimates.The obtained pruned split sets are then combined to a final merged set, which contains as much as possible true components.An interim estimate is then calculated based on the final merged set using least-square algorithm.An estimated support set is obtained by retaining  indices corresponding to the largest magnitude entries in the interim estimate, based on which a final estimate is generated.The iterations repeat if the  2 norm (magnitude) of the calculated residual is less than a threshold .

Detailed Procedures of Splitting Matching Pursuit Method.
The detailed procedure of the SMP method is listed as follows.Algorithm 1 (the SMP method).Input.Sensing matrix Φ, measurement vector , parameters  and , and threshold  to halt the iterations.
Output.The estimated signal  out .Initialization (1)  0 = 0, where  0 indicates the initial estimated support set and 0 denotes empty set. is the estimated support set with length .
Iteration.At the th iteration, go through the following steps.
(2) Support Merging.Each newly identified split set is united with the estimated support set from the previous iteration, resulting in a sequence of merged split sets as where    denotes the th merged split set at the th iteration, and  −1 denotes the estimated support set at the ( − 1)th iteration.
(3) Estimation.The proposed algorithm then solves a least square problem to approximate the nonzero entries of the target signal on each merged split set (   ,  = 1, . . ., ) and sets other entries as zero, resulting in a sequence of split estimates as where    denotes the th split estimate of the target signal at the th iteration.
where    denotes the th pruned split set at the th iteration.
(5) Split Sets Merging.The pruned split sets are merged to form a final merged set,   , as (6) Estimation.An interim estimate of the original signal is calculated based on the final merged set   using least-square algorithm as where   denotes the interim estimate of the original signal at the th iteration.
where   denotes the estimated support set at the th iteration.(8) Estimation.Obtain the final estimate at each iteration, based on   using least-square algorithm as where    denotes the final estimate at the th iteration.
where   denotes the residual at the th iteration.(10) If ‖  ‖ 2 > , perform the correlation calculation   = Φ *   , and then go to step (1) of the ( + 1)th iteration; otherwise, set  out =    and quit the iteration.

Convergency Analysis.
Here, we will discuss the convergency of the proposed SMP method.
Theorem 2 (Theorem 2.1 in [8]).Let  ∈ R  be a -sparse signal, and let its corresponding measurement be  = Φ +  ∈ R  .If the sampling matrix satisfies the restricted isometry property (RIP) with constant In particular, where V denotes the unrecoverable energy in the signal.
Proposition 3. Let  ∈ R  be a -sparse signal, and let its corresponding measurement be  = Φ +  ∈ R  .If the sampling matrix satisfies the RIP with constant then the proposed SMP algorithm is guaranteed to recover  from  via a finite number of iterations.
Proof.The proving process is very similar to that of Theorem 2 (Theorem 2.1 in [8]).The CoSaMP algorithm selects the 2 largest components of the correlation vector at each iteration, while the proposed SMP method selects the  largest components of the correlation vector.We can obtain the similar results through replacing 2 with  in the derivation process in [8].

Parameters Setting.
Here, we will discuss how to set values of the number of split sets  and the length of the split set .
Proof.Considering the monotonicity of   : for any two integers  ≤   ,   ≤    according to Lemma 1 in [7].So we have   <  4 provided that  < 4.
Proof.As the estimated support set,   contains at least  elements, resulting in  ≤ .
According to Propositions 3∼5, in order to guarantee a perfect recovery of the -sparse vector  from  via a finite number of iterations, we have  < 4 ≤ 4, resulting in  < 4. Since  is set large enough to expand the searching area,  is set as 3 in the proposed method.We then have  < (4/3) according to  < 4.Propositions 3∼5 are based on Theorem 2, which has a rigid setting for  4 .We can relax the range of  to [, 2).As a result, we need not know the exact value of  and can select an integer randomly from [, 2) based on an estimated value of .

Complexity Analysis.
The process of complexity analysis is similar to that of [7].In each iteration, the correlation maximization procedure requires  computations in general, while the cost of computing the projections is of the order of ((2) 2  + (2) 2 ), which could be approximated as ( 2 ) when  is chosen as a small value (e.g. 3 in the simulation setup).As a result, the total cost of computing is of the order of ( +  2 ), which is comparable to the SP/CoSaMP algorithm.

Convergence Speed Analysis.
Here, we will discuss the convergence speed of the proposed SMP method.
Theorem 6 (Theorem 8 in [7]).The number of iterations ( SP  ) of the SP algorithm is upper bounded by where Proposition 7. The number of iterations ( SMP  ) of the SMP algorithm is upper bounded by Proof.The SP algorithm selects the  largest components of the correlation vector at each iteration, while the proposed SMP method selects the  ( ≥  according to Proposition 5) largest components of the correlation vector.The searching area is thus expanded to at least  times of that of SP algorithm.The probability of finding the true components at one iteration is increased to at least  times of that of SP algorithm.As a result, the number of iterations is decreased to 1/ of that of the SP algorithm by average according to Theorem 6.

Simulation Results and Analysis
In this section, a simple example is firstly carried out to verify the performance of the proposed SMP method in reconstructing zero-one binary and Gaussian signals.The proposed method is then extended cope with the DOA estimation problem in phased-array radar system.

A Simple Example.
The simulations are carried out to compare the accuracy of different reconstruction algorithms empirically.The proposed SMP algorithm is compared with some popular greedy algorithms (including OMP, ROMP, StOMP, SP, CoSaMP, and BAOMP algorithms) and the convex optimization algorithm (BP method).
In the simulation setup, a signal sparsity level  is chosen such that  ≤ /2 given the length of the original signal  ( = 256) and the length of the measurement vector  ( = 128).An  ×  sampling matrix Φ is randomly generated from the standard i.i.d.Gaussian ensemble.A support set  of size  is selected uniformly at random, and the original sparse signal vector  is chosen as either Gaussian signal or zero-one signal [7].The estimate of the original signal, x, is computed based on the measurement vector  generated through  = Φ.In the experiments, OMP uses  iterations, StOMP and BP methods use the default settings (OMP, StOMP, and BP tools use SparseLab [13]), and ROMP and SP methods use the parameters given in [6,7], respectively.The proposed SMP method use  = 3,  = 84,  max = ,  = 10 −5 , as the input parameters.
The signal sparsity level  is varied from 0 to /2.For each fixed , five hundred Monte Carlo simulations are carried out for each algorithm.The reconstruction is considered to be exact when the  2 norm of the difference between the original signal  and the reconstructed one x is smaller than 10 −5 , that is, ‖ − x‖ 2 < 10 −5 .The frequency of exact reconstruction () is used to evaluate the reconstruction performance of the different methods, which is defined as where  denotes the number of exact reconstructions for each algorithm given a fixed , and  MC denotes the number of Monte Carlo simulations.The frequency of exact reconstruction is also adopted by the SP method to evaluate the reconstruction performance [7].Figures 2 and 3 show the reconstruction results for binary zero-one and Gaussian sparse signals, respectively.We only present the results of the SP algorithm since the SP and CoSaMP algorithms are almost the same with different deviation process, and both obtain the same simulation results.As can be seen in Figure 2, for binary zero-one sparse signal which is a difficult case for OMP-type methods, the performance of the proposed SMP method is much better than all other OMP-type methods and comparable to the BP minimization method.Of particular interest is the sparsity level at which the recovery rate drops below 100%, that is, the critical sparsity defined in [7].It could be seen from Figure 2 that the proposed method is with the largest critical sparsity (39), which exceeds that of the BP method (36).For the Gaussian sparse signal, as shown in Figure 3, the proposed method also gives the comparable performance to the BP method.

DOA Estimation Based on Splitting Matching Pursuit
Method.In this section, the proposed SMP method is extended to solve the DOA estimation problem in phasedarray radar system.The signal model for DOA estimation in phased-array radar system is represented in a standard compressed sensing form in Section 3.2.1,where the sparse radar scene is abstracted as a sparse signal.And the simulation results are shown in Section 3.2.2, which shows that the proposed method succeeds in identifying multiple targets in a sparse radar scene.

Signal Model for DOA Estimation and Sparse Representation.
Assume that a phased-array radar system consists of half wavelength spaced uniform linear arrays (ULAs).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 [14]).We assume that the other parameters like range and Doppler frequency have been isolated before by appropriate processing.
The ULA of the phased-array radar system consists of  antennas, which are used to emit the transmitted signal ().The  × 1 received complex vector of array observations is defined as () = [ 1 (), . . .,   ()]  .Assuming a hypothetical target located at a DOA angle of  in the far field, the received complex vector of array observations can be written as where () is the reflection coefficient of the hypothetical target, and () is an  × 1 complex Gaussian noise vector.() is the  × 1 steering vector, which is defined as where  is the distance between the elements of the arrays and  denotes wavelength.Assuming  targets are observed with reflection coefficients {  }  =1 and DOA angles {  }  =1 , the  × 1 received complex vector of array observations can be written as where () is an  × 1 complex Gaussian noise vector.Equation (24) could be rewritten as where () = [( 1 )( 2 ) ⋅ ⋅ ⋅ (  )] is an  ×  steering matrix, and (, ) = ()[( 1 )( 2 ) ⋅ ⋅ ⋅ (  )]  denotes a  × 1 reflection vector.
Since the radar scene is generally in practice sparse, compressed sensing is a valid candidate for estimating the DOA angles for multiple targets.To do so, the DOA angle plane is divided into  fine grids, each cell generally with the same size Δ.The th grid represents the DOA angle of  0 + ( − 1)Δ, where  0 is the initial angle of the DOA plane.The steering matrix and reflection vector in (25) are extended to obtain the  ×  extended steering matrix Φ and the  × 1 extended reflection vector , which are defined as Φ = [( 0 )( 0 +Δ) ⋅ ⋅ ⋅ ( 0 + ( − 1)Δ)( 0 + ( − 1)Δ)] and  = ()[( 0 )( 0 + Δ) ⋅ ⋅ ⋅ ( 0 + ( − 1)Δ)( 0 + ( − 1)Δ)]  .Since small number of grids are occupied by the targets,  is a sparse vector with the th element defined as () = ()( 0 + ( − 1)Δ) if the th grid is occupied by the target; otherwise, () = 0.As a result, the  × 1 received complex vector of array observations  could be written as where  is an  × 1 complex Gaussian noise vector.Though in (26) the radar vectors and matrices are complex valued in contrary to the original compressed sensing environment, it is easy to transfer it to real variables according to [15,16].
Discussion.In [17,18], it is assumed that the discretized step is small enough so that each target falls on some specific grid point.However, no matter how finely the parameter space is gridded, the sources may not lie in the center of the grid cells, and consequently there is mismatch between the assumed and the actual bases for sparsity.The sensitivity of compressed sensing to mismatch between the assumed and the actual sparsity bases is studied in [19].The effect of basis mismatch is analyzed on the best -term approximation error, and some achievable bounds for the  1 error of the best -term approximation are provided.The readers can refer to [19] for a detailed analysis on the influence of the griding operations on the estimation performance.

Simulation Results.
In this section, an example about DOA estimation is provided based on the phased-array radar system, which consists of half wavelength spaced uniform linear arrays (ULAs).The number of transmit/receive antennas is 20.The antennas transmit independent orthogonal quadrature phase shift keyed (QPSK) waveforms and the carrier frequency is 8.62 GHz.The SNR of the measurement noise is set to a fixed value (20 dB).The range of the DOA plane is [0 ∘ , 90 ∘ ], which is divided into 30 cells with the initial angle ( 0 ) and angle interval (Δ) equaling 0 ∘ and 3 ∘ , respectively.A maximum of  = 512 snapshots are considered at the receive node.Targets may appear at directions represented by DOA angles.The proposed SMP method is compared to several popular greedy algorithms, for example, the OMP, ROMP, StOMP, SP, CoSaMP and BAOMP algorithms, and the convex optimization algorithm, BPDN method.The proposed method is further compared to three commonly used methods in DOA estimation, for example, the Capon, APES, and GLRT methods.The compressed sensing based methods (the SMP method, the greedy algorithms, and the BPDN method) use one snapshot only, and the Capon, APES, and GLRT methods use 512 snapshots each.Five hundred Monte Carlo simulations are carried out, and in each trial four targets locate randomly within the DOA range of [0 ∘ , 90 ∘ ], and the corresponding reflection coefficients are set as {  = 1,  = 1, . . ., 4}.
The average reconstruction error is adopted to evaluate the reconstruction performance of the methods, which is defined as where   denotes the reconstruction error at the th Monte Carlo simulation, which is defined as where   and   estimate represent the true and estimated signal representing the sparse radar scene at the th Monte Carlo simulation, respectively.The average root mean square error (RMSE) is also adopted to evaluate the DOA estimation performance of the methods, which is defined in [20].
The results in Table 1 show that the proposed SMP method is with the smallest average reconstruction error and average RMSE.The proposed method succeeds in identifying multiple targets in a sparse radar scene, outperforming other OMP-type methods.It also obtains more precise estimation of DOA angle using one snapshot compared with the traditional estimation methods such as Capon, APES, and GLRT based on 512 snapshots.

Conclusion
We have presented a novel SMP method for sparse signal reconstruction in compressed sensing.The proposed method expands the searching area and increases the probability of finding the true components at one iteration.It also does not require the sparsity level  to be known in prior.The proposed method is then extended to more practical condition, for example, the direction of arrival (DOA) estimation problem in phased-array radar system using compressed sensing.Numerical simulations show that the proposed method succeeds in identifying multiple targets in a sparse radar scene, outperforming other OMP-type methods.The proposed method also obtains more precise estimation of DOA angle using one snapshot compared with the traditional estimation methods such as Capon, APES, and GLRT based on hundreds of snapshots.

Figure 1 :
Figure 1: Description of reconstruction procedures of the SMP method.
The vector   |   is composed of the entries of   indexed by  ∈   , and   | (  )  is composed of the entries of   indexed by  ∈ (  )  .(7) Pruning.Obtain the estimated support set by retaining  indices corresponding to the largest magnitude entries in the vector   |   , as   = {indices of  largest magnitude entries in         } , The vector    |   is composed of the entries of    indexed by  ∈   , and    | (  )  is composed of the entries of    indexed by  ∈ (  )  .The matrix Φ   consists of the columns of Φ with indices  ∈   .(9) Residual calculation: