Augmented Lagrange Based on Modified Covariance Matching Criterion Method for DOA Estimation in Compressed Sensing

A novel direction of arrival (DOA) estimation method in compressed sensing (CS) is presented, in which DOA estimation is considered as the joint sparse recovery from multiple measurement vectors (MMV). The proposed method is obtained by minimizing the modified-based covariance matching criterion, which is acquired by adding penalties according to the regularization method. This minimization problem is shown to be a semidefinite program (SDP) and transformed into a constrained quadratic programming problem for reducing computational complexity which can be solved by the augmented Lagrange method. The proposed method can significantly improve the performance especially in the scenarios with low signal to noise ratio (SNR), small number of snapshots, and closely spaced correlated sources. In addition, the Cramér-Rao bound (CRB) of the proposed method is developed and the performance guarantee is given according to a version of the restricted isometry property (RIP). The effectiveness and satisfactory performance of the proposed method are illustrated by simulation results.


Introduction
Direction of arrival (DOA) estimation of multiple narrowband sources has been an interesting research topic in array signal processing. Its applications span many fields including radar, communication systems, and medical imaging [1,2]. Many effective algorithms are proposed for DOA estimation, which are mainly classified into three categories. The beamforming algorithms such as MVDR [3] and MEM [4] obtain a nonparametric spatial spectrum by optimizing the filter weights. The subspace algorithms such as MUSIC [5] and ESPRIT [6] and their derivatives [7,8] exploit the orthogonality of signal subspace and noise subspace for DOA estimation. The subspace fitting algorithms including Maximum Likelihood (ML) [9] and Weighted Subspace Fitting (WSF) [10] solve a multidimensional nonlinear optimization problem to obtain a precise estimation, but a good initialization is required to ensure global convergence and computational complexity is very high. All these algorithms focus on two important issues: resolution (i.e., the ability to correctly resolve two closely spaced sources) and precision (i.e., the deviation from true DOA) which are considered to be the theoretical bases of evaluating certain algorithms. Many high resolution methods suffer from serious performance degradation in the scenarios with low SNR, small number of snapshots, and closely spaced corrected sources. More recently, many research applications involving compressed sensing (CS) [11], especially DOA estimation [12,13], have become more and more popular in the signal processing. Moreover, the restricted isometry property (RIP) based on the perfect theoretical framework given with modern probability theory by Donoho [11] and Candés et al. [14] provides the performance guarantee in CS. Hence, in this paper, DOA estimation is posed as a joint sparse recovery problem where we recover jointly sparse sources from multiple measurement vectors (MMV) under CS framework.
CS is an emerging area which can break through the limit of Nyquist sampling theorem. On one hand, CS can simultaneously capture and store compressed or sparse source at a rate much lower than that of Nyquist sampling. On the other hand, it can recover the original source using nonadaptive linear projection measurements onto a suitable measurement matrix which satisfies the RIP. In CS, the joint sparse recovery aims to find a common support shared by unknown sparse 2 The Scientific World Journal vectors from MMV, which is obtained by the measurement matrix. Support denotes the indices of the nonzero elements in the unknown sparse vectors. A sparse solution can be obtained as long as the support is determined.
CS theory has been widely applied to DOA estimation according to the source sparsity, which results from the fact that there are much fewer source directions than all potential directions in the spatial domain. The DOA estimation methods in CS are attractive since they have much better estimation performance than conventional estimation methods. In [15], Gürbüz et al. firstly formulate the DOA estimation problem under CS framework in the time domain. Wang et al. [16] propose a new architecture to estimate DOA by exploiting compressive sampling in the spatial domain. Stoica et al. [17] make full use of covariance matching criterion and present a semiparametric/sparse estimation method and its derivative called LIKES [18] for the separable model. Malioutov et al. [19] present the 1 -SVD algorithm for DOA estimation which combines the SVD of the data matrix with a sparse recovery method based on 1 -norm minimization. A new class of subspace-base algorithms such as compressive MUSIC (CS-MUSIC) [20] and subspace-augmented MUSIC (SA-MUSIC) [21] is proposed in recent years.
The RIP and various modified versions of it have been used as a foundation of performance guarantee [21][22][23][24] for the joint sparse recovery. The performance guarantee of MUSIC based on joint sparse recovery is given for identifying the unique support in a favorable case [25]. However, in the unfavorable case where the number of measurement vectors is smaller than the sparsity or the covariance matrix tends to lose rank due to correlated sources, performance guarantee fails. Lee et al. [21] propose a new performance guarantee in terms of a version of the RIP under such unfavorable conditions. The performance guarantees of other methods such as greedy algorithms and convex relaxation have been developed to find the sparse solution in [23,24]. However, the guarantees of these methods cannot be simply extended to the MMV case to obtain a better bound for the sparse solution.
In this paper, we propose a novel augmented Lagrange based on modified covariance matching criterion method called AULMC for DOA estimation in CS. This method utilizes the minimization of the modified covariance matching criterion which is acquired by using regularization method to add penalties in order to obtain the stable sparse parameter estimation, especially in the low sparsity case. This minimization problem is shown to be a semidefinite program (SDP) and transformed into the constrained quadratic programming problem for the sake of reducing computational complexity. The augmented Lagrange function is formed to solve this problem by the use of augmented Lagrange method. AULMC has a number of advantages over the other methods. For example, it provides more precise estimation and higher resolution in the scenarios with low SNR, small number of snapshots, and closely spaced correlated sources, and it does not need priori knowledge about the number of sources or to choose the regularization parameter of the 1 -optimization framework which is very difficult to select in the DOA estimation. In addition, we give a detailed derivation process of the closed-form expression for the Cramér-Rao bound (CRB) of the new method and discuss an explicit condition that guarantees performance of the new method. This performance guarantee is given in terms of weaker version of the RIP which is referred to as weak-1 RIP. Simulation results illustrate the performance of the proposed method.
It is worth noting that covariance matching criterion has been used for DOA estimation [26]. In [26], a sparse iterative covariance-base estimation method, abbreviated as SPICE, is proposed. Our approach is different from SPICE because it utilizes modified covariance matching criterion which can guarantee the stability of solution even if the source sparsity is rather low. In the future work, we will focus on the application of our approach to data-driven design [27][28][29][30]. Now we briefly summarize the contributions of this work as follows.
(i) A modified covariance matching criterion is proposed by adding penalties according to the regularization method.
(ii) The original SDP problem is transformed into the constrained quadratic programming problem. The motivation to transform the problem is that it can reduce computational complexity.
(iii) Augmented Lagrange based on modified covariance matching criterion method is devised to solve the resulting programming problem.
(iv) The CRB and performance guarantee of the new method are given in detail.
The rest of this paper is organized as follows. In Section 2, we formulate the DOA estimation problem. Section 3 presents a modified covariance matching criterion. A novel augmented Lagrange based on modified covariance matching criterion method for DOA estimation is described in detail in Section 4; the performance of which is analyzed in Section 5. The performance of the proposed method is evaluated by simulations in Section 6. Conclusions are provided in Section 7.

Problem Formulation
Consider narrow-band sourceŝ1,̂2, . . . ,̂impinging on a sensor array that consisted of omnidirectional sensors from far field. At time instant , the array received source can be given by where n( ) ∈ C ×1 denotes a noise term, ∈ Ω is the unknown direction, of the th source where Ω denotes the entire spatial range and a( ) is the × 1 steering vector. Although the DOA estimation of the single snapshot, which is a typical single measurement vector (SMV) model, has its value, the number of snapshots is larger than one in the most practical applications. Correspondingly, the multiple snapshots model is a MMV model.
Hence, the multiple snapshots model can be written as the following sparse form: where A = [a(̃1) a(̃2) ⋅ ⋅ ⋅ a(̃)] is the × manifold matrix corresponding to all the potential directions which is also referred to as an overcomplete dictionary in CS.
is a -sparse vector since it has at most nonzero elements in elements, and is defined as sparsity where the operator [⋅] denotes transport. {s( )} =1 are jointly -sparse if they share a common support. Hence, the matrix S = [s( 1 ) s( 2 ) ⋅ ⋅ ⋅ s( )] ∈ C × has no more than nonzero rows in order to be called row -sparse. The MMV problem is that of identifying the row support of the unknown matrix S from the matrix Y ∈ C × that consisted of MMV which is given by with a common measurement matrix Φ of the size × with < where is the number of nonadaptive linear projection measurements, such as random Gaussian matrix or random partial Fourier matrix, and noise matrix N.

Modified Covariance Matching Criterion
In this section, the modified covariance matching criterion is developed in the CS scenario. A conventional covariance matrix of the compressed measurement source with the size × is given by where R = [s( )s * ( )] is a × covariance matrix of the sparse source whose off-diagonal elements denote the source correlation and diagonal elements denote the source powers. Since the powers are one-to-one corresponding to all the potential directions and our focus is on the source angle parameter estimation, R can be reduced to a diagonal matrix R (̃) = diag (̃1̃2 ⋅ ⋅ ⋅̃). According to (5), the measurement matrix can change the covariance matrix of the noise to render the noise colored even if the noise is white. Therefore, this adverse factor must be considered before recovering jointly sparse sources. In addition, the operators (⋅) * and [⋅] denote conjugate transpose and expectation, respectively.
Since ΦR Φ * is a positive definite Hermitian matrix, a prewhitened method is given by the Cholesky decomposition. Let B be the Cholesky factor that satisfies where B ∈ C × is an upper triangular matrix of positive line. Hence, a prewhitened process is implemented by multiplying y( ) by B in order to obtain a pure source whose covariance matrix is written as where C = [c 1 , c 2 , . . . , c ] and I is an identity matrix of the size × . It is important to note that the unknown covariance matrix of the noise is transformed into the known identity matrix after the prewhitened process. Therefore, the prewhitened method improves the robustness to the noise. Then, the covariance matrix of the compressed measurement source denoising is realized by The parameter̃can be estimated by a class of the covariance matching estimation techniques (COMET) based on covariance matching criterion [31]. This parameter estimation method is well known to be a large-sample approximation of ML method and provide a more attractive solution than ML estimator.
The principle of COMET is that of using the right data to minimize its data model by the weighted least-squares (WLS) method. However, the lower the source sparsity is, the more likely it is to be ill-posed for the covariance matrix estimation error meaning that the optimal solution obtained directly by minimizing the conventional covariance matching criterion is instable. Hence, we employ regularization method to add penalties in order to sufficiently exploit prior knowledge to reduce the solution space for determining the stable optimal solution. The modified covariance matching criterion is proposed as the following form: where parameter ≥ 0 controls the solution smoothness (guarantee precision), parameter ≥ 0 controls the solution scale (guarantee sparsity), the inverses of the matrix Γ, the sample covariance matrixR, and R(̃) are assumed to exist, and the matrix P −1 is the inverse of the covariance matrix of the residuals,̃= Γ vec(R − R(̃)). SinceR is equal to R(̃) as → ∞, it follows from [32] that vec(R − R(̃)) satisfies the asymptotic normal distribution with mean zero and variance −1 R (̃) ⊗R. In addition, the operators Tr[⋅], vec(⋅), and ⊗ denote matrix trace, column stacking operation, 4 The Scientific World Journal and Kronecker product, respectively. Then, the matrix P can be given by We consider the normalized data model in (9) and choose the regularization parameters = = based on perceptual criterion [33]. By substituting (10) into (9), we have By the properties of vec, ⊗, and Tr, the data model can be further simplified to where ‖ ⋅ ‖ denotes the Frobenius norm for matrices or the Euclidean norm for vectors. The data model in (12) is considered to be the modified covariance matching criterion. It can be seen from (12) that a positive definite Hermitian matrix is added to the covariance matrix estimation error by applying penalties according to regularization method in order to guarantee the stability of solution.

DOA Estimation
In this section, we will utilize the minimization of the modified covariance matching criterion to estimate parameter̃. Let̃= be the optimal solution of̃in the structure model of (8). Our goal is to utilize the modified covariance matching criterion for an estimate that is asymptotically equal tõ. By the properties of the trace and the Hermitian matrix, a derivation process is shown as follows, where we omit the dependence oñfor notational convenience: Due to we can deduce that the minimization of is equal to the minimization of ℎ: Then, we will demonstrate that the minimization of ℎ in (16) with respect to {̃} =1 is an SDP. To prove this fact, One assumesR The ℎ in (16) can be rewritten as By the Schur complement, let { } =1 be auxiliary variables satisfying The equivalent form of this minimization problem is expressed as It is easy to see that (20) is an SDP [34]. Many software packages can solve an SDP, but solving (16) as an SDP is not a good choice because SDP solvers have generally rather high computational complexity for the values of , , and in the DOA estimation. To solve it effectively, we transform it into the constrained quadratic programming problem for reducing computational complexity, as described next.
Since a consistent estimation of̃is given by (15), we can reformulate the minimization of ℎ in (16) as the following constrained minimization by the Schur inequality of the trace: The Scientific World Journal is a -dimensional vector with = −1/2 , = 1, 2, . . . , . By substituting (8) into (21), the objective function in (21) can be rewritten as Based on the following equation: where ⊙ denotes the Schur-Hadamard product, M and N are both × matrices, T(d) = diag( 1 2 ⋅ ⋅ ⋅ ), and d = [ 1 2 ⋅ ⋅ ⋅ ] , the objective function (22) can be written as where Hence, the minimization problem in (21) is transformed into the following form: which is a typical constrained quadratic programming problem. It is well known that an important class of methods for solving the constrained quadratic programming problem is to form the auxiliary function. To solve (25), we form the following augmented Lagrange function with respect to (25): where is the asymptotical solution of the Lagrange multiplier of (25) and is a penalty factor. By setting the gradient and the Hessian matrix of (26) with respect to z to zero, we have where (z, ) = (z) + * Π(z) and q(z) = ∇Π * (z) = [∇Π 1 (z) ∇Π 2 (z) ⋅ ⋅ ⋅ ∇Π (z)]. One assumes that b (z, ) = ∇ 2 zz (z, ) + ∑ =1 Π (z) ∇ 2 Π (z) .
Applying the Newton method, we obtain For notational convenience, we assume that q = q(z ), ∇ = ∇ (z ), Π = Π(z ), b = b(z , ), d = z +1 − z and we get the following equation by (29) Assuming that the inverse of b exists, by left-multiplying (30) by (1/ )q * b −1 , we have It follows from (31) that where satisfies By substituting (32) into (30), we have Both the multiplier factor and the penalty factor are determined with difficulty for utilizing the augmented Lagrange function to solve the constrained quadratic programming problem. Hence, an updated sequence for the multiplier factor is given in terms of Proposition 1.

Proposition 1.
One assumes that (z) is the optimal solution of the problem min ∈C ‖q(z) + ∇ (z)‖ 2 . Then, the following equation holds: Proof. See Appendix A.
It can be deduced from Proposition 1 that + is referred to as the next iteration of . We apply a heuristic update sequence for the penalty factor to achieve a stable solution. If the th iterative solution z is closer to the feasible region than the previous solution z −1 , the penalty factor is decreased. Inversely, we increase the penalty factor when z is not closer to the feasible region.
The specific steps of solving by the augmented Lagrange method are given as follows.
(2) We use Armijo search method to find the maximum of satisfying (3) Set z +1 = z + d and update the multiplier factor and penalty factor, respectively: (4) Set = + 1 and return to step (1).

Cramér-Rao Bound.
In this subsection, the closed-form expression for the CRB of the proposed method with complex white Gaussian noise after the prewhitened process is illustrated. The bound of the noise variance estimation can be computed separately as CRB = 1/ (see [35]). The remaining parameters consist of an unknown vector = [̃s ] . It is not an easy task to get the CRB of the unknown parameters. However, fortunately, [36,37] have provided a critical inspiration for the derivation in this paper. First, the likelihood function is given by ( y s ( ) ,̃) = 1 Then, the partial derivatives of (39) with respect tõ, s ( ) = Re{s( )} and s ( ) = Im{s( )} are given by where S = diag(s( )) and U = [ a(̃1)/̃1 a(̃2)/ 2 ⋅ ⋅ ⋅ a(̃)/̃]. Following [35,37], we can obtain the Fisher information matrix (FIM) as follows: where It is well known that It can be deduced from (41) and (43) that where H = BΦ is a × matrix and (⋅) + denotes pseudoinverse. Note that the CRB in CS is affected not only by the conventional factors, for example, SNR, array structure and signal relative location, but also by the measurement matrix.

Performance Guarantee.
In CS, the RIP has been deeply studied for the joint sparse recovery by minimizing the 1 norm. We say that the matrix C ∈ C × obeys the RIP of the order if there exists a constant ∈ (0, 1) satisfying Therefore, all submatrices of C with L columns are uniformly well conditioned. The restricted isometry constant (RIC) of the order L, described as (C), is the smallest that satisfies (45) and (C) satisfies where is a -dimensional subsupport of = [1, 2, . . . , ] and C denotes the submatrix of C with columns indexed by . Note that the condition satisfying the RIP is so demanding that its applications are limited. Therefore, we should make use of a new version of the RIP, which is called the weak-1 RIP [38], to control the size of the recovery error. The weak-1 RIP is given in the following form: for all k supported on , where the cardinality of the set is +1. If the matrix C satisfies the weak-1 RIP, it can be deduced that 0 ≤ ≤ +1 (C ), where (C ) denotes the th largest eigenvalue of C . The corresponding weak-1 RIC is given by In this paper, when the estimation quality is imperfect, especially in the unfavorable case, the support is no longer identified by the algebraic property of R. Hence, a new performance guarantee is given in terms of the weak-1 RIP in the following proposition.
Proof. See Appendix B.
It follows from Proposition 2 that the performance guarantee of the proposed method requires a mild condition in the unfavorable case. However, in the favorable case, the performance guarantee only requires a much milder condition, +1 (R; 0 ) > 0, which is an algebraic condition.

Simulation Results
In this section, the performance of the proposed method is illustrated by several simulation results and compared with that of CS-MUISC, SPICE, and CRB under the condition that the number of sources is unknown. We consider the spatial signal impinging on the uniform linear array (ULA) with interspacing /2 where denotes the wavelength of source. In the ULA case, the steering vector corresponding to the DOA equal to is given by where the number of the array elements is set to be = 8.
In the simulation, the average root mean square error (RMSE) of the DOA estimation with 50 Monte Carlo runs is defined as the significant performance index: where is the estimate of in the th run. The resolution of the grid is closely related to the precision of the DOA estimation. A coarse grid can lead to poor precision, but a too fine grid increases computational complexity. Therefore, an adaptive grid refinement method is used to 8 The Scientific World Journal balance the tradeoff between precision and computational complexity. In the simulation, we make a coarse grid with 1 ∘ step in the range of −90 ∘ to 90 ∘ and perform a local fine grid in the vicinity of locations obtained by using the coarse grid.
In the first simulation, we display the superimposed spatial spectra of three algorithms in 10   as follows: the spatial spectrum obtained by CS-MUSIC suffers from severe interference at the true directions, especially at two most closely spaced correlated sources whose bias is clearly seen from the insert. Two most closely spaced correlated sources can be resolved by SPICE (note that peaks in the spatial spectrum are much larger than 2, but they are cut off at 2 to use the same scale as the other figures in Figure 1), but SPICE can yield false peaks and slight bias in the vicinity of the correlated sources and uncorrelated sources, respectively. The proposed method AULMC yields a nearly ideal spatial spectrum and provides a precise estimation for all the sources. In summary, AULMC outperforms CS-MUSIC and SPICE in terms of the spatial spectrum in the scenario with low SNR, small number of snapshots, and closely spaced correlated sources.
We analyze the RMSE of three algorithms under different conditions in the second simulation. The source model is the same as the first simulation. Figure 2 shows the RMSE as a function of SNR of all the algorithms and CRB in 50 Monte Carlo runs for the fixed number of snapshots 50, whereas the RMSE versus number of snapshots is shown in Figure 3   and 3, we can draw the conclusions that the RMSE of AULMC is smaller than those of other two algorithms and AULMC has the more significant performance advantages than the other two algorithms, especially in the scenarios with low SNR or small number of snapshots. One possible explanation is that AULMC can give the stable estimation in every Monte Carlo run. It can be also seen that the RMSE is close to the CRB with the increase of SNR and the number of snapshots.
In Figure 4, we display the relation between the RMSE and angel separation of correlated sources which can illustrate the resolving ability. Let two correlated sources at angles 20 ∘ and 20 ∘ +Δ , where the step of Δ is 1 ∘ , be impinged on the ULA. The SNR is 13 dB and the number of snapshots is 100. It can be seen from Figure 4 that when angle separation is 2 ∘ , AULMC fails; however, AULMC can still provide a precise estimation as long as the angle separation is no less than 3 ∘ and has higher resolution than the other two algorithms.
Finally, the computation time of different algorithms versus number of snapshots is shown in Table 1 for comparing the efficiency of these algorithms and SDP solver. Two correlated sources impinge on the ULA at 20 ∘ and 25 ∘ . The SNR is fixed at 13 dB. The computation time is obtained by the MATLAB 7.8 (R2009a) on a 2.8 GHz 4 GB PC. For AULMC, the computation time is mainly spent on the iterations of augmented Lagrange.
It can be seen from Table 1 that the computation time of SDP solver is the longest, and although the computation time of AULMC is longer than that of other two algorithms, it is comparable. Moreover, it is worth noting that the performance of AULMC is much better than that of CS-MUSIC or SPICE.

Conclusion
A novel augmented Lagrange based on modified covariance matching criterion method for DOA estimation is proposed in CS. It is proved that the problem of minimizing the modified covariance matching criterion is an SDP, which can be transformed into the constrained quadratic programming problem solved by the augmented Lagrange method. A detailed derivation for the CRB and a theoretical performance guarantee for identifying the support are provided. Simulation results show that AULMC outperforms CS-MUSIC and SPICE in terms of the spatial spectrum and has more precise estimation as well as higher resolution, especially in the scenarios with low SNR, small number of snapshots, and closely spaced correlated sources.