# Superresolution 2D DOA Estimation for a Rectangular Array via Reweighted Decoupled Atomic Norm Minimization

Authors: Ming-Ming Liu; Chun-Xi Dong; Yang-Yang Dong; Guo-Qing Zhao
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This paper proposes a superresolution two-dimensional (2D) direction of arrival (DOA) estimation algorithm for a rectangular array based on the optimization of the atomicl0 norm and a series of relaxation formulations. The atomic l0 norm of the array response describes the minimum number of sources, which is derived from the atomic norm minimization (ANM) problem. However, the resolution is restricted and high computational complexity is incurred by using ANM for 2D angle estimation. Although an improved algorithm named decoupled atomic norm minimization (DAM) has a reduced computational burden, the resolution is still relatively low in terms of angle estimation. To overcome these limitations, we propose the direct minimization of the atomic l0 norm, which is demonstrated to be equivalent to a decoupled rank optimization problem in the positive semidefinite (PSD) form. Our goal is to solve this rank minimization problem and recover two decoupled Toeplitz matrices in which the azimuth-elevation angles of interest are encoded. Since rank minimization is an NP-hard problem, a novel sparse surrogate function is further proposed to effectively approximate the two decoupled rank functions. Then, the new optimization problem obtained through the above relaxation can be implemented via the majorization-minimization (MM) method. The proposed algorithm offers greatly improved resolution while maintaining the same computational complexity as the DAM algorithm. Moreover, it is possible to use a single snapshot for angle estimation without prior information on the number of sources, and the algorithm is robust to noise due to its iterative nature. In addition, the proposed surrogate function can achieve local convergence faster than existing functions.


---

## Body

## 1. Introduction

The problem of two-dimensional (2D) direction of arrival (DOA) estimation is encountered in various signal processing applications, such as 2D angle estimation using 2D arrays in multiple-input and multiple-output (MIMO) communication systems. Although numerous algorithms for this purpose have emerged in the past three decades [1–5], there is still an urgent need for fast and superresolution algorithms that can utilize bits of snapshots, especially for massive systems.Conventional superresolution algorithms such as MUSIC and ESPRIT [6–9] are based on subspace techniques. Considering that these methods either require a heavy computational load for the 2D spectrum search or are suitable only for arrays that exhibit geometric invariance based on a number of samples that is larger than the number of sensors, we instead focus on algorithms that rely on the popular theory of compressed sensing (CS), which makes it possible to use a smaller number of snapshots, or even a single snapshot, to estimate the DOA, even when the sources are correlated. However, another challenge arises due to grid mismatch [10, 11].Recently, the mathematical concept of continuous compressed sensing (CCS) was introduced by Candès and Fernandes-Granda [12]. CCS no longer follows traditional grid-based models, which require linear independence among the grids, but rather relies on a superresolution concept based on the total variation norm; thus, it fundamentally avoids the problem of grid mismatch and further facilitates the application of CS in parameter estimation. Subsequently, Chandrasekaran et al. [13] and Tang et al. [14] developed an atomic norm metric and presented an atomic norm minimization (ANM) problem in the positive semidefinite (PSD) form utilizing the Vandermonde structure of the atoms. The theory was further developed for the estimation of one-dimensional (1D) parameters by Yang et al. [15], who also proposed the reweighted atomic norm minimization (RAM) algorithm to enhance the resolution [16]. In addition, the metric was extended for higher-dimensional parameter estimation by Xu et al. [17] and the authors [18, 19]. However, the computational complexity is enormous even for 2D applications. To combat this issue, Tian et al. [20] presented a decoupled atomic norm minimization (DAM) algorithm in which the two-level Toeplitz matrix is divided into two one-level matrices, thereby improving the calculation efficiency. However, the resolution is limited because the minimum angular interval is inversely proportional to the number of rows or columns in the array [12, 19, 20].To address the above issues, we propose a low-complexity superresolution optimization algorithm for 2D DOA estimation with a rectangular array. In essence, the proposed algorithm adopts an extended atomicl0 norm that avoids the constraint on the angular intervals, and we argue that the problem can be translated into a decoupled rank minimization problem in the PSD form, whereafter since the rank minimization problem is difficult to solve directly, a novel sparse surrogate for it is presented, inspired by the works [16, 21]. Then, the new concave objective function can be iteratively optimized by the well-known majorization-minimization (MM) method. Our algorithm eventually takes the form of a weighted DAM problem in each iteration; thus, we name it reweighted decoupled atomic norm minimization (RDAM). A series of proofs are presented in this paper to illustrate the rationality of the proposed method, and corresponding simulations are also provided to show its validity. Our proposed algorithm enhances the resolution to a great extent without increasing the computational burden, and it requires neither prior information on the number of sources nor the noise level which can be obtained in the iterative process. In addition, our method is more robust to noise than the DAM algorithm due to its iterative nature.


## 2. Signal Model

Suppose thatK far-field narrowband sources {sk(t)}k=1K are impinging on a uniform rectangular array (URA) composed of N×M sensors, where N and M are the numbers of elements in the x and z directions, respectively. For the sake of subsequent discussion, we adopt a decoupled observation model for the array output [22]. As shown in Figure 1, we redefine the azimuth θ of the received signal as the angle between the signal and the yoz plane instead of using the traditional definition; the definition of the elevation φ remains unchanged. Correspondingly, the single-snapshot array observation model without noise is(1)Xo=∑k=1KskaNθkaMHφkwhere Xo is the array response, with dimensions of N×M, and sk=skejϑk is the k-th received source, with sk denoting the amplitude of the source and ϑk denoting the phase, which takes values in the range 0,2π (the frequency part of each source is ignored for simplicity). Moreover, aN(θk)=[ξ1k,ξ2k,…,ξNk]T and aM(φk)=[η1k,η2k,…,ηMk]T are the 1D steering vectors in the x and z directions, respectively. Suppose that the array spacing is set to half the wavelength; then, we have ξnk=exp⁡(-j2π·0.5n-1sin⁡θk) and ηmk=exp⁡(j2π·0.5m-1sin⁡φk). Note that ηmk is the conjugate of the regular vector.Figure 1
Signal model.Our goal in this paper is to estimateθk and φk from the observed data Xo, and we assume that the following conditions hold:(1)aN(θk) and aM(φk) have the Vandermonde structure.(2)0.5sin⁡θk∈[-1/2,1/2] and 0.5sin⁡φk∈[-1/2,1/2].(3) The phases{ϑk}k=1K are either i.i.d. and uniformly drawn from the complex unit circle or independent random variables with mean Eϑk=0.The first two conditions, which obviously hold for array models, are the basic requirements for the algorithm subsequently proposed in this paper, and the last one guarantees that the solution of our proposed algorithm is exact with a probability of at least1-δ, where δ is a very small value [14, 19, 23].


## 3. Proposed Reweighted Decoupled Atomic Norm Minimization (RDAM) Algorithm

### 3.1. Decoupled Rank Minimization Problem

The atomicl0 norm was proposed in [14] as a measure of the sparsest number of spikes with unknown positions that could compose the observed signal in the frequency domain. In this paper, we employ this norm to express the number of sources with unknown directions composing the array output in the spatial domain and implement 2D DOA estimation by optimizing the atomic representation.More precisely, we merge the signal phase information into the steering vectors and redefine them as follows:(2)aNθk,12ϑkn=exp⁡-j2π·0.5n-1sin⁡θk+12ϑk,n=1,…,N(3)aMφk,12ϑkm=exp⁡j2π·0.5m-1sin⁡φk12ϑk,m=1,…,MRegarding aN(θk,1/2ϑk)aMH(φk,1/2ϑk) as an atom, we define the atom set B as(4)B=aNθ,12ϑaMHφ,12ϑ,θ∈0,2π,φ∈0,2π,ϑ∈0,2π=Aθ,φ,ϑ,θ∈0,2π,φ∈0,2π,ϑ∈0,2πThen, the atomic l0 norm induced by B is given by(5)XB,0=infK:X=∑k=1KskAθk,φk,ϑk,Aθ,φ,ϑ∈BThis is an optimization problem with the objective of finding the fewest possible number of sources given the observed data Xo, which can be represented as a positive combination of the atoms {A(θk,φk,ϑk)}. Inspired by the idea proposed in [14] and previous research on rank minimization [24–26], we will demonstrate that the above norm can be represented as the sum of the ranks of two 1D Toeplitz matrices satisfying a PSD constraint and propose a decoupled rank minimization problem as depicted by the following proposition.Proposition 1.
LetXo∈CN×M be an array observed data that has the following atomic form and satisfies the conditions given in the previous section:(6)Xo=∑k=1KskAθk,φk,ϑkThe corresponding atomic l0 norm minimization problem can be expressed as(7)minXXB,0s.t.X=XoThen, (7) is equivalent to the optimum value of the following decoupled rank minimization problem:(8)minuθ,uφ,X12rankTuθ+12rankTuφs.t.TuφXHXTuθ≥0,X=Xowhere T(uθ) and T(uφ) denote two 1D Toeplitz matrices and uθ and uφ are their first rows, respectively.Proof.
The trivial case ofX=0 is unnecessary to discuss. For X≠0, let r1★ and r2★ be the optimal values of rank(T(uθ)) and rank(T(uφ)), respectively; then, the corresponding optimal value of (8) can be expressed as 1/2r1★+1/2r2★. Furthermore, let the optimal XB,0 be denoted by K★.
On one side, we will show that1/2r1★+1/2r2★≤K★ given problem (7). Suppose that XB,0=K★<min⁡N,M can be realized in the arbitrary atomic decomposed form X=∑k=1K★skaN(θk,1/2ϑk)aMH(φk,1/2ϑk), where {sk}k=1K is a set of positive real numbers. Let uθ=∑kskaN(θk,0) and uφ=∑kskaM(φk,0); then, the homologous Toeplitz matrices can be represented as(9)Tuθ=∑kskaNθk,0aNHθk,0=∑kskaNθk,12ϑkaNHθk,12ϑk≥0and(10)Tuφ=∑kskaMφk,0aMHφk,0=∑kskaMφk,12ϑkaMHφk,12ϑk≥0Thus, we obtain the following PSD matrix:(11)TuφXHXTuθ=∑kK★skaMφk,12ϑkaNφk,12ϑkaMφk,12ϑkaNφk,12ϑkH≥0;i.e., the relational expressions r1★≤K★ and r2★≤K★ hold. Hence, 1/2r1★+1/2r2★≤K★.
On the other side, given problem (8), we will try to verify that K★≤1/2r1★+1/2r2★ in the case of r1★<N and r2★<M (the full-rank case is trivial). Given the PSD condition, and letting uθ★ and uφ★ compose the optimal solution to (8), the following can be derived in accordance with the Schur Complement Lemma [27, 28]:(12)Tuφ★XHXTuθ★≥0⟹Tuφ★≥0IM-Tuφ★Tuφ★+XH=0(13)Tuθ★XXHTuφ★≥0⟹Tuθ★≥0XHIN-Tuθ★Tuθ★+=0In addition, note that T(uθ★) and T(uφ★) have Vandermonde decompositions of the forms T(uθ★)=VθDθVθH and T(uφ★)=VφDφVφH [14], where Vθ=[aN(θ1,0),…,aN(θr1★,0)], Vφ=[aM(φ1,0),…,aM(φr2★,0)], and both Dθ (with dimensions of r1★×r1★) and Dφ (with dimensions of r2★×r2★) are diagonal matrices with positive real elements. Therefore, (12) implies that X lies in the range of Vθ, which means that X can be composed of no more than r1★ atoms; i.e., K★≤r1★. Similarly, (13) leads to the conclusion that K★≤r2★. Thus, K★≤1/2r1★+1/2r2★.
Based on the above proposition, we can straightforwardly optimize the decoupled rank minimization problem (8) under the PSD constraint to implement the optimization of (7). It is evident from the above proof that θk and φk are encoded in the Toeplitz matrices T(uθ★) and T(uφ★), respectively. Once the optimal solution (uθ★,uφ★,X★) is obtained, θk and φk can be determined from the Vandermonde decompositions of these two matrices; there are many off-the-shelf computational methods such as the dual root-finding method [12, 14, 15] and the matrix pencil method [29, 30]. In addition, supposing θk and φk are acquired, the simple and straightforward method of angle pairing provided in [20] can be applied. Therefore, our task is simply to optimize problem (8) for (uθ★,uφ★,X★).


### 3.2. A Novel Sparse Surrogate for Rank Minimization

The essential objective of the decoupled rank minimization problem (8) discussed above is to recover two low-rank matrices within the PSD feasible set. Unfortunately, due to the difficulty of solving this problem (which is similar to that of solving the original optimization problem), there is an urgent need to seek new metrics to express the objective function. The main idea is to find one simpler surrogate function that not only can effectively simulate the ranks of the two Toeplitz matrices but also can simultaneously encourage a low-rank solution. In most papers that discuss the rank minimization problem [16, 27, 31, 32], the log-det function f1(T(u)), which is expressed as follows, is generally used:(14)f1Tu=ln⁡Tu+ε1I=ln⁡1ε1Tu+I+ln⁡ε1However, in this paper, we consider a novel concave surrogate function for T(u)∈CQ×Q, given by (15)f2Tu=-ε2trTu+ε2I-1where I represents the identity matrix. The regularization parameter ε2>0 enables the function f2(T(u)) to approximate rank(T(u)) as it varies and prevents the function from taking a value of -∞ when T(u) is rank-deficient. Moreover, note that the parameter ε1 does the same job with f1(T(u)).From the properties of matrix eigenvalues and trace, it follows that(16)eigTu+ε2I-1=1λi+ε2,i=1,…,Qand(17)f2Tu=-∑ε2λi+ε2=-∑1λi/ε2+1,i=1,…,Qwhere {λi}i=1Q is the set of eigenvalues of T(u). Equation (17) implies that the essence of function f2(T(u)) is to operate on {λi}i=1Q; thus, to better understand the rank-approximation property of f2(T(u)) on a PSD definition domain, we denote a new subfunction by p2(λi)=-(λi/ε2+1)-1+1 for any eigenvalue λi≥0, which is an inverse function obtained through a certain translation and scaling process and has the features of monotonicity and concavity. Note that the constant term 1 is added merely for consistency with the l0 norm; it is equivalent to adding an identity matrix to f2(T(u)) and does not influence the outcome. Obviously, we have the identity p2(0)=0; at λi>0 and near the zero point, p2(λi) has an increasingly steep slope of 1/ε2 as ε2→0, which equivalently imposes a large penalty on near-zero λi and ensures that the values of p2(λi) at these λi become much smaller. In addition, p2(λi) tends towards a constant value of 1 for a sufficiently large λi. The above trend illustrates that as ε2→0, p2(λi) favours the separation of λi at different levels, thus becoming a relaxation of the l0 norm λil0, which has an infinite slope near the zero point and a constant value of 1 except at the zero point. To provide an intuitive illustration, we plot p2(λi) along different ε2’s at λi≥0 in Figure 2. Meanwhile, the curves of the l0 norm, the l1 norm, and the log-det subfunction represented by p1(λi)=ln⁡(λi/ε1+1) are provided for comparison. Suitable proportionality constants are set to ensure that p2(λi) and p1(λi) are zero at λi=0 and have a constant value of 1 at λi=1 without modifying their properties. As shown in Figure 2, both p2 and p1 approach the l0 norm as the regularization parameters decrease; the smaller ε2 is, the steeper the slope is, which offers a basis for the selection of ε2 during both initialization and iteration, as discussed in Section 3.3. Moreover, the function p2 is more likely to perform better since it has a much faster speed of approximation; i.e., p2 is deemed to be more sparsity-enhancing.Figure 2
An illustration of the sparsity-enhancing properties of p1 and the proposed p2, compared with thel0 norm and l1 norm.The approximate relationship betweenp2(λi) and λil0 provides an explanation for why rank(T(u)) can be relaxed by f2(T(u)). In particular, since each larger λi corresponds to a value of p2(λi) that tends towards a constant value of 1 as ε2→0 and the addition operation does not affect the monotonicity and concavity of a function, the sum of all p2(λi) becomes a relaxation of the number of larger eigenvalues, which is considered to be the rank of the Hermitian matrix Tu or the number of sources; i.e., the sparse surrogate f2(T(u)) can be regarded as a relaxation of rank(T(u)). Furthermore, the proposed function f2(T(u)) leads to better performance than f1(T(u)) in encouraging low-rank solutions, as seen by comparing the curves for p2 with the curves for p1; stronger interpretations will be presented in later subsections.


### 3.3. RDAM via the Majorization-Minimization (MM) Algorithm

Using the surrogate functionf2(T(u)), the decoupled rank minimization problem (8) can be relaxed as follows:(18)minuθ,uφ,X-12εθtrTuθ+εθI-1-12εφtrTuφ+εφI-1s.t.TuφXHXTuθ≥0,X=XoHere, the parameters εθ and εφ are deemed to be mutually independent. In this subsection, we aim to minimize the objective function in (18), denoted by F(T(uθ),T(uφ)), on the PSD cone PSD(uθ,uφ). Considering that F(T(uθ),T(uφ)) is a sum of two concave functions, the MM algorithm [33, 34] can be adopted to obtain a local optimum. This algorithm follows an iterative procedure consisting of two steps, i.e., majorization and minimization, in each iteration.Given an initial point{T(uθ,0),T(uφ,0)}∈PSD(uθ,uφ), a series of feasible points {T(uθ,j),T(uφ,j)}j∈N are generated. At each feasible point, the majorization step requires the construction of a tractable function majorizing F(T(uθ),T(uφ)) while ensuring that its value at the current point {T(uθ,j),T(uφ,j)} is equal to F(T(uθ,j),T(uφ,j)). Here, we use the first-order Taylor expansion for this purpose. Since both terms of F(T(uθ),T(uφ)) possess the same structure, we consider only the first term as an example, for which we have the following first-order approximate expression:(19)εθtrTuθ,j+εθI-1+trεθTuθ,j+εθI-2Tuθ-Tuθ,j=εθtr(Tuθ,j+εθI-2Tuθ+cθ,jwhere cθ,j is a constant. The derivation is presented in Appendix B. Thus, F(T(uθ),T(uφ)) has the following upper bound:(20)FTuθ,Tuφ≤εθ2tr(Tuθ,j+εθI-2Tuθ+εθ2tr(Tuφ,j+εφI-2Tuφ+Cj=ΔGTuθ,Tuφ:Tuθ,j,Tuφ,jwhere Cj=cθ,j+cφ,j is constant. It is evident that GTuθ,j,Tuφ,j:Tuθ,j,Tuφ,j=FTuθ,j,Tuφ,j. In addition, many other construction methods [33] (e.g., second-order Taylor expansion, inequality, and the Schur complement) could be used; these possibilities are left for future work.Then, in the minimization step, the following minimization problem is solved, ignoring the constantCj, to update T(uθ,j+1) and T(uφ,j+1):(21)minuθ,j+1,uφ,j+1,Xj+112tr(εθTuθ,j+εθI-2Tuθ+12tr(εθTuφ,j+εφI-2Tuφs.t.TuφXHXTuθ≥0,X=Xjwhere the optimization problem is a semidefinite programming (SDP) that can be solved in polynomial time. Notice that this problem is equivalent to a weighted DAM problem in the PSD form [20], with weighting factors of Wθ,j=εθ(T(uθ,j)+εθI)-2 and Wφ,j=εθ(T(uφ,j)+εφI)-2, and that the proposed algorithm is iterative in nature; therefore, our algorithm is called RDAM. In addition, we have the following inequality: (22)FTuθ,j+1,Tuφ,j+1≤GTuθ,j+1,Tuφ,j+1:Tuθ,j,Tuφ,j+Cj≤GTuθ,j,Tuφ,j:Tuθ,j,Tuφ,j+Cj=FTuθ,j,Tuφ,j;i.e., {F(T(uθ,j),T(uφ,j))}j∈N is a nonincreasing sequence, implying that the iterative process makes the objective function F(T(uθ),T(uφ)) repeatedly smaller, and it can be expected to converge to a local minimum F(T(uθ★),T(uφ★)), as proven in detail in [21, 33]. Furthermore, since the algorithm is locally convergent, suitable initial values are indispensable for guaranteeing convergence to the desired local optimum. One immediate possibility is to use the DAM solution for initialization. In particular, let Wθ,0 and Wφ,0 be identity matrices; then, we can obtain the initial values T(uθ,1) and T(uθ,1) by solving a DAM problem. Next, let us emphasize the selection of εθ,1 and εφ,1, both of which influence the outcome. As argued in Section 3.2, F(T(uθ),T(uφ)) is the sum of two subfunctions, denoted by p2(λθ,i)=-(λθ,i/εθ+1)-1 and p2(λφ,i)=-(λφ,i/εφ+1)-1, for θ and φ, respectively. For consistency with the curves in Figure 2, we set εθ=aλθ,max and εφ=bλφ,max. Then, both the factors a and b can be viewed as equivalent to the parameter ε2 in this figure, and we have λθ,i/λθ,max∈[0,1] and λφ,i/λφ,max∈[0,1]. Now, let us consider the choice of a as an example. Let larger eigenvalues, corresponding to the signals, be denoted by λθ,si, while smaller eigenvalues, corresponding to the noise, are denoted by λθ,ni; thus, we should be able to ensure a relatively large difference in the values of p2(λθ,i) between the points λθ,si/λθ,max and λθ,ni/λθ,max by tuning the initial factor a. For example, for points λθ,si/λθ,max and λθ,ni/λθ,max that are close to each other, we can select a small a value, such as 0.1 or smaller; then, there will be a sharp slope such that the subfunction values will be distinctly separated. In contrast, for the relatively distant points, a large a value, such as 1, will be sufficient to achieve good estimation performance. The same analysis applies for the initialization of the factor b.Finally, it is instructive to explore the potential implications of the RDAM algorithm in combination with the atomic norm [13, 14]. Notice that the concept of the 1D weighted atom set proposed in [16] by Yang et al. can be extended to a new 2D weighted atom set Bw as follows:(23)Bw=wθθaNθ,12ϑwφ∗φaMHφ,12ϑ,θ,φ∈0,2π×0,2π=Awθ,φ,ϑ,θ,φ∈0,2π×0,2πwhere {wθ(θ),wφ∗(φ)} denotes a pair of nonnegative weights on the 2D space angles. Accordingly, the weighted atom norm is given by(24)XBw=inf⁡t>0:X∈tconvAw=inf⁡∑kskw:X=∑kskwAwθk,φk,ϑk,Awθ,φ,ϑ∈Bw=inf⁡∑kskw:X=∑kskwwθθkwφ∗φkaNθk,12ϑkaMHφk,12ϑk=inf⁡∑kskwθθkwφ∗φk:X=∑kskAθk,φk,ϑk,Aθk,φk,ϑk∈BThen, a 2D weighted ANM problem in decoupled PSD form can be obtained based on the following proposition.Proposition 2.
LetXBw,j+1 denote the (j+1)th 2D weighted atomic norm. The following equivalence relation holds given the data Xj∈CN×M recovered from the jth iteration.(25)min⁡XBw,j+1=minuθ,j+1,uφ,j+112trWθ,jTuθ+12trWφ,jTuφs.t.TuφXHXTuθ≥0,X=Xjwhere(26)wθ,jθ=1aNHθWθ,jaNθandwφ,jφ=1aMHφWφ,jaMφ.This implies that iteratively minimizing the weighted atomic normXBw induced by the 2D weighted atom set Bw serves the same purpose as the proposed RDAM algorithm; i.e., RDAM can also be called decoupled reweighted atom norm minimization (DRAM). More significantly, our proposed trace function f2(T(u)) gives rise to quadratic weighting factors, which greatly accelerates the speed of local convergence and further enhances the preference for low-rank solutions in (18) compared with the first-order factors generated by the log-det function f1(T(u)) [16].


### 3.4. Another Interpretation from the Perspective of Conventional Direction of Arrival (DOA) Estimation Methods

The proposed RDAM algorithm can also be interpreted from the perspective of traditional DOA estimation methods. Suppose that the signals are independent of each other and that all steering vectors have been normalized, without exception. Then, we can redefine two decoupled covariance matrices for the noiseless observation data as follows:(27)R~θ≜XXH=∑k=1KaNθkaMHφksktsk∗taMφkaNHθk=∑k=1Kskt2aNθk,12ϑkaNHθk,12ϑk(28)R~φ≜XHX=∑k=1KaMφkaNHθksktsk∗taNθkaMHφk=∑k=1Kskt2aMφk,12ϑkaMHφk,12ϑkFrom the poof of Proposition 1, we know that the Toeplitz matrices T(uθ) and T(uφ) are exactly equivalent to the two decoupled covariance matrices above except that the coefficients include a squared difference. The whole RDAM process, therefore, can be viewed as the process of recovering R~θ and R~φ with the regulators εθ and εφ. Furthermore, the two atomic weight values in (24) are related to the Capon algorithm; in particular, wθ-2(θ) and wφ-2(φ) happen to be the 1D modified cost function of the Capon algorithm. In other words, we utilize the inverse of the revised power spectrum to weight the atoms, with the aim of penalizing directions without sources. Moreover, note that the atomic weight values produced by the log-det function f1(T(u)) correspond to the conventional cost function; this comparison, on the other hand, shows that our surrogate f2(T(u)) performs better.Remarkably, however, our algorithm demonstrates superiority beyond 2D DOA estimation with respect to conventional subspace methods such as MUSIC and ESPRIT. The 2D MUSIC method not only requires prior information on the number of sources but also incurs a heavy computational load for the 2D spectrum search; by contrast, the proposed algorithm can directly estimate continuous-valued angles through two decompositions, and it merely needs prior knowledge of the noise level (moreover, even this is unnecessary, as shown by a later analysis). Meanwhile, the 2D ESPRIT method is suitable only for an array with the property of so-called geometric invariance, whereas the proposed algorithm can be applied to a sparse rectangular array (SRA) with missing array outputs or only a few elements of a URA and still achieves better estimation performance under certain conditions [12, 15, 19, 35]. In addition, far fewer snapshots than the array size are required for exact estimation, whereas this is not true of the subspace methods.


## 3.1. Decoupled Rank Minimization Problem

The atomicl0 norm was proposed in [14] as a measure of the sparsest number of spikes with unknown positions that could compose the observed signal in the frequency domain. In this paper, we employ this norm to express the number of sources with unknown directions composing the array output in the spatial domain and implement 2D DOA estimation by optimizing the atomic representation.More precisely, we merge the signal phase information into the steering vectors and redefine them as follows:(2)aNθk,12ϑkn=exp⁡-j2π·0.5n-1sin⁡θk+12ϑk,n=1,…,N(3)aMφk,12ϑkm=exp⁡j2π·0.5m-1sin⁡φk12ϑk,m=1,…,MRegarding aN(θk,1/2ϑk)aMH(φk,1/2ϑk) as an atom, we define the atom set B as(4)B=aNθ,12ϑaMHφ,12ϑ,θ∈0,2π,φ∈0,2π,ϑ∈0,2π=Aθ,φ,ϑ,θ∈0,2π,φ∈0,2π,ϑ∈0,2πThen, the atomic l0 norm induced by B is given by(5)XB,0=infK:X=∑k=1KskAθk,φk,ϑk,Aθ,φ,ϑ∈BThis is an optimization problem with the objective of finding the fewest possible number of sources given the observed data Xo, which can be represented as a positive combination of the atoms {A(θk,φk,ϑk)}. Inspired by the idea proposed in [14] and previous research on rank minimization [24–26], we will demonstrate that the above norm can be represented as the sum of the ranks of two 1D Toeplitz matrices satisfying a PSD constraint and propose a decoupled rank minimization problem as depicted by the following proposition.Proposition 1.
LetXo∈CN×M be an array observed data that has the following atomic form and satisfies the conditions given in the previous section:(6)Xo=∑k=1KskAθk,φk,ϑkThe corresponding atomic l0 norm minimization problem can be expressed as(7)minXXB,0s.t.X=XoThen, (7) is equivalent to the optimum value of the following decoupled rank minimization problem:(8)minuθ,uφ,X12rankTuθ+12rankTuφs.t.TuφXHXTuθ≥0,X=Xowhere T(uθ) and T(uφ) denote two 1D Toeplitz matrices and uθ and uφ are their first rows, respectively.Proof.
The trivial case ofX=0 is unnecessary to discuss. For X≠0, let r1★ and r2★ be the optimal values of rank(T(uθ)) and rank(T(uφ)), respectively; then, the corresponding optimal value of (8) can be expressed as 1/2r1★+1/2r2★. Furthermore, let the optimal XB,0 be denoted by K★.
On one side, we will show that1/2r1★+1/2r2★≤K★ given problem (7). Suppose that XB,0=K★<min⁡N,M can be realized in the arbitrary atomic decomposed form X=∑k=1K★skaN(θk,1/2ϑk)aMH(φk,1/2ϑk), where {sk}k=1K is a set of positive real numbers. Let uθ=∑kskaN(θk,0) and uφ=∑kskaM(φk,0); then, the homologous Toeplitz matrices can be represented as(9)Tuθ=∑kskaNθk,0aNHθk,0=∑kskaNθk,12ϑkaNHθk,12ϑk≥0and(10)Tuφ=∑kskaMφk,0aMHφk,0=∑kskaMφk,12ϑkaMHφk,12ϑk≥0Thus, we obtain the following PSD matrix:(11)TuφXHXTuθ=∑kK★skaMφk,12ϑkaNφk,12ϑkaMφk,12ϑkaNφk,12ϑkH≥0;i.e., the relational expressions r1★≤K★ and r2★≤K★ hold. Hence, 1/2r1★+1/2r2★≤K★.
On the other side, given problem (8), we will try to verify that K★≤1/2r1★+1/2r2★ in the case of r1★<N and r2★<M (the full-rank case is trivial). Given the PSD condition, and letting uθ★ and uφ★ compose the optimal solution to (8), the following can be derived in accordance with the Schur Complement Lemma [27, 28]:(12)Tuφ★XHXTuθ★≥0⟹Tuφ★≥0IM-Tuφ★Tuφ★+XH=0(13)Tuθ★XXHTuφ★≥0⟹Tuθ★≥0XHIN-Tuθ★Tuθ★+=0In addition, note that T(uθ★) and T(uφ★) have Vandermonde decompositions of the forms T(uθ★)=VθDθVθH and T(uφ★)=VφDφVφH [14], where Vθ=[aN(θ1,0),…,aN(θr1★,0)], Vφ=[aM(φ1,0),…,aM(φr2★,0)], and both Dθ (with dimensions of r1★×r1★) and Dφ (with dimensions of r2★×r2★) are diagonal matrices with positive real elements. Therefore, (12) implies that X lies in the range of Vθ, which means that X can be composed of no more than r1★ atoms; i.e., K★≤r1★. Similarly, (13) leads to the conclusion that K★≤r2★. Thus, K★≤1/2r1★+1/2r2★.
Based on the above proposition, we can straightforwardly optimize the decoupled rank minimization problem (8) under the PSD constraint to implement the optimization of (7). It is evident from the above proof that θk and φk are encoded in the Toeplitz matrices T(uθ★) and T(uφ★), respectively. Once the optimal solution (uθ★,uφ★,X★) is obtained, θk and φk can be determined from the Vandermonde decompositions of these two matrices; there are many off-the-shelf computational methods such as the dual root-finding method [12, 14, 15] and the matrix pencil method [29, 30]. In addition, supposing θk and φk are acquired, the simple and straightforward method of angle pairing provided in [20] can be applied. Therefore, our task is simply to optimize problem (8) for (uθ★,uφ★,X★).


## 3.2. A Novel Sparse Surrogate for Rank Minimization

The essential objective of the decoupled rank minimization problem (8) discussed above is to recover two low-rank matrices within the PSD feasible set. Unfortunately, due to the difficulty of solving this problem (which is similar to that of solving the original optimization problem), there is an urgent need to seek new metrics to express the objective function. The main idea is to find one simpler surrogate function that not only can effectively simulate the ranks of the two Toeplitz matrices but also can simultaneously encourage a low-rank solution. In most papers that discuss the rank minimization problem [16, 27, 31, 32], the log-det function f1(T(u)), which is expressed as follows, is generally used:(14)f1Tu=ln⁡Tu+ε1I=ln⁡1ε1Tu+I+ln⁡ε1However, in this paper, we consider a novel concave surrogate function for T(u)∈CQ×Q, given by (15)f2Tu=-ε2trTu+ε2I-1where I represents the identity matrix. The regularization parameter ε2>0 enables the function f2(T(u)) to approximate rank(T(u)) as it varies and prevents the function from taking a value of -∞ when T(u) is rank-deficient. Moreover, note that the parameter ε1 does the same job with f1(T(u)).From the properties of matrix eigenvalues and trace, it follows that(16)eigTu+ε2I-1=1λi+ε2,i=1,…,Qand(17)f2Tu=-∑ε2λi+ε2=-∑1λi/ε2+1,i=1,…,Qwhere {λi}i=1Q is the set of eigenvalues of T(u). Equation (17) implies that the essence of function f2(T(u)) is to operate on {λi}i=1Q; thus, to better understand the rank-approximation property of f2(T(u)) on a PSD definition domain, we denote a new subfunction by p2(λi)=-(λi/ε2+1)-1+1 for any eigenvalue λi≥0, which is an inverse function obtained through a certain translation and scaling process and has the features of monotonicity and concavity. Note that the constant term 1 is added merely for consistency with the l0 norm; it is equivalent to adding an identity matrix to f2(T(u)) and does not influence the outcome. Obviously, we have the identity p2(0)=0; at λi>0 and near the zero point, p2(λi) has an increasingly steep slope of 1/ε2 as ε2→0, which equivalently imposes a large penalty on near-zero λi and ensures that the values of p2(λi) at these λi become much smaller. In addition, p2(λi) tends towards a constant value of 1 for a sufficiently large λi. The above trend illustrates that as ε2→0, p2(λi) favours the separation of λi at different levels, thus becoming a relaxation of the l0 norm λil0, which has an infinite slope near the zero point and a constant value of 1 except at the zero point. To provide an intuitive illustration, we plot p2(λi) along different ε2’s at λi≥0 in Figure 2. Meanwhile, the curves of the l0 norm, the l1 norm, and the log-det subfunction represented by p1(λi)=ln⁡(λi/ε1+1) are provided for comparison. Suitable proportionality constants are set to ensure that p2(λi) and p1(λi) are zero at λi=0 and have a constant value of 1 at λi=1 without modifying their properties. As shown in Figure 2, both p2 and p1 approach the l0 norm as the regularization parameters decrease; the smaller ε2 is, the steeper the slope is, which offers a basis for the selection of ε2 during both initialization and iteration, as discussed in Section 3.3. Moreover, the function p2 is more likely to perform better since it has a much faster speed of approximation; i.e., p2 is deemed to be more sparsity-enhancing.Figure 2
An illustration of the sparsity-enhancing properties of p1 and the proposed p2, compared with thel0 norm and l1 norm.The approximate relationship betweenp2(λi) and λil0 provides an explanation for why rank(T(u)) can be relaxed by f2(T(u)). In particular, since each larger λi corresponds to a value of p2(λi) that tends towards a constant value of 1 as ε2→0 and the addition operation does not affect the monotonicity and concavity of a function, the sum of all p2(λi) becomes a relaxation of the number of larger eigenvalues, which is considered to be the rank of the Hermitian matrix Tu or the number of sources; i.e., the sparse surrogate f2(T(u)) can be regarded as a relaxation of rank(T(u)). Furthermore, the proposed function f2(T(u)) leads to better performance than f1(T(u)) in encouraging low-rank solutions, as seen by comparing the curves for p2 with the curves for p1; stronger interpretations will be presented in later subsections.


## 3.3. RDAM via the Majorization-Minimization (MM) Algorithm

Using the surrogate functionf2(T(u)), the decoupled rank minimization problem (8) can be relaxed as follows:(18)minuθ,uφ,X-12εθtrTuθ+εθI-1-12εφtrTuφ+εφI-1s.t.TuφXHXTuθ≥0,X=XoHere, the parameters εθ and εφ are deemed to be mutually independent. In this subsection, we aim to minimize the objective function in (18), denoted by F(T(uθ),T(uφ)), on the PSD cone PSD(uθ,uφ). Considering that F(T(uθ),T(uφ)) is a sum of two concave functions, the MM algorithm [33, 34] can be adopted to obtain a local optimum. This algorithm follows an iterative procedure consisting of two steps, i.e., majorization and minimization, in each iteration.Given an initial point{T(uθ,0),T(uφ,0)}∈PSD(uθ,uφ), a series of feasible points {T(uθ,j),T(uφ,j)}j∈N are generated. At each feasible point, the majorization step requires the construction of a tractable function majorizing F(T(uθ),T(uφ)) while ensuring that its value at the current point {T(uθ,j),T(uφ,j)} is equal to F(T(uθ,j),T(uφ,j)). Here, we use the first-order Taylor expansion for this purpose. Since both terms of F(T(uθ),T(uφ)) possess the same structure, we consider only the first term as an example, for which we have the following first-order approximate expression:(19)εθtrTuθ,j+εθI-1+trεθTuθ,j+εθI-2Tuθ-Tuθ,j=εθtr(Tuθ,j+εθI-2Tuθ+cθ,jwhere cθ,j is a constant. The derivation is presented in Appendix B. Thus, F(T(uθ),T(uφ)) has the following upper bound:(20)FTuθ,Tuφ≤εθ2tr(Tuθ,j+εθI-2Tuθ+εθ2tr(Tuφ,j+εφI-2Tuφ+Cj=ΔGTuθ,Tuφ:Tuθ,j,Tuφ,jwhere Cj=cθ,j+cφ,j is constant. It is evident that GTuθ,j,Tuφ,j:Tuθ,j,Tuφ,j=FTuθ,j,Tuφ,j. In addition, many other construction methods [33] (e.g., second-order Taylor expansion, inequality, and the Schur complement) could be used; these possibilities are left for future work.Then, in the minimization step, the following minimization problem is solved, ignoring the constantCj, to update T(uθ,j+1) and T(uφ,j+1):(21)minuθ,j+1,uφ,j+1,Xj+112tr(εθTuθ,j+εθI-2Tuθ+12tr(εθTuφ,j+εφI-2Tuφs.t.TuφXHXTuθ≥0,X=Xjwhere the optimization problem is a semidefinite programming (SDP) that can be solved in polynomial time. Notice that this problem is equivalent to a weighted DAM problem in the PSD form [20], with weighting factors of Wθ,j=εθ(T(uθ,j)+εθI)-2 and Wφ,j=εθ(T(uφ,j)+εφI)-2, and that the proposed algorithm is iterative in nature; therefore, our algorithm is called RDAM. In addition, we have the following inequality: (22)FTuθ,j+1,Tuφ,j+1≤GTuθ,j+1,Tuφ,j+1:Tuθ,j,Tuφ,j+Cj≤GTuθ,j,Tuφ,j:Tuθ,j,Tuφ,j+Cj=FTuθ,j,Tuφ,j;i.e., {F(T(uθ,j),T(uφ,j))}j∈N is a nonincreasing sequence, implying that the iterative process makes the objective function F(T(uθ),T(uφ)) repeatedly smaller, and it can be expected to converge to a local minimum F(T(uθ★),T(uφ★)), as proven in detail in [21, 33]. Furthermore, since the algorithm is locally convergent, suitable initial values are indispensable for guaranteeing convergence to the desired local optimum. One immediate possibility is to use the DAM solution for initialization. In particular, let Wθ,0 and Wφ,0 be identity matrices; then, we can obtain the initial values T(uθ,1) and T(uθ,1) by solving a DAM problem. Next, let us emphasize the selection of εθ,1 and εφ,1, both of which influence the outcome. As argued in Section 3.2, F(T(uθ),T(uφ)) is the sum of two subfunctions, denoted by p2(λθ,i)=-(λθ,i/εθ+1)-1 and p2(λφ,i)=-(λφ,i/εφ+1)-1, for θ and φ, respectively. For consistency with the curves in Figure 2, we set εθ=aλθ,max and εφ=bλφ,max. Then, both the factors a and b can be viewed as equivalent to the parameter ε2 in this figure, and we have λθ,i/λθ,max∈[0,1] and λφ,i/λφ,max∈[0,1]. Now, let us consider the choice of a as an example. Let larger eigenvalues, corresponding to the signals, be denoted by λθ,si, while smaller eigenvalues, corresponding to the noise, are denoted by λθ,ni; thus, we should be able to ensure a relatively large difference in the values of p2(λθ,i) between the points λθ,si/λθ,max and λθ,ni/λθ,max by tuning the initial factor a. For example, for points λθ,si/λθ,max and λθ,ni/λθ,max that are close to each other, we can select a small a value, such as 0.1 or smaller; then, there will be a sharp slope such that the subfunction values will be distinctly separated. In contrast, for the relatively distant points, a large a value, such as 1, will be sufficient to achieve good estimation performance. The same analysis applies for the initialization of the factor b.Finally, it is instructive to explore the potential implications of the RDAM algorithm in combination with the atomic norm [13, 14]. Notice that the concept of the 1D weighted atom set proposed in [16] by Yang et al. can be extended to a new 2D weighted atom set Bw as follows:(23)Bw=wθθaNθ,12ϑwφ∗φaMHφ,12ϑ,θ,φ∈0,2π×0,2π=Awθ,φ,ϑ,θ,φ∈0,2π×0,2πwhere {wθ(θ),wφ∗(φ)} denotes a pair of nonnegative weights on the 2D space angles. Accordingly, the weighted atom norm is given by(24)XBw=inf⁡t>0:X∈tconvAw=inf⁡∑kskw:X=∑kskwAwθk,φk,ϑk,Awθ,φ,ϑ∈Bw=inf⁡∑kskw:X=∑kskwwθθkwφ∗φkaNθk,12ϑkaMHφk,12ϑk=inf⁡∑kskwθθkwφ∗φk:X=∑kskAθk,φk,ϑk,Aθk,φk,ϑk∈BThen, a 2D weighted ANM problem in decoupled PSD form can be obtained based on the following proposition.Proposition 2.
LetXBw,j+1 denote the (j+1)th 2D weighted atomic norm. The following equivalence relation holds given the data Xj∈CN×M recovered from the jth iteration.(25)min⁡XBw,j+1=minuθ,j+1,uφ,j+112trWθ,jTuθ+12trWφ,jTuφs.t.TuφXHXTuθ≥0,X=Xjwhere(26)wθ,jθ=1aNHθWθ,jaNθandwφ,jφ=1aMHφWφ,jaMφ.This implies that iteratively minimizing the weighted atomic normXBw induced by the 2D weighted atom set Bw serves the same purpose as the proposed RDAM algorithm; i.e., RDAM can also be called decoupled reweighted atom norm minimization (DRAM). More significantly, our proposed trace function f2(T(u)) gives rise to quadratic weighting factors, which greatly accelerates the speed of local convergence and further enhances the preference for low-rank solutions in (18) compared with the first-order factors generated by the log-det function f1(T(u)) [16].


## 3.4. Another Interpretation from the Perspective of Conventional Direction of Arrival (DOA) Estimation Methods

The proposed RDAM algorithm can also be interpreted from the perspective of traditional DOA estimation methods. Suppose that the signals are independent of each other and that all steering vectors have been normalized, without exception. Then, we can redefine two decoupled covariance matrices for the noiseless observation data as follows:(27)R~θ≜XXH=∑k=1KaNθkaMHφksktsk∗taMφkaNHθk=∑k=1Kskt2aNθk,12ϑkaNHθk,12ϑk(28)R~φ≜XHX=∑k=1KaMφkaNHθksktsk∗taNθkaMHφk=∑k=1Kskt2aMφk,12ϑkaMHφk,12ϑkFrom the poof of Proposition 1, we know that the Toeplitz matrices T(uθ) and T(uφ) are exactly equivalent to the two decoupled covariance matrices above except that the coefficients include a squared difference. The whole RDAM process, therefore, can be viewed as the process of recovering R~θ and R~φ with the regulators εθ and εφ. Furthermore, the two atomic weight values in (24) are related to the Capon algorithm; in particular, wθ-2(θ) and wφ-2(φ) happen to be the 1D modified cost function of the Capon algorithm. In other words, we utilize the inverse of the revised power spectrum to weight the atoms, with the aim of penalizing directions without sources. Moreover, note that the atomic weight values produced by the log-det function f1(T(u)) correspond to the conventional cost function; this comparison, on the other hand, shows that our surrogate f2(T(u)) performs better.Remarkably, however, our algorithm demonstrates superiority beyond 2D DOA estimation with respect to conventional subspace methods such as MUSIC and ESPRIT. The 2D MUSIC method not only requires prior information on the number of sources but also incurs a heavy computational load for the 2D spectrum search; by contrast, the proposed algorithm can directly estimate continuous-valued angles through two decompositions, and it merely needs prior knowledge of the noise level (moreover, even this is unnecessary, as shown by a later analysis). Meanwhile, the 2D ESPRIT method is suitable only for an array with the property of so-called geometric invariance, whereas the proposed algorithm can be applied to a sparse rectangular array (SRA) with missing array outputs or only a few elements of a URA and still achieves better estimation performance under certain conditions [12, 15, 19, 35]. In addition, far fewer snapshots than the array size are required for exact estimation, whereas this is not true of the subspace methods.


## 4. Further Analysis

### 4.1. Superresolution

The DAM algorithm requires the following angular interval conditions [20] to guarantee that directions can be estimated from the observed data X∗ with probability 1-δ, where δ is the estimation precision.(29)Δθ=mini≠j⁡sin⁡θi-sin⁡θj≥1.19N-1/4(30)Δφ=mini≠j⁡sin⁡φi-sin⁡φj≥1.19M-1/4Here, Δθ and Δφ are the wrap-around distances on the complex unit circle. These conditions mean that the resolutions in terms of θ and φ are inversely proportional to the numbers of sensors by row and column, respectively. By contrast, the proposed algorithm (RDAM) has no such restriction on the angular intervals. As the analysis in Section 3 shows, the quadratic inverse of the solution from the previous iteration, denoted by Wθ,j=εθ(T(uθ,j)+εθI)-2 and Wφ,j=εθ(T(uφ,j)+εφI)-2, is utilized as the weighting factors for the objective function (21) in the current iteration; thus, large weights are applied to small eigenvalues, while small weights are applied to large eigenvalues. Therefore, minimizing (21) is equivalent to penalizing small eigenvalues, making them much smaller while making large ones relatively larger. Thus, it becomes far easier to pick out the eigenvalues corresponding to signals. In other words, the iterative nature of the algorithm allows it to overcome the limits of the minimal angular interval, and, hence, the resolution is greatly improved. In addition, the simulations presented in the next section similarly indicate that, given a proper initialization, our algorithm yields a very high resolution.


### 4.2. Computational Complexity

Another merit of RDAM is its low complexity. As mentioned in [18], the vector-based atomic norm minimization (vector-ANM) method has a computational complexity as high as o((NM)3.5log⁡(1/δ)) due to its PSD constraint of size (NM+1)×(NM+1). However, in the DAM method [20], the two-level Toeplitz matrix is divided into two one-level Toeplitz matrices, thus reducing the computational complexity to o((N+M)3.5log⁡(1/δ)).Fortunately, our proposed RDAM algorithm has the same complexity as DAM in a single iteration; thus the computational efficiency depends merely on the number of iterations. Since such a reweighted iterative algorithm can converge within a few steps [21], the computational cost of RDAM is no greater than d·o((N+M)3.5log⁡(1/δ)), where d denotes the number of iterations; i.e., the complexity is comparable to that of DAM but far below the heavy load of vector-ANM. An empirical run-time comparison is also provided based on the simulations reported below.


### 4.3. The Noisy Case

In fact, DOA estimation is typically implemented in a noisy environment; therefore, we consider the following observation model with additive noiseN, which involves independent entries drawn from N(0,σ2):(31)Yo=X+N=∑k=1KskaNθkaMHφk+σ2IThen, an iterative denoising optimization problem is obtained as follows:(32)minuθ,j+1,uφ,j+1,Xj+1trWθ,jTuθ+trWφ,jTuφ+λYo-XF2s.t.TuφXHXTuθ≥0where λ is a regularization factor that keeps a balance between data fitting and a low-rank solution. When λ≥EM2, where ·2 represents the spectral norm, the soft threshold denoising algorithm has the following asymptotic mean square error (MSE) [36]:(33)1NM2X-XoF2=oσ2KNMIn this paper, we directly optimize the dual problem [37] of (32) as follows:(34)minAj+1,Zj+1,Vj+12RtrYoHV+2ηVFs.t.AVHVZ≥0∑n=1N-lZn,n+l=∑n=1N-lWθ,jn,n+l,l=0,…,N-1∑m=1M-kAm,m+k=∑m=1M-kWφ,jm,m+k,k=0,…,M-1where η≥Yo-XF. In practice, we can set an initial value for η and then gradually decrease it by a suitable small factor in each iteration; as shown by the simulations reported below, this approach results in better performance.


## 4.1. Superresolution

The DAM algorithm requires the following angular interval conditions [20] to guarantee that directions can be estimated from the observed data X∗ with probability 1-δ, where δ is the estimation precision.(29)Δθ=mini≠j⁡sin⁡θi-sin⁡θj≥1.19N-1/4(30)Δφ=mini≠j⁡sin⁡φi-sin⁡φj≥1.19M-1/4Here, Δθ and Δφ are the wrap-around distances on the complex unit circle. These conditions mean that the resolutions in terms of θ and φ are inversely proportional to the numbers of sensors by row and column, respectively. By contrast, the proposed algorithm (RDAM) has no such restriction on the angular intervals. As the analysis in Section 3 shows, the quadratic inverse of the solution from the previous iteration, denoted by Wθ,j=εθ(T(uθ,j)+εθI)-2 and Wφ,j=εθ(T(uφ,j)+εφI)-2, is utilized as the weighting factors for the objective function (21) in the current iteration; thus, large weights are applied to small eigenvalues, while small weights are applied to large eigenvalues. Therefore, minimizing (21) is equivalent to penalizing small eigenvalues, making them much smaller while making large ones relatively larger. Thus, it becomes far easier to pick out the eigenvalues corresponding to signals. In other words, the iterative nature of the algorithm allows it to overcome the limits of the minimal angular interval, and, hence, the resolution is greatly improved. In addition, the simulations presented in the next section similarly indicate that, given a proper initialization, our algorithm yields a very high resolution.


## 4.2. Computational Complexity

Another merit of RDAM is its low complexity. As mentioned in [18], the vector-based atomic norm minimization (vector-ANM) method has a computational complexity as high as o((NM)3.5log⁡(1/δ)) due to its PSD constraint of size (NM+1)×(NM+1). However, in the DAM method [20], the two-level Toeplitz matrix is divided into two one-level Toeplitz matrices, thus reducing the computational complexity to o((N+M)3.5log⁡(1/δ)).Fortunately, our proposed RDAM algorithm has the same complexity as DAM in a single iteration; thus the computational efficiency depends merely on the number of iterations. Since such a reweighted iterative algorithm can converge within a few steps [21], the computational cost of RDAM is no greater than d·o((N+M)3.5log⁡(1/δ)), where d denotes the number of iterations; i.e., the complexity is comparable to that of DAM but far below the heavy load of vector-ANM. An empirical run-time comparison is also provided based on the simulations reported below.


## 4.3. The Noisy Case

In fact, DOA estimation is typically implemented in a noisy environment; therefore, we consider the following observation model with additive noiseN, which involves independent entries drawn from N(0,σ2):(31)Yo=X+N=∑k=1KskaNθkaMHφk+σ2IThen, an iterative denoising optimization problem is obtained as follows:(32)minuθ,j+1,uφ,j+1,Xj+1trWθ,jTuθ+trWφ,jTuφ+λYo-XF2s.t.TuφXHXTuθ≥0where λ is a regularization factor that keeps a balance between data fitting and a low-rank solution. When λ≥EM2, where ·2 represents the spectral norm, the soft threshold denoising algorithm has the following asymptotic mean square error (MSE) [36]:(33)1NM2X-XoF2=oσ2KNMIn this paper, we directly optimize the dual problem [37] of (32) as follows:(34)minAj+1,Zj+1,Vj+12RtrYoHV+2ηVFs.t.AVHVZ≥0∑n=1N-lZn,n+l=∑n=1N-lWθ,jn,n+l,l=0,…,N-1∑m=1M-kAm,m+k=∑m=1M-kWφ,jm,m+k,k=0,…,M-1where η≥Yo-XF. In practice, we can set an initial value for η and then gradually decrease it by a suitable small factor in each iteration; as shown by the simulations reported below, this approach results in better performance.


## 5. Numerical Simulations

In this section, we carry out numerical simulations conducted to illustrate the performance of the proposed superresolution algorithm (RDAM) in 2D DOA estimation for a rectangular array. In all experiments, the value ofK was not given a priori, and a random phase (uniformly distributed between 0 and 2π), associated with a constant magnitude of 1, was set for each received signal. In addition, all SDPs involved were implemented using CVX [38], and the matrix pencil method was adopted for the Vandermonde decomposition of the Toeplitz matrices. In addition, the RDAM algorithm was initialized as discussed in Section 3, and we not only selected suitable initial factors a and b, as required, but also iteratively tuned them by decreasing their values by a factor of 0.5 in each subsequent iteration.First, let us present an intuitive illustration of the higher resolution of RDAM compared with DAM. Suppose thatK=5 narrowband signals with directions of (10.1°,56.12°), (16.52°,65.08°), (25.01°,40.25°), (34.07°,32.25°), and (40.15°,47.01°) impinge on a rectangular array of size N×M=10×10. Note that both max⁡{Δθ} and max⁡{Δφ} are no larger than 10°. Figure 3(a) shows the DAM estimation result, while Figure 3(b) shows the performance of RDAM with initial factors a=b=1; all results were obtained through 100 Monte Carlo experiments. It is apparent that, in the noiseless case, all the angles are estimated correctly and stably by our algorithm, whereas DAM yields unsatisfactory outcomes. Moreover, RDAM is even robust to changes in the signal phases, possibly extending to certain special fields.Figure 3
An intuitive illustration of 2D DOA estimation. (a) DAM algorithm; (b) proposed algorithm.
(a) (b)Second, let us compare the proposed algorithm with three other methods, namely, DAM, vector-ANM, and vector-based reweighted atomic norm minimization (vector-RAM) [18], to provide a more reliable verification of the superresolution of our algorithm. In particular, let K=2 and N=M=10. We fix the direction of one of the two signals as (30°,60°), while the other signal is radially separated from the fixed signal by azimuth and elevation spacings that both range from 2° to 20°. Each interval value is allowed to randomly fluctuate within a range of ±0.025°. Let a=b=1, and let the number of iterations of RDAM be four; the same parameter settings are considered for vector-RAM. Figure 4 compares the results of 500 Monte Carlo experiments each for the azimuth θ and the elevation. It is clear that the proposed algorithm can distinguish spacing angles that differ by 3° in both dimensions, while the other algorithms cannot. Moreover, it achieves better performance than vector-RAM, which is a reweighted algorithm based on stretched vectors, under the same initial conditions. In addition, we list the run-times of these four methods in Table 1, from which we can see that, for both vector-based methods, the computational complexity increases dramatically as the size of the array increases. By contrast, the proposed algorithm and the DAM algorithm both have very low complexity regardless of the array size.Table 1
Run-times of four methods.
Algorithm N = M = 10 N = M = 15 N = M = 20 vector-ANM 3.7192 30.8332 293.0253 vector-RAM 11.4040 140.8327 1276.8 DAM 0.3402 0.3855 0.7388 RDAM 1.3968 1.8500 3.1428Figure 4
Advantages in terms of resolution for two angles, respectively. (a) Azimuthθ; (b) elevation φ.
(a) (b)Furthermore, we examine the iterative process of RDAM to investigate the hypothesis that the proposed surrogate functionf2(T(u)) enables the number of iterations to be reduced, i.e., more effectively encourages low-rank solutions and enables a faster convergence speed. Accordingly, let K=3 and N=M=10, and consider randomly selected signal directions of (34.35°,67.03°), (38.09°,76.04°), and (44.12°,70.56°). Initial factors of 0.01 and 0.1 are selected for the functions f1(T(u)) and f2(T(u)), respectively; under these conditions, the subfunctions have similar approximation properties, as seen from Figure 2. The iterative processes using the two surrogate functions are shown in Figure 5, where the first two subplots in both (a) and (b) present the changes in the eigenvalues for θ and φ, respectively, and the third subplot shows the estimation performance in each iteration. It is obvious that the previous function f1(T(u)) requires four iterations wihle the proposed function f2(T(u)) needs only three iterations, which indicates that the proposed function has a faster estimation speed.Figure 5
Advantages of convergence speed: (a) log-det function; (b) proposed function.
(a) (b)Finally, we carry out 200 Monte Carlo experiments considering the case of very closely spaced sources to illustrate the antinoise performance of the proposed algorithm relative to the Cramér-Rao bound (CRB). The array response with i.i.d. additive Gaussian noise satisfyingN(0,σ2) was obtained; the signal-to-noise ratio (SNR) was 10log⁡σ-2 since the powers of all signals were set to a constant value of 1. As before, we set N=M=10, and the directions of the two sources were set to (10.10°,40.01°) and (13.56°,43.12°). Figure 6 shows a graph of the average root mean square error (RMSE) with respect to the SNR with the parameters set to a=b=1 and the initial noise level set to NMσ2. It is worth noting that the noise level was reduced by a factor of 0.8 in the first weighted iteration for both our RDAM algorithm and the vector-RAM algorithm. Obviously, the performance of the proposed algorithm is close to the CRB at higher SNRs, although its antinoise performance may be worse below SNR=15dB; the main reason for the performance degradation is that, at lower SNRs, the probability of efficient estimation decreases so that a very small amount of deteriorating estimates is obtained, which influences the performance of the curve. However, our algorithm still shows more robust performance than the others in the absence of noise.Figure 6
Noise performance: RMSE versus SNR.


## 6. Conclusion

In this paper, we proposed a superresolution decoupled reweighted atomic norm minimization (RDAM) algorithm for 2D DOA estimation based on the popular theory of continuous compressed sensing (CCS). By directly employing the atomicl0 norm and converting the optimization problem into a decoupled rank minimization problem, we ensured that the computational complexity is greatly reduced. On the other side, we presented a novel surrogate function and formulated a new concave optimization problem, overcoming the shortcomings of existing methods in terms of resolution and also accelerating the convergence speed. Our algorithm is even suitable for application to sparse rectangular array (SRA), although there is a need for better denoising processing, which remains to be studied in future work.



---

*Source: 6797168-2019-07-08.xml*