Adaptive matching pursuit with constrained total least squares

Compressive sensing (CS) can effectively recover a signal when it is sparse in some discrete atoms. However, in some applications, signals are sparse in a continuous parameter space, e.g., frequency space, rather than discrete atoms. Usually, we divide the continuous parameter into finite discrete grid points and build a dictionary from these grid points. However, the actual targets may not exactly lie on the grid points no matter how densely the parameter is grided, which introduces mismatch between the predefined dictionary and the actual one. In this article, a novel method, namely adaptive matching pursuit with constrained total least squares (AMP-CTLS), is proposed to find actual atoms even if they are not included in the initial dictionary. In AMP-CTLS, the grid and the dictionary are adaptively updated to better agree with measurements. The convergence of the algorithm is discussed, and numerical experiments demonstrate the advantages of AMP-CTLS.


Introduction
A new class of techniques called Compressed Sampling or Compressive Sensing (CS) has been widely used recently, due to the fact that CS techniques have shown good performance in different areas such as signal processing, communication and statistics; see, e.g., [1].Generally, CS finds the sparsest vector x from measurements y = Φx, where Φ is often referred to as dictionary with more columns than rows, and each column of the dictionary is called an atom or a basis.
Matching Pursuit (MP) is a set of popular greedy approaches to Compressive Sensing.The basic idea is to sequentially find the support set of x and then project on the selected atoms.The atoms selected in the support set are mainly determined by correlations between atoms and the regularized measurements [2].
These MP methods [2][3][4][5][6][7] do not consider the off-grid problem in grid-based CS approaches.In some applications of CS, such as harmonic retrieval and radar signal processing (e.g.range profiling [8,9], direction of arrival estimate [10][11][12]), we usually divide a continuous parameter space into discrete grid points to generate the dictionary.For example, in harmonic retrieval, frequency space is divided and dictionary is a discrete Fourier transform (DFT) matrix.The off-grid problem emerges when the actual frequencies are placed off the predefined grid.The mismatch between the predefined and actual atom can lead to performance degradation in sparse recovery (e.g., [13][14][15]).
The grid misalignment problem in CS has recently received growing interest.The sensitivity of CS to the mismatch between the predefined and actual atoms is studied in [13]; however, the focus of that paper is mainly on mismatch analysis rather than development of an algorithm.Cabrera et al. [16] and Zhu et al. [14] respectively provided an Iterative Re-Weighted (IRW)-based and a Lasso-based method to recover an unknown vector considering the atom misalignment, whereas we focus on MP methods in this paper.Compared with IRW or Lasso, MP methods greedily find the support set and greatly reduce the dimension of the CS problem; thus, they have an advantage in computability.Gabriel [17] proposed best basis compressive sensing in a tree-structured dictionary, but some dictionaries (e.g., DFT matrix) do not possess a tree structure.
To alleviate the off-grid problem in Matching Pursuit, we developed Adaptive Matching Pursuit with Constrained Total Least Squares (AMP-CTLS).In AMP-CTLS, we model the grid as an unknown parameter, and adaptively search for the best one.We choose harmonic retrieval to demonstrate the performance of AMP-CTLS.The algorithm can also be applied to jointly estimate range and velocity in randomized step frequency (RSF) radar.Note that in the RSF scenario range-velocity estimation is hard to be directly solved by subspace-based methods, e.g.Capon's method, MUSIC and ESPRIT [18].Since only one snapshot data is available in RSF radar, to obtain the covariance matrix these subspace-based methods need to apply smoothing method, which requires uniform and linear condition [18].However, this condition is not satisfied in the case of random frequency model.This paper is structured as follows.Section 2 introduces grid-based CS and outlines the procedures of AMP-CTLS.In Section 3 and 4, we discuss the implementation of AMP-CTLS in harmonic retrieval and RSF radar, respectively.In Section 5, numerical examples are presented to illustrate merits of AMP-CTLS.denotes the expectation of a random variable.

Grid-based CS and the AMP-CTLS Algorithm
The signal model of grid-based CS is introduced in Subsection 2.1.We combine the greedy idea of MP methods and the Constrained Total Least Squares (CTLS) technique [20], and thus produce AMP-CTLS to alleviate the off-grid problem.In AMP-CTLS, the grid is cast as an unknown parameter, and is jointly estimated together with x.In Subsection 2.2, the framework of AMP-CTLS is given.Subsection 2.3 is dedicated to the Iterative Joint Estimator (IJE) algorithm, which is implemented in AMP-CTLS.In the IJE algorithm, the CTLS technique is used, which is presented in Subsection 2.4.Subsection 2.5 summarizes the entire procedure of AMP-CTLS.In Subsection 2.6, the convergence of IJE is analyzed.

Grid-based CS
CS promises efficient recovery of sparse signals.In many applications, signals are sparse in a continuous parameter space rather than finite discrete atoms.Usually, we divide the continuous parameter into discrete grid points and cast the problem as a grid-based CS model: where y ∈ C M×1 and w ∈ C M×1 are measurement vector and white Gaussian noise (WGN) vector, respec- , and the mapping g → Φ is known.For example, to recover a frequency sparse signal, we grid the frequency space into discrete frequency points g Φ is a DFT matrix, of which the mth-row, nth-column element is exp j2π n N m .However, the signal is only sparse in the DFT atoms if all of the sinusoids are exactly at the pre-defined grid points [13].In some cases, no matter how densely we grid the frequency space, the sinusoids could be off-grid, which saps the performance of CS methods [13].

Main Idea of AMP-CTLS
The off-grid problem usually emerges because we do not often have enough priori knowledge to generate a perfect grid to guarantee that all of the signals exactly lie on grid points.Thus, we cast the grid as an unknown parameter, and search for the best grid g as well as the sparsest x by solving the optimum problem: x, ĝ = arg min where η is the noise power.Equ. ( 2) is similar to that used in traditional MP methods [2][3][4][5][6][7], except that we recover x and simultaneously estimate the grid.In most cases, solving (2) is a complex non-linear optimum problem.In this paper, an iterative method is introduced.
AMP-CTLS inherits the greedy idea from MP methods, which use correlations to iteratively find the support set.In each iteration, one or more atoms are added into the support set.Suppose the support set is obtained as Λ (k) after the kth iteration, and denote the corresponding grid points as g Λ .In traditional MP methods [2][3][4][5][6][7], x Λ is estimated by solving a least squares problem.In AMP-CTLS, considering the off-grid problem, we jointly search for x Λ and the best grid points in the neighboring continuous region of g (k) Λ via (3), in which we minimize norm of the residual error, which is defined as We develop the Iterative Joint Estimator (IJE) algorithm to solve (3), which is detailed in ensuing subsection.

IJE Algorithm
It is difficult to find an analytical solution to (3).The IJE algorithm is devised to seek a numerical solution.
Given initial grid points g Λ (0), IJE searches for the best grid points g Λ in the neighborhood of g Λ (0).The mismatch of the grid is denoted as ∆g Λ = g Λ − g Λ (0) = ∆g 1 , . . ., ∆g |Λ| T .IJE includes three steps: calculate the estimation of the mismatch, ∆g Λ ; update the grid with ∆g Λ ; and estimate x Λ with projection onto the new grid points.These three steps are executed iteratively to pursue more accurate results.To distinguish from iterations in search for the support set in (3), we denote l as the counter of loops in IJE; thus, IJE is expressed as follows: x Λ (l + 1) = arg min In ( 4), CTLS technique is applied to simultaneously search for the mismatch ∆g Λ and x Λ , and ∆g Λ (l) and x CTLS are the results.C CTLS denotes the penalty function of CTLS, which is detailed in Subsection 2.4.
Since ( 6) is a linear least squares problem, the closed-form solution is The loops are terminated when the norm of residual error is scarcely reduced.

CTLS Technique
Traditional MP methods [2][3][4][5][6][7] apply least squares to calculate amplitudes of x Λ after finding the support set.When there are off-grid signals, mismatches occur in the dictionary; thus, we replace the least squares model with total least squares (TLS) criterion, which is appropriate to deal with the fitting problem when perturbations exist in both the measurement vector and in the dictionary [21].Since the dictionary mismatches are constrained by errors of grid points, we introduce the Constrained Total Least Squares (CTLS) technique [20] in AMP-CTLS to jointly estimate the grid point errors and x Λ , i.e. solving (4).It has been proved that CTLS is a constrained space state maximum likelihood estimator [20].
Suppose that we obtain the estimate of grid points as g Λ (l) after lth IJE iteration.Assume that the mismatch ∆g Λ is significantly small; thus we can approximate the perfect dictionary Φ(g Λ ) as a linear combination of the mismatch ∆g with Taylor expansion: where and o(•) denotes higher order terms.For simplicity, in this subsection we ignore the iteration counter in the notations, and R i g Λ (l) , Φ g Λ (l) are respectively simplified as R i , Φ Λ .Neglect o ∆g 2 i and the signal model in ( 1) is replaced by: CTLS models ∆g Λ as an unknown random perturbation vector.The grid misalignment and the noise vector are combined into a (M + |Λ|)-dimensional vector v = (∆g Λ ) T , w T T , and CTLS aims at minimizing It has been proved that CTLS is a constrained space state maximum likelihood estimator if v is a WGN vector [20].Thus, we first whiten v. Assume that ∆g Λ is independent of w.The covariance matrix The variance of white noise w is σ 2 w .We denote an unknown normalized vector u ∈ C (M+|Λ|)×1 as (11); thus, u is a WGN vector.
Minimize the penalty function and ( 4) is detailed as follows: The constraint condition (13) can be rewritten as: where where The equivalence between ( 13) and ( 14) is proved as follows: When W x is of full-row rank, the optimum problem (12,14) are equivalent to (17 -19), which has been proved in [20]. (17) It is quite difficult to obtain analytical solution to (17).A complex version of Newton method is developed in [20], which is presented in Appendix A. Initial value of x Λ required in Newton's method for (17) can be given as: ∆g Λ is extracted from û via ∆g Λ = D −1 0 N u, thus ( 4) is solved.The sketch of CTLS is given in Algorithm 1.As the authors' best knowledge, the convergence guarantees for this Newton method are still open question.
ALGORITHM 1: The CTLS Technique 1) Input the dictionary Φ Λ and all the coefficient matrices R i .
2) Compose W x and solve (17) with the initial value given in (20).

Sketch of AMP-CTLS
Similarly to traditional MP methods [2][3][4][5][6][7], AMP-CTLS first greedily finds the support set.Then AMP-CTLS adaptively optimizes the grid points indexed in the support set.In this paper, we imitate the greedy approach of OMP, in which only one atom is added to the support set in each iteration.If the number of atoms is known, terminate the iterations when the Cardinality of the support set reaches the pre-specified number.If it is not known, we can apply some other successfully used stopping criterions, e.g.norm of residual is below a threshold [22].A sketch of AMP-CTLS is presented in Algorithm 2.
ALGORITHM 2: The AMP-CTLS Algorithm 1) Divide the continuous parameter f into grid point ĝ(0) ; input the sparsity level K or residual threshold δ.Set the support set Λ (0) = ∅, and the residual error r (0) = y.

6) Update the residual error r
Λ , and set the elements of x not indexed in Λ to 0.

Convergence of the IJE Algorithm
Here, we analyze convergence of the IJE algorithm.Assume that the mapping g Λ → Φ (g Λ ) is linear, which means and G should be a constant matrix.
Proposition.If the measurement y is perturbed by WGN and ( 21) is obeyed, IJE monotonically reduces values of the penalty function in (3).The estimates of x Λ and g Λ satisfy: Proof.Define a penalty function as follows: thus, ∆g Λ (l) and x CTLS are obtained by solving which is the same as (4), for σ 2 w is a constant.Thus, it is satisfied that Substitute ( 5), ( 21) into f p ∆g Λ (l), x CTLS , and note that C −1 g is a positive definite matrix; thus, where the last inequality is taken from (6).The inequalities in ( 25) and ( 26) are transformed to equalities if and only if ∆g Λ (l) = 0.
For simplicity, we assume that the transform Φ(g Λ ) is linear.In some practical applications like harmonic retrieval, linearity is not strictly guaranteed.However, when atom mismatch ∆g is significantly small, the higher order errors due to Taylor expansion (8) are ignorable, and ( 21) is approximately satisfied.Numerical examples are performed in Section 5 which demostrate the convergence of the proposed algorithm in the case of harmonic retrieval.

Application in the Harmonic Retrieval
In this section, we apply AMP-CTLS in harmonic retrieval.In Subsection 3.1, the signal model of harmonic retrieval is presented and adverse effects of MP approaches [2][3][4][5][6][7] in harmonic retrieval is discussed.In Subsection 3.2, we detail the implementation of AMP-CTLS in harmonic retrieval.

Signal Model of Harmonic Retrieval
Consider a complex sinusoidal signal where y m is the mth measurement, and w m is the mth noise, m = 0, 1, . . ., M − 1.There are K sinusoids, and amplitude α k , frequency f k of the kth sinusoid are unknown parameters.When the sinusoids are sparse, i.e., K ≪ M , harmonic retrieval problem can be solved by grid-based CS approaches.Divide the digital . When all frequencies are exactly at grid points, rewrite (27) as where g n is the frequency of the nth grid point and x n = α k , the kth sinusoid is present at nth grid point, 0, no sinusoid is present at nth grid point.
Rewrite (28) in matrix form as where the mth-row, nth-column element of Φ is of the form φ(m, n) = exp(j2πg n m).Apply CS methods to seek the sparsest solution of (30).Then, estimates of the frequencies and amplitudes are obtained with the indices and magnitudes of nonzero coefficients in x, respectively.The sparsest solution can be obtained with computational MP methods, which greedily minimize the ℓ 0 norm.It can also be obtained by minimizing the ℓ 1 norm [23], the quasi-norm [24,25] or the ℓ p≤1 p-norm-like diversity [26].
We focus on MP methods in this paper for the high computation efficiency.However, conventional MP methods [2][3][4][5][6][7] suffer from performance degradation if the frequency space is not perfectly grided.When the frequency is sparsely divided, sinusoids may lie off the grid points, and accuracy of frequency estimates is limited by the gap between neighboring grid points.MP methods iteratively search for the sinusoids.
If an off-grid sinusoid emerges, the energy of this sinusoid can not be totally canceled and performs as an interference in the next iterations.The leakage of the energy may mask the weak sinusoids.On the other hand, if the frequency space is densely divided, correlations between atoms are enhanced [27], which also reduces the performance of MP methods.Especially in those MP methods that select multiple atoms into the support set in a single iteration, e.g.CoSaMP, SP, ROMP and StOMP, highly correlated atoms could be chosen in a same iteration, which impairs the numerical stability of projection onto the adopted atoms.

Harmonic Retrieval with AMP-CTLS
The AMP-CTLS algorithm can be applied for harmonic retrieval.AMP-CTLS adaptively finds the atoms and recovers the sinusoids.In those MP approaches with constant predefined atoms, frequency estimates are discrete values, depending on grid points.In AMP-CTLS frequency estimates are continuous, since estimates of the grid misalignments are continuous.In this subsection, we adjust two steps of AMP-CTLS presented in Section 2 to better fit the harmonic retrieval problem.
Calculate the R matrix in (9).According to ( 9), the mth-row, ith-column element of the R i is expressed as follows: Elements in other columns are all zeros.
Adjust the grid-updating formula in (5).In CTLS as presented in Subsection 2.4, the grid misalignment ∆g Λ is assumed to be a complex vector; therefore, the estimate ∆g Λ is complex.However, frequency grid points are restrained to be real, so regularization ∆g Λ = (∆g Λ ) * should be added to (12) in the case of harmonic retrieval.Unfortunately, the solver becomes complex, which is derived in Appendix B. For simplicity, ( 5) is replaced with (32) to approximatively update the grid points: 4 Application in RSF Radar AMP-CTLS can also be applied in randomized step frequency (RSF) radar.RSF radar can improve the range-velocity resolution and avoid range-velocity coupling problems [28,29].However, RSF radar suffers from the sidelobe pedestal problem, which results in small targets being masked by noise-like components due to dominant targets [29].Our problem of interest is to recover small targets.When the observed scene is sparse, i.e.only few targets exists, we can use sparse recovery to exploit the sparseness [9].AMP-CTLS relieves the sidelobe pedestal problem in RSF radar and recovers small targets well.
Correlation-matrix-based spectral analysis methods, e.g.MUSIC, ESPRIT [18], are hard to directly utilized in range-Doppler estimation in RSF scheme.Since only one snapshot of radar data is available and radar echoes from different scatterers are coherent, smoothing technique is invoked to obtain a full rank correlation matrix [30].Smoothing method requires that the array is uniform and linear [18].However, in RSF radar, the echoes are determined by a random permutation of integers, see (34); thus, the uniform and linear condition is not satisfied, which restricts application of correlation-matrix-based methods.We discuss a specific example of RSF radar, in which the waveform is a monotone pulse signal and the frequency of the mth pulse is f 0 + C m δf , m = 0, 1, . . ., M − 1, where f 0 is carrier frequency and δf is frequency step size.C m is a random permutation of integers from 0 to M − 1.The mth echo of radar can be expressed as (see [8,28,29]): where w m is noise in the mth echo.K denotes the number of targets and k denotes kth target.α k , p k and q k are to be learned.α k presents the scattering intensity.
where the mth-row, (c + dC)th-column element of Φ (p, q) ∈ C M×CD is s m (p c , q d ).
AMP-CTLS is implemented to solve (35).First, we find the support set Λ and then use IJE and CTLS to adjust the grid points, though CTLS described in Subsection 2.4 requires modification.The grid misalignment vector consists of two parts: Assume that ∆p Λ and ∆q Λ are independent of each other and of the noise.The covariance matrix of ∆g Λ is where In the case of RSF radar u = (D p ∆p Λ ) T , (D q ∆q Λ ) T , 1 σw w T T ∈ C (2|Λ|+M)×1 , and ) , where Since p and q are both real, formula ( 32) is used to update grid points, in which the imaginary parts of ∆p and ∆q estimates are abandoned.

Simulations
Numerical results are provided to illustrate the performance of the new algorithm.In all examples, the noise is additive Gaussian white noise.

Accuracy of AMP-CTLS
We compare the accuracy of AMP-CTLS with standard OMP.We assume that there is a single sinusoid in the measurements of form (27), where α = 1 and the signal to noise ratio SNR = α 2 /σ 2 = 5 dB, where σ 2 is the variance of noise.The number of measurements M is 32.The frequency of sinusoid is varied between two adjoining frequency grid points.The mean square errors (MSEs) of frequency estimates are calculated.
The MSEs are compared with the corresponding Cramer-Rao lower bound (CRB) [31].The frequency is uniformly divided into m grid points in xMP m (OMP M , OMP 10M , CoSaMP M , etc.) and into M points in AMP-CTLS.AMP-CTLS is configured as follows: IJE loops no more than 14 times; the normalization factors in (11) are D = I/(σ ∆f ), σ ∆f = 0.005 and σ w = 1.As shown in Fig. 1, MSEs of AMP-CTLS are close to the CRB and lower than those of OMP, except when the sinusoid is in the vicinity of the grid point.

Convergence of AMP-CTLS
We first discuss the convergence speed of the proposed IJE algorithm in noise-free case.Suppose that the sinusoid is located at f = 9.5/M , M = 32.Other conditions are the same as described in Subsection 5.1.
In the l-th iteration of IJE, we can obtain a grid point ĝ(l) with ( 5) and residual error r(l) = y − Φ (ĝ) x(l) after (6).we calculate the norm of residual r(l) 2 and the grid error |ĝ(l) − f |, and normalize the results with r(0) 2 and |ĝ(0) − f |, respectively.As shown in Fig. 2, both the residual error and the grid error converge fast (about five steps) to 0 in noiseless case.The purpose of what follows is to discuss the feasible zone of the initial grid points in noisy circumstance.
In Section 2.6, the convergence analysis of IJE is based on the assumption that the transform Φ is linear.This is only approximately satisfied in harmonic retrieval when the higher order terms of Taylor expansion (8) is ignorable, which means that the grid points indexed in the support set are required to be close to the actual frequencies.We assign SNR = 5 dB and the initial frequency grid point as g(0) = 9/M .The true frequency of the sinusoid varies from 9/M to 11/M .As shown in Fig. 2, when the distance (normalized by 1/M ) between the true frequency and the initial grid point is less than 0.7, the initial grid is adjusted to be close to the actual value, and MSEs of the frequency estimates converge to CRB.When the distance is greater than 1, the AMP-CTLS curve is close to the initial distance, which means that AMP-CTLS fails to improve the initial grid, because errors of Taylor expansion cannot be ignored and affect convergence of the algorithm.

Input of Sparsity
In Subsection 5.1 and 5.2, we assume that the sparsity K, i. e. the number of modes, is known and we use K to terminate AMP-CTLS, while a priori sparsity is not obligatory.When K is unknown, we can use norm of residual error r = y − Φ (g Λ ) x Λ as termination criterion.Furthermore, AMP-CTLS does not seriously rely on the given sparsity K ′ , and the performance is slightly affected when K ′ > K. Suppose there are three sinusoids denoted as Si 1 , Si 2 , Si 3 , where 3 /σ 2 = 10 dB.In AMP-CTLS, frequency is uniformly grided to 2M points, and other  We calculate means of the final residual norm r (K ′ ) 2 versus K ′ and present the results in Fig. 4. When all of the sinusoids have been chosen into the support set and K ′ ≥ K, energy of the sinusoids are canceled thoroughly and only noise exist in the residual.The norm of residual error becomes small and is slowly reduced along with K ′ .The results illustrate that we can use threshold of values or decrease rate of the norm of residual to end AMP-CTLS loops.
Spurious sinusoids emerge when K ′ > K. Denote the amplitude estimates by α1 , α2 , . . ., αK ′ in descend order of magnitudes and their counterparts of frequency estimates by f1 , f2 , . . ., fK ′ .MSEs of f 3 estimates versus K ′ are presented in Fig. 5, which indicates that accuracy of frequency estimates of Si 3 is slightly affected (< 2 dB) by K ′ .
We also calculate |α K+1 /α K | as measurement of the level of spurious sinusoids.Fig. 6 presents the results and shows that the ratios |α K+1 /α K | stay at low level (< 0.2) and are not sensitive to K ′ .Fig. 7 presents MSEs of f 3 estimates versus SNR at different K ′ .Noise variance σ 2 is altered such that SNR 3 varies.The MSEs converges to CRB at high SNR (SNR 3 > 2 dB) when K ′ = K = 3.The results of K ′ = 6 are close to those of K ′ = 3.

Recovering Small Sinusoids
We compare the performance on recovering weak sinusoids of AMP-CTLS with CS methods, e.g.OMP and, CoSaMP, and conventional spectral analysis methods, e.g.ESPRIT and root MUSIC [18].Suppose there CoSaMP iterates 50 times.In both ESPRIT and root MUSIC, the model orders are set as K, and the covariance matrix orders are M/2 according to [32].ESPRIT and root MUSIC output frequency estimates and the corresponding magnitudes are obtained by projection on these frequencies.AMP-CTLS is configured the same as mentioned in Subsection 5.3.
Some intuitive results are presented in Fig. 8.The Si 3 is recovered by the AMP-CTLS algorithm and is masked via other tested algorithms.CoSaMP 100M is also tested, but the results are not displayed because the amplitude estimates are too large (> 1000), which is caused by projection onto the ill-conditional matrix consisting of highly correlated atoms.In OMP 2M , the sinusoids are not exactly at the grid points, so the energies of Si 1 and Si 2 cannot be totally canceled in the beginning two iterations, and the leakage of the energies masks the smallest signal Si 3 .In OMP 100M , all sinusoids are placed at the grid points, and Si 1 and Si 2 are better recovered than in OMP 2M , but energy leakage of dominant sinusoids still exists.In AMP-CTLS, the grid points are adaptively adjusted to match the sinusoids, so the algorithm is less sensitive to grid mismatch and can achieve better performance than OMP and CoSaMP even if the frequency space is sparsely divided.Since there is only one snapshot data, smoothing method [30] is used in ESPRIT and root MUSIC to estimate the covariance matrix.This reduces the performance of frequency estimation.

Range-Velocity Joint Estimate in RSF Radar
In this subsection, we discuss merits of AMP-CTLS in recovering small targets with RSF radar.Suppose there are three targets: two large targets T 1 and T 2 and a small target T 3 .The number of measurements M is 32.The scattering intensities are α 1 = α 2 = 10, α 3 = 1, and the ranges and the velocities are set such that the p, q parameters are p 1 = 10.1/M , p 2 = 10.7/M , p 3 = 20/M , q 1 = 19.4/M, q 2 = 10.2/M and q 3 = 15.2/M .AMP-CTLS is configured as follows: the p, q spaces are both uniformly divided into M grid points; the normalization factors are D p = D q = I/σ ∆ , σ ∆ = 0.025, σ w = 1; and the IJE algorithm iterates fewer than 14 times.In OMP m , both the p and q spaces are uniformly divided into m grid points.Note that all of the targets lie on the grid points in OMP 10M .We focus on the results of recovering the weakest

MSEs of f3
Estimates target T 3 .Change the noise covariance σ 2 ; thus, the signal to noise ratio SNR 3 = α 2 3 /σ 2 varies.Calculate MSEs of p 3 , q 3 parameters.As shown in Fig. 9, the MSEs with AMP-CTLS are lower than those with OMP and converge to the CRB when the SNR 3 is no less than 2 dB.The difference between these MSEs of AMP-CTLS at high SNR and CRB is less than 0.5 dB.

DoA Estimation
In this subsection, AMP-CTLS is compared with the Lasso-based TLS method WSS-TLS [14] on direction of arrival (DoA) estimation.The goal is estimating DoA of plane waves from far-field, narrowband sources with uniform linear array of antennas [14].We focus on the single-snapshot case.Suppose the antenna

Section 6
is dedicated to a brief conclusion.Notations: (•) H denotes conjugate transpose matrix; (•) T transpose matrix; (•) * conjugate matrix; (•) † pseudo-inverse matrix; I L /0 L the L × L identity/zero matrix; • 2 the ℓ 2 norm; {•} denotes a set; | • | the absolute value of a complex number or the cardinality of a set; (•) Λ denotes elements/columns indexed in the set Λ of a vector/matrix; supp(•) is the support set of a vector, that is, the indices of the nonzero elements in the vector; Re(•) the real part of a complex number; ⊗ denotes the right Kronecker product [19]; and E[•] p k ∈ [0 1) and q k ∈ [0 1) are determined by range and radial velocity of the kth target, respectively.Note that in (34) the echo is simultaneously related to the sequence m and the random integer C m .Divide p space into C grid points p c = c/C, c = 0, 1, . . ., C − 1. Divide q space into D grid points p d = d/D, d = 0, 1, . . ., D − 1. Rewrite (33) as:

Figure 1 :
Figure 1: MSEs of the frequency estimates obtained from 1000 independent Monte-Carlo trials via OMP and AMP-CTLS.The distance to the nearest grid point is normalized by 1/M .CRB denotes Cramer-Rao bound.

1
Iteration counter l of IJE Normalized Residual Norm r(l) 2 Normalized Frequency Error | g(l) − f|

Figure 2 :
Figure 2: The normalized errors at different IJE iterations in noiseless case.r(l) 2 denotes the norm of residual error after l-th iteration.f is actual frequency of the sinusoid and ĝ(l) is the estimate of grid point.In the plots, r(l) 2 and |ĝ(l) − f | are normalize with r(0) 2 and |ĝ(0) − f |, respectively.

Figure 3 :
Figure 3: MSEs of the frequency estimates versus distances between the initial grid point and the actual values.The distance is normalized by 1/M .The estimates are obtained from 1000 independent Monte-Carlo trials.'Initial distance (dB)' in the legend denotes the square differences between the frequency grid point and the true ones.

Figure 5 :
Figure 5: MSEs of f 3 estimates versus given sparsity K ′ obtained from 10000 independent Monte-Carlo trials.The results are compared with corresponding CRB.

Figure 7 :
Figure 7: MSEs of f 3 estimates versus SNR at different given sparsity K ′ obtained from 1000 independent Monte-Carlo trials.The results are compared with CRB.
(b) Comparison with ESPRIT, root MUSIC

Figure 8 : 6 Conclusion
Figure 8: Amplitude and frequency estimates of the sinusoids obtained with OMP, CoSaMP, ESPRIT, root MUSIC and AMP-CTLS from 15 independent Monte-Carlo trials.The estimates are compared with the true parameters.Subscripts of OMP and CoSaMP denote numbers of grid points.Note that in OMP 100M , all of the sinusoids exactly lie on the grid points.The results of conventional MP methods and spectral analysis methods are respectively shown in (a) and (b) for better readability.