Direction of Arrival Estimation in Elliptical Models via Sparse Penalized Likelihood Approach

Abstract

In this paper, an $l_1$-penalized maximum likelihood (ML) approach is developed for estimating the directions of arrival (DOAs) of source signals from the complex elliptically symmetric (CES) array outputs. This approach employs the $l_1$-norm penalty to exploit the sparsity of the gridded directions, and the CES distribution setting has a merit of robustness to the uncertainty of the distribution of array output. To solve the constructed non-convex penalized ML optimization for spatially either uniform or non-uniform sensor noise, two majorization-minimization (MM) algorithms based on different majorizing functions are developed. The computational complexities of the above two algorithms are analyzed. A modified Bayesian information criterion (BIC) is provided for selecting an appropriate penalty parameter. The effectiveness and superiority of the proposed methods in producing high DOA estimation accuracy are shown in numerical experiments.

1. Introduction

Estimating the directions of arrival (DOAs) of a number of far-field narrow-band source signals is an important problem in signal processing. Many DOA estimation methods were proposed early on, such as multiple signal classification (MUSIC) [1], estimation of signal parameters via rotation invariance techniques (ESPRIT) [2] and their variants [3,4]. Many of them work well when having accurate estimates of the array output covariance matrix and source number. In scenarios with sufficient array snapshots and a moderately high signal-to-noise (SNR), the array output covariance matrix and source number can be accurately estimated.

Recently, some sparse DOA estimation methods are popularly proposed based on sparse constructions of the array output model. They are applicable if the number of sources is unknown and many of them are effective when with a limited number of data snapshots. In [5], the sparse DOA estimation methods are categorized into three groups: the on-grid, the off-grid and the grid-less. The on-grid method, which is widely studied and straightforward to implement, assumes that the true DOAs are on a predefined grid [6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22]. The off-grid method also uses a prior grid but does not constrain the DOA estimates to be on this grid, while it introduces more unknown parameters to be estimated and complicates the algorithm [23,24,25,26,27,28]. The grid-less method directly operates in the entire direction domain without a pre-specified grid but, currently, is developed mainly for linear array and encounters rather large computational burden [29,30,31]. Although the on-grid methods may induce the grid mismatch, it is still attractive due to its easy accessibility in a wide range of applications (see, e.g., [28,32,33]). To alleviate the grid mismatch, methods for selecting an appropriate prior grid are proposed in [9,34].

During the past two decades, various on-grid sparse techniques, such as linear least-square, covariance-based and maximum-likelihood (ML) type methods, were researched. The linear least-square methods [6,9,11,12], minimizing the $l_2$-norm of the noise residual of the deterministic model, enforce the sparsity by constraining the $l_p$-norm of the signal vector for $p\in [0,1]$ (note that $l_p$-norm for $p\in [0,1)$ is defined as the same form as the one for $p\ge 1$ but is not truly a norm in mathematics). The covariance-based methods, such as the sparse iterative covariance-based estimation (SPICE) method [13,14] and its variants [15,16], and those proposed in [17,18,19,20,21,22], are derived from different covariance-fitting criteria and use the $l_p$-norm penalty to enforce the sparsity. The ML type methods, including the sparse Bayesian learning (SBL) methods [7,8,10] and the likelihood-based estimation of sparse parameters (LIKES) method [15,16], are deduced under the assumption that the array output signal is multivariate Gaussian distributed. The SBL methods use different prior distributions to model the sparsity. The LIKES method appears to be not explicitly utilizing the sparsity, but provides sparse estimates.

The ML principle is generally believed to be statistically more sound than the covariance-fitting principle [15]. The existing on-grid ML methods are developed on the Gaussian distribution assumption that are not satisfied in many signal processing practical applications. The complex elliptically symmetric (CES) distributions, containing the t-distribution, the K-distribution, the Gaussian distribution and so on, can be used to characterise the Gaussian and non-Gaussian random variables. Particularly, the heavy-tailed data, which usually appear in the field of radar and array signal processing [35,36], can be modeled by some CES distributions.

The $l_p$-norm penalty with $p\in [0,1]$ is well-known for its contributions in deriving sparse solutions [37]. In particular, the convex $l_1$-norm penalty, also known as the least absolute shrinkage and selection operator (LASSO) penalty, is easy to handle and widely used in various applications. Generally, a sparse and accurate estimate of the unknown sparse parameter can be expected when minimizing the negative likelihood plus an $l_p$-norm penalty.

For estimating the DOAs from the CES distributed array outputs, we provide an $l_1$-norm penalized ML method based on a sparse reconstruction of the array output covariance matrix. The characteristics and advantages of our method are explained as follows:

• Our method is a sparse on-grid method based on the penalized ML principle, and is designed especially for the CES random output signals. The Gaussian distribution and many heavy-tailed and light-tailed distributions are included in the class of the CES distributions, and it is worth mentioning that the heavy-tailed output signals that are non-Gaussian are common in the field of signal processing [35,36].
• The sparsity of the unknown sparse vector is enforced by penalizing its $l_1$-norm. When the penalty parameter becomes zero, the proposed method is actually the general ML DOA estimation method that is applicable to the scenarios with CES random output signals.
• Two penalized ML optimization problems are formulated for spatially uniform and non-uniform white sensor noise, respectively. Since it is difficult to solve the two non-convex penalized ML optimizations globally, two majorization-minimization (MM) algorithms having different iterative procedures are developed for seeking the optimal solutions locally for each of them. Some discussions on the computational complexities of the above two algorithms are provided. In addition, the optimal penalty parameter is suggested.
• The proposed methods are evaluated numerically in scenarios with Gaussian and non-Gaussian output signals. Particularly, the performance gains originated from the added $l_1$-norm penalty are numerically demonstrated.

The remainder of this paper is organized as follows. Section 2 introduces a sparse CES data model of the array output. In Section 3, a sparse penalized ML method is developed to estimate the DOAs. For solving the proposed $l_1$-penalized ML optimizations that are non-convex, algorithms in the MM framework are developed in Section 4. Section 5 numerically shows the performance of our method in Gaussian and non-Gaussian scenarios. Finally, some conclusions are given in Section 6.

Notations

The notation ${\mathbb{C}}^m\left({\mathbb{R}}^m\right)$ denotes the set of all m-dimensional complex- (real-) valued vectors, and ${\mathbb{C}}^{m\times n}\left({\mathbb{R}}^{m\times n}\right)$ denotes the set of all $m\times n$ complex- (real-) valued matrices. ${\mathbf{1}}_p$ and ${\mathbf{0}}_p$ are the p-dimensional vectors with all elements equal to 1 and 0, respectively. ${\mathbf{I}}_p$ is the $p\times p$ identity matrix. ${\parallel ·\parallel }_1$ and ${\parallel ·\parallel }_2$ denote the $l_1$-norm and $l_2$-norm of a vector, respectively. The superscripts ${(·)}^T$ and ${(·)}^H$, respectively, denote the transpose and the conjugate transpose of a vector or matrix. The imaginary unit is denoted as ı defined by ${ı}^2=-1$.

For a vector $\mathbf{x}={[x_1,\dots ,x_p]}^T\in {\mathbb{C}}^p$, define element-wise square root function $sqrt\left(\mathbf{x}\right)={[\sqrt{|x_1|},\dots ,\sqrt{|x_p|}]}^T$ and element-wise division operation $\mathbf{x}\oslash \mathbf{y}={[x_1/y_1,\dots ,x_p/y_p]}^T$ for a vector $\mathbf{y}={[y_1,\dots ,y_p]}^T\in {\mathbb{C}}^p$ with non-zero elements. $max\left(\mathbf{x}\right)=max\{x_1,\dots ,x_p\}$, ${\mathbf{x}}_{[a:b]}={[x_a,\dots ,x_b]}^T$ for $1\le a\le b\le p$, $\mathbf{x}>0(\ge 0)$ means $x_i>0(\ge 0)$ for $i=1,\dots ,p$, $[\mathbf{x};\mathbf{z}]$ denotes the stacked vector of the column vectors $\mathbf{x}$ and $\mathbf{z}$. $Diag\left(\mathbf{x}\right)$ denotes a square diagonal matrix with the elements of vector $\mathbf{x}$ on the main diagonal, and ${(Diag\left(\mathbf{x}\right))}^{-1}$ denotes the diagonal matrix with main diagonal elements $1/x_i(i=1,\dots ,p)$ by taking $1/0=\infty$.

For a square matrix $\mathbf{X}\in {\mathbb{C}}^{p\times p}$, ${\mathbf{X}}_{i,j}$ denotes the $(i,j)$th entry of $\mathbf{X}$, $\mathbf{X}>0(\ge 0)$ means that $\mathbf{X}$ is Hermitian and positive definite (semidefinite), $tr\left(\mathbf{X}\right)$ denotes the trace of $\mathbf{X}$, and $diag\left(\mathbf{X}\right)$ denotes a column vector of the main diagonal elements of $\mathbf{X}$.

2. Problem Formulation

Consider the problem of estimating the DOAs of $k_0$ narrow-band signals impinging on an array of m sensors.

Given a set of grid points

$$\{{\theta }_1,\cdots ,{\theta }_k\},$$

(1)

we assume that the true $k_0$ ($k_0\ll k$) DOAs, respectively, denoted as ${\xi }_1,\cdots ,{\xi }_{k_0}$, take values in it. The array output measurement at the instant t, denoted as $\mathbf{x}\left(t\right)\in {\mathbb{C}}^m$, can be modeled as

$$\mathbf{x}\left(t\right)=\mathbf{A}\left(\mathbf{\theta }\right)\mathbf{s}\left(t\right)+\mathbf{v}\left(t\right),$$

(2)

where

• $\mathbf{\theta }={[{\theta }_1,\dots ,{\theta }_k]}^T$,
• $\mathbf{A}\left(\mathbf{\theta }\right)=[\mathbf{a}\left({\theta }_1\right),\dots ,\mathbf{a}\left({\theta }_k\right)]\in {\mathbb{C}}^{m\times k}$ is the known array manifold matrix with $\mathbf{a}\left({\theta }_i\right)$ being the steering vector corresponding to ${\theta }_i$, $i=1,\dots ,k$;
• $\mathbf{s}\left(t\right)={[s_1\left(t\right),\dots ,s_k\left(t\right)]}^T\in {\mathbb{C}}^k$ is the source signal vector at the time instant t, in which $s_i\left(t\right)$ is the unknown random signal from a possible source at ${\theta }_i$ and then is zero if ${\theta }_i$ is not in the true DOA set $\{{\xi }_1,\cdots ,{\xi }_{k_0}\}$, $i=1,\dots ,k$;
• $\mathbf{v}\left(t\right)={[v_1\left(t\right),\dots ,v_m\left(t\right)]}^T\in {\mathbb{C}}^m$ is the noise vector impinging on the sensor array at the time instant t.

Some necessary statistical assumptions are made as follows:

• The possible source signals $s_1\left(t\right),\dots ,s_k\left(t\right)$ are uncorrelated and zero-mean at any time instant t.
• The noise components $v_1\left(t\right),\dots ,v_m\left(t\right)$ are uncorrelated, zero-mean, and independent of $s_1\left(t\right),\dots ,s_k\left(t\right)$ for any time instant t.
• The $n(n>m)$ snapshots $\mathbf{x}\left(1\right),\dots ,\mathbf{x}\left(n\right)$ of sensor array signals are independent and identically distributed from a CES distribution.

Note that the zero-mean assumptions above are common in the signal processing literature [5,35]. The CES distribution setting of the array output $\mathbf{x}\left(t\right)$ enables us to effectively process the Gaussian, heavy-tailed or light-tailed array snapshots, because the class of CES distributions [38] includes the Gaussian distribution, the t-distribution, the K-distribution and so on.

For the simplicity of the notations, we denote $\mathbf{A}\left(\mathbf{\theta }\right)$ as $\mathbf{A}$ and $\mathbf{a}\left({\theta }_i\right)$ as ${\mathbf{a}}_i$ for $i=1,\dots ,k$ in the following. Under the above assumptions, we can find that the array output $\mathbf{x}\left(t\right)$ in Equation (2) at any time instant t has mean zero and covariance matrix

$$\begin{array}{c}\hfill \mathbf{R}=\mathbb{E}[\mathbf{x}\left(t\right)\mathbf{x}{\left(t\right)}^H]=\mathbf{A}\mathbf{P}{\mathbf{A}}^H+\mathbf{V},\end{array}$$

(3)

where

$$\begin{array}{c}\hfill \mathbf{P}=\mathbb{E}[\mathbf{s}\left(t\right)\mathbf{s}{\left(t\right)}^H]=Diag\left(\mathbf{p}\right)\end{array}$$

(4)

is the (unknown) source signal covariance matrix with signal power vector $\mathbf{p}={[p_1,\dots ,p_k]}^T$, and

$$\begin{array}{c}\hfill \mathbf{V}=\mathbb{E}\left[\mathbf{v}\left(t\right)\mathbf{v}{\left(t\right)}^H\right]=Diag\left(\mathbf{\sigma }\right)\end{array}$$

(5)

is the (unknown) noise covariance matrix with the noise power vector $\mathbf{\sigma }={[{\sigma }_1^2,\dots ,{\sigma }_m^2]}^T$. The matrix $\mathbf{R}$ can be rewritten as

$$\begin{array}{c}\hfill \mathbf{R}=[\mathbf{A},{\mathbf{I}}_m]Diag\left([\mathbf{p};\mathbf{\sigma }]\right){[\mathbf{A},{\mathbf{I}}_m]}^H.\end{array}$$

(6)

For any $i=1,\dots ,k$, the signal power $p_i=0$ if ${\theta }_i$ is not in the set of the true DOAs. Therefore, the true DOAs can be identified by the locations of nonzero elements of the power vector $\mathbf{p}$. In the following, the DOA estimation problem is formulated as a problem of estimating the locations of nonzero elements of the power vector $\mathbf{p}$.

3. Sparse DOA Estimation

The array output $\mathbf{x}\left(t\right)$ in Equation (2) is CES distributed with mean zero and covariance matrix $\mathbf{R}$, and then the normalized random vector

$$\begin{array}{c}\hfill \mathbf{y}\left(t\right)\triangleq \mathbf{x}\left(t\right)/{\parallel \mathbf{x}\left(t\right)\parallel }_2,\end{array}$$

(7)

which actually refers to the angle of the array output vector $\mathbf{x}\left(t\right)$, has a complex angular central Gaussian (ACG) distribution [35] with the probability density function

$$\begin{array}{c}\hfill p(\mathbf{y}\left(t\right);\mathbf{R})\propto det{\left(\mathbf{R}\right)}^{-1}{(\mathbf{y}{\left(t\right)}^H{\mathbf{R}}^{-1}\mathbf{y}\left(t\right))}^{-m}.\end{array}$$

(8)

Denote

$$\begin{array}{c}\hfill L_0\left(\mathbf{R}\right)=logdet\left(\mathbf{R}\right)+\frac{m}{n}\sum _{t=1}^{n}log(\mathbf{x}{\left(t\right)}^H{\mathbf{R}}^{-1}\mathbf{x}\left(t\right)),\end{array}$$

(9)

and then the negative log-likelihood of $\mathbf{y}\left(t\right)=\mathbf{x}\left(t\right)/{\parallel \mathbf{x}\left(t\right)\parallel }_2,t=1,\dots ,n$, becomes

$$\begin{array}{c}\hfill n{L}_0\left(\mathbf{R}\right)-m\sum _{t=1}^n{log(\parallel \mathbf{x}\left(t\right)\parallel }_2^2)+c_1,\end{array}$$

(10)

where $c_1$ is a constant.

3.1. Spatially Non-Uniform White Noise

In the case of spatially non-uniform white sensor noise, not all noise variances ${\sigma }_i^2$, $i=1,\cdots ,m$, are equal. Assuming ${\sigma }_m^2>0$, we denote

$$\begin{array}{cc}\hfill \mathbf{r}& =[\mathbf{p};{\mathbf{\sigma }}_{[1:m-1]}]/{\sigma }_m^2,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \mathbf{W}& =\mathbf{B}Diag\left(\mathbf{r}\right){\mathbf{B}}^H+{\mathbf{J}}_m,\hfill \end{array}$$

(12)

where $\mathbf{B}\in {\mathbb{C}}^{m\times (k+m-1)}$ is the first $(k+m-1)$ columns of $[\mathbf{A},{\mathbf{I}}_m]$ and ${\mathbf{J}}_m$ is an $m\times m$ matrix with the $(m,m)$-entry being 1 and the other entries being 0. Since

$$\begin{array}{cc}\hfill \mathbf{W}& =\mathbf{R}/{\sigma }_m^2,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill L_0\left(a\mathbf{R}\right)& =L_0\left(\mathbf{R}\right),\forall a>0\hfill \end{array}$$

(14)

and that the locations of nonzero elements of the sparse vector ${\mathbf{r}}_{[1:k]}$ identify the true DOAs, we formulate the DOA estimation problem as solving the penalized likelihood optimization problem

$$\underset{\mathbf{r}\in \mathrm{\Omega }}{argmin}{L}_0\left(\mathbf{W}\right)+\lambda {\parallel {\mathbf{r}}_{[1:k]}\parallel }_1,$$

(15)

where

$$\mathrm{\Omega }=\{\mathbf{t}\in {\mathbb{R}}^{k+m-1}\mid \mathbf{t}\ge 0,\mathbf{B}Diag\left(\mathbf{t}\right){\mathbf{B}}^H+{\mathbf{J}}_m>0\}$$

(16)

and $\lambda \ge 0$ is pre-specified. The $l_1$-norm penalty term $\lambda \parallel {\mathbf{r}}_{[1:k]}{\parallel }_1$ would help deduce a sparse solution, and the penalty parameter $\lambda$ controls the sparsity level [37].

Remark 1.

We explain why the DOAs are not estimated by solving the plausible $l_1$-penalized ML optimization problem

$$\begin{array}{c}\hfill \underset{\mathbf{p},\mathbf{\sigma }}{argmin}{L}_0\left(\mathbf{R}\right)+\lambda {\parallel \mathbf{p}\parallel }_1.\end{array}$$

(17)

Recalling Equation (6), we find for any ${\lambda }_1>0$ and ${\lambda }_2>0$,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill L_0\left(\mathbf{R}\right)+{\lambda }_1{\parallel \mathbf{p}\parallel }_1=L_0\left([\mathbf{A},{\mathbf{I}}_m]Diag\left(\left[\frac{{\lambda }_1}{{\lambda }_2}\mathbf{p};\frac{{\lambda }_1}{{\lambda }_2}\mathbf{\sigma }\right]\right){[\mathbf{A},{\mathbf{I}}_m]}^H\right)+{\lambda }_2{∥\frac{{\lambda }_1}{{\lambda }_2}\mathbf{p}∥}_1.\end{array}\end{array}$$

(18)

Thus, $[\mathbf{p},\mathbf{\sigma }]$ is a locally optimal solution of Equation (17) with $\lambda ={\lambda }_1$ if and only if $\left[\frac{{\lambda }_1}{{\lambda }_2}\mathbf{p};\frac{{\lambda }_1}{{\lambda }_2}\mathbf{\sigma }\right]$ is a locally optimal solution of Equation (17) with $\lambda ={\lambda }_2$. This means if we estimate $\mathbf{p}$ by Equation (17), the parameter λ cannot work theoretically in adjusting the sparsity level of the estimate of $\mathbf{p}$. Instead of considering Equation (17), we formulate the optimization (Equation (15)) for the DOA estimation. In Equation (15), noticing the constant matrix ${\mathbf{J}}_m$ in Equation (12), we find that different values of λ would result in solutions with different sparsity levels.

3.2. Spatially Uniform White Noise

In the case of spatially uniform white sensor noise, ${\sigma }_1=\cdots ={\sigma }_m=\sigma$. Assuming ${\sigma }^2>0$, we denote

$$\begin{array}{cc}\hfill \mathbf{q}& =\mathbf{p}/{\sigma }^2,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \mathbf{Q}& =\mathbf{A}Diag\left(\mathbf{q}\right){\mathbf{A}}^H+{\mathbf{I}}_m.\hfill \end{array}$$

(20)

Estimating the DOAs means identifying the locations of nonzero elements of the vector $\mathbf{q}$. Considering $\mathbf{Q}=\mathbf{R}/{\sigma }^2$, we can estimate $\mathbf{q}$ through solving the penalized likelihood optimization problem

$$\underset{\mathbf{q}\ge 0}{argmin}{L}_0\left(\mathbf{Q}\right)+\lambda {\parallel \mathbf{q}\parallel }_1$$

(21)

with $\lambda \ge 0$, where the term ${\lambda \parallel \mathbf{q}\parallel }_1$ plays the same role as the penalty term in Equation (15).

Note that the number of unknown parameters to be estimated is k in the case of spatially uniform noise in contrast to $k+m-1$ in the case of spatially non-uniform noise.

4. DOA Estimation Algorithms

In this section, we provide methods to solve the optimization problems in Equations (15) and (21) with $\lambda$ fixed. As Equations (15) and (21) are non-convex, it is generally hard to give their globally optimal solutions. Based on the MM framework [39,40], we develop algorithms to find the locally optimal solutions of Equations (15) and (21).

A function $f\left(x\right)$ is said to be majorized by a function $g\left(x\right|x_0)$ at $x_0$, if $f\left(x\right)\le g\left(x\right|x_0)$ for all x and $f\left(x_0\right)=g\left(x_0\right|x_0)$. In the MM framework, the problem $arg{min}_{x}f\left(x\right)$ can be solved through iteratively solving $x_{u+1}=arg{min}_{x}g\left(x\right|x_u)$.

When solving the problems in Equations (15) and (21) in the MM framework, majorizing functions can be constructed based on the following two inequalities. For any positive definite matrices $\mathbf{U}\in {\mathbb{C}}^{m\times m}$ and ${\mathbf{U}}_u\in {\mathbb{C}}^{m\times m}$ [40],

$$logdet\left(\mathbf{U}\right)\le logdet\left({\mathbf{U}}_u\right)+tr({\mathbf{U}}_u^{-1}\mathbf{U})-m$$

(22)

and

$$log(\mathbf{x}{\left(t\right)}^H{\mathbf{U}}^{-1}\mathbf{x}\left(t\right))\le log(\mathbf{x}{\left(t\right)}^H{\mathbf{U}}_u^{-1}\mathbf{x}\left(t\right))+\frac{\mathbf{x}{\left(t\right)}^H{\mathbf{U}}^{-1}\mathbf{x}\left(t\right)}{\mathbf{x}{\left(t\right)}^H{\mathbf{U}}_u^{-1}\mathbf{x}\left(t\right)}-1,$$

(23)

where both equalities are achieved at $\mathbf{U}={\mathbf{U}}_u$.

4.1. Algorithms for Spatially Non-Uniform White Noise

In this subsection, using two different majorizing functions of the objection function in Equation (15), we develop two different MM algorithms named MM1 and MM2 to solve Equation (15).

4.1.1. MM1 Algorithm

Denote

$${\mathbf{W}}_u=\mathbf{B}Diag\left({\mathbf{r}}_u\right){\mathbf{B}}^H+{\mathbf{J}}_m,$$

(24)

where ${\mathbf{r}}_u\in \mathrm{\Omega }$. Replacing the $\mathbf{U}$ and ${\mathbf{U}}_u$ in Equations (22) and (23) by $\mathbf{W}$ and ${\mathbf{W}}_u$, respectively, we have for any $\mathbf{W}>0,{\mathbf{W}}_u>0$,

$$\begin{array}{cc}\hfill L_0\left(\mathbf{W}\right)& \le tr\left({\mathbf{W}}_u^{-1}\mathbf{W}\right)+\frac{m}{n}\sum _{t=1}^n\frac{\mathbf{x}{\left(t\right)}^H{\mathbf{W}}^{-1}\mathbf{x}\left(t\right)}{\mathbf{x}{\left(t\right)}^H{\mathbf{W}}_u^{-1}\mathbf{x}\left(t\right)}+c_2\hfill \end{array}$$

(25)

$$\begin{array}{c}=tr\left({\mathbf{W}}_u^{-1}\mathbf{W}\right)+tr\left({\mathbf{M}}_u{\mathbf{W}}^{-1}\right)+c_2\hfill \end{array}$$

(26)

where $c_2$ is a constant, the equality in Equation (25) is achieved at $\mathbf{W}={\mathbf{W}}_u$, and

$${\mathbf{M}}_u=\frac{m}{n}\sum _{t=1}^n\frac{\mathbf{x}\left(t\right)\mathbf{x}{\left(t\right)}^H}{\mathbf{x}{\left(t\right)}^H{\mathbf{W}}_u^{-1}\mathbf{x}\left(t\right)}.$$

(27)

Denote the sum of the first two terms on the right side of Equation (26) and the $l_1$-norm penalty term in Equation (15) as $g_1(\mathbf{r}\mid {\mathbf{r}}_u,\lambda )$, and then due to Equation (12), $g_1(\mathbf{r}\mid {\mathbf{r}}_u,\lambda )$ can be rewritten as

$$g_1(\mathbf{r}\mid {\mathbf{r}}_u,\lambda )=tr\left({\mathbf{M}}_u{\mathbf{W}}^{-1}\right)+{\mathbf{w}}_u^T\mathbf{r}$$

(28)

with

$$\begin{array}{cc}\hfill {\mathbf{w}}_u& =diag\left({\mathbf{B}}^H{\mathbf{W}}_u^{-1}\mathbf{B}\right)+[\lambda {\mathbf{1}}_k;{\mathbf{0}}_{m-1}].\hfill \end{array}$$

(29)

From Equation (25), $g_1(\mathbf{r}\mid {\mathbf{r}}_u,\lambda )+c_2$ is found to be a majorizing convex function of the objective function in Equation (15). Based on the MM framework, the problem in Equation (15) can be solved by iteratively solving the convex optimization problem

$${\mathbf{r}}_{u+1}=\underset{\mathbf{r}\in \mathrm{\Omega }}{argmin}{g}_1(\mathbf{r}\mid {\mathbf{r}}_u,\lambda ).$$

(30)

Proposition 1 

(The MM1 algorithm for Equation (15)). The sequence $\left\{{\mathbf{r}}_u\right\}$ generated by

$${\mathbf{r}}_{u+1}=\underset{\mathbf{r}\ge 0}{argmin}{g}_1(\mathbf{r}\mid {\mathbf{r}}_u,\lambda )$$

(31)

with any initial value ${\mathbf{r}}_0\in \mathrm{\Omega }$ converges to a locally optimal solution of the optimization problem in Equation (15).

Proof. 

Since $g_1(\mathbf{r}\mid {\mathbf{r}}_u,\lambda )$ tends to $+\infty$ as $\mathbf{W}$ goes to the boundary of the positive semidefinite cone, the optimization problems in Equations (30) and (31) are equivalent.

This proposition follows by the convergence property of the MM algorithm [40] and the fact that $L_0\left(\mathbf{W}\right)\to +\infty$ with probability 1 as $\mathbf{W}$ tends to the boundary of the positive semidefinite cone [41].  □

Due to Equation (11), ${\mathbf{r}}_0={[{\left(r_1\right)}_0,\dots ,{\left(r_{k+m-1}\right)}_0]}^T$ required in the iteration in Equation (31) can be the one obtained by inserting the initial estimate of $[\mathbf{p};\mathbf{\sigma }]$ presented in [13]. Specifically,

$$\begin{array}{c}\hfill {\left(r_i\right)}_0=\frac{{\mathbf{b}}_i^H\widehat{\mathbf{R}}{\mathbf{b}}_i}{{\widehat{\mathbf{R}}}_{m,m}{\parallel {\mathbf{b}}_i\parallel }_2^4}>0,i=1,\dots ,k+m-1,\end{array}$$

(32)

where ${\mathbf{b}}_i$ is the ith column of the matrix $\mathbf{B}$ and

$$\begin{array}{c}\hfill \widehat{\mathbf{R}}=\frac{1}{n}\sum _{t=1}^n\mathbf{x}\left(t\right)\mathbf{x}{\left(t\right)}^H.\end{array}$$

(33)

To solve the optimization problem in Equation (31) globally, we develop two available solvers: the coordinate descent (CD) and SPICE-like solvers. Denote $\mathbf{r}={[r_1,\cdots ,r_{k+m-1}]}^T$ in the following.

Proposition 2 

(The CD solver for Equation (31)). The sequence of $\mathbf{r}$ generated by iterating

$$r_i=\left(\left(\sqrt{\frac{{\beta }_i}{{\alpha }_i}}-1\right)\frac{1}{{\gamma }_i}\right)\delta ({\beta }_i-{\alpha }_i),i=1,\cdots ,k+m-1$$

(34)

with any initial value $\mathbf{r}\ge 0$ converges to the globally optimal solution ${\mathbf{r}}_{u+1}$ of the problem in Equation (31), where $\delta (·)$ is the indicator function of the interval $(0,\infty )$,

$$\begin{array}{cc}\hfill {\alpha }_i& ={\left(w_i\right)}_u,\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\beta }_i& ={\mathbf{b}}_i^H{\mathbf{N}}_i^{-1}{\mathbf{M}}_u{\mathbf{N}}_i^{-1}{\mathbf{b}}_i,\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill {\gamma }_i& ={\mathbf{b}}_i^H{\mathbf{N}}_i^{-1}{\mathbf{b}}_i,\hfill \end{array}$$

(37)

with ${\left(w_i\right)}_u$ being the ith element of the vector ${\mathbf{w}}_u$, and 

$${\mathbf{N}}_i=\sum _{l\ne i}{r}_l{\mathbf{b}}_l{\mathbf{b}}_l^H+{\mathbf{J}}_m.$$

(38)

Proof. 

By the general analysis of the CD method in [42], it is easy to find that the convex problem in Equation (31) can be globally solved by iterating

$$r_i=\underset{r_i\ge 0}{argmin}{g}_1(\mathbf{r}\mid {\mathbf{r}}_u,\lambda ),i=1,\dots ,k+m-1$$

(39)

until convergence.

Due to

$$\mathbf{W}=\mathbf{B}Diag\left(\mathbf{r}\right){\mathbf{B}}^H+{\mathbf{J}}_m=r_i{\mathbf{b}}_i{\mathbf{b}}_i^H+{\mathbf{N}}_i$$

(40)

and

$${(r_i{\mathbf{b}}_i{\mathbf{b}}_i^H+{\mathbf{N}}_i)}^{-1}={\mathbf{N}}_i^{-1}-\frac{r_i{\mathbf{N}}_i^{-1}{\mathbf{b}}_i{\mathbf{b}}_i^H{\mathbf{N}}_i^{-1}}{1+r_i{\mathbf{b}}_i^H{\mathbf{N}}_i^{-1}{\mathbf{b}}_i},$$

(41)

the iteration in Equation (39) can be reformulated as

$$r_i=\underset{r_i\ge 0}{argmin}{\alpha }_i{r}_i-\frac{{\beta }_i{r}_i}{1+{\gamma }_i{r}_i},i=1\dots ,k+m-1.$$

(42)

In addition, ${\alpha }_i>0,{\beta }_i>0,{\gamma }_i>0$ for $i=1,\dots ,k+m-1$. Solving the convex optimization problems in Equation (42) by the gradient method, we conclude that the equations in Equation (42) are equivalent to those in Equation (34). □

As the solver in Proposition 2 is derived by the CD method, we call it the CD solver. The MM1 algorithm in Proposition 1 with the CD solver is specially named the MM1-CD algorithm.

For clarity, detailed steps of the MM1-CD algorithm for Equation (15) are presented in Algorithm 1.

Note that, in the CD solver, for $i=1,\dots ,k$, we have

$${\beta }_i=\frac{m}{n}\sum _{t=1}^n{\zeta }_i{\zeta }_i^H$$

(43)

where

$${\zeta }_i=\frac{{\mathbf{a}}_i^H{\mathbf{N}}_i^{-1}\mathbf{x}\left(t\right)}{\sqrt{\mathbf{x}{\left(t\right)}^H{\mathbf{W}}_u^{-1}\mathbf{x}\left(t\right)}}$$

(44)

with ${\mathbf{a}}_i$ being the ith column of the matrix $\mathbf{A}$. The ${\beta }_i$ can be interpreted as the correlation between the signal from a possible source at the grid ${\theta }_i$ and the array responses $\mathbf{x}\left(1\right),\dots ,\mathbf{x}\left(n\right)$. From the indicator function $\delta ({\beta }_i-{\alpha }_i)$ in Equation (34), we find that it is more likely to force to zeros the powers of the assumed signals that are less correlated with the array responses.

Proposition 3 

(The SPICE-like solver for Equation (31)). The sequence of $\mathbf{r}$ generated by iterating

$$\mathbf{r}=sqrt(diag(Diag\left(\mathbf{r}\right){\mathbf{B}}^H{\mathbf{W}}^{-1}{\mathbf{M}}_u{\mathbf{W}}^{-1}\mathbf{B}Diag\left(\mathbf{r}\right))\oslash {\mathbf{w}}_u)$$

(45)

with any initial value $\mathbf{r}>0$ converges to the globally optimal solution ${\mathbf{r}}_{u+1}$ of problem in Equation (31).

Algorithm 1 The MM1-CD algorithm for spatially non-uniform white noise.

1: Give ${\mathbf{r}}_0$ and $ϵ$, 2: $u=0$, 3: repeat 4:     ${\mathbf{W}}_u=\mathbf{B}Diag\left({\mathbf{r}}_u\right){\mathbf{B}}^H+{\mathbf{J}}_m$, 5:     Calculate ${\mathbf{M}}_u$ and ${\mathbf{w}}_u$, 6:     $d=0$, 7:     ${\mathbf{r}}^0={\mathbf{r}}_u$, $\mathbf{r}={\mathbf{r}}^0$, 8:     repeat 9:         for $i=1:k+m-1$ do 10:            ${\mathbf{N}}_i={\sum }_{l\ne i}{r}_l{\mathbf{b}}_l{\mathbf{b}}_l^H+{\mathbf{J}}_m$, 11:            Calculate ${\alpha }_i$, ${\beta }_i$ and ${\gamma }_i$, 12:            if ${\beta }_i\le {\alpha }_i$ then 13:                $r_i=0$, 14:            else 15:                $r_i=(\sqrt{{\beta }_i/{\alpha }_i}-1)/{\gamma }_i$, 16:            end if 17:         end for 18:         $d=d+1$, 19:         ${\mathbf{r}}^d={[r_1,\dots ,r_{k+m-1}]}^T$, 20:     until $\parallel {\mathbf{r}}^d-{\mathbf{r}}^{d-1}{\parallel }_2/{\parallel {\mathbf{r}}^{d-1}\parallel }_2<ϵ$, 21:     $u=u+1$, 22:     ${\mathbf{r}}_u={\mathbf{r}}^d$, 23: until$\parallel {\mathbf{r}}_u-{\mathbf{r}}_{u-1}{\parallel }_2/{\parallel {\mathbf{r}}_{u-1}\parallel }_2<ϵ$.

Proof. 

It is obvious that Equation (31) and the SPICE criterion in [15] are with similar forms. By the same way as the SPICE criterion is analyzed in [15], we find that the globally optimal solution of the problem in Equation (31) is identical to the minimizer $\mathbf{r}$ of the optimization problem

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{\mathbf{r}\ge 0,\mathbf{E}\in {\mathbb{C}}^{(k+m)\times n}}{min}& tr({\mathbf{E}}^H{(Diag\left([\mathbf{r};1]\right))}^{-1}\mathbf{E})+{\mathbf{w}}_u^T\mathbf{r}\hfill \\ \hfill s.t.& [\mathbf{A},{\mathbf{I}}_m]\mathbf{E}={\mathbf{X}}_u,\hfill \end{array}\end{array}$$

(46)

where ${\mathbf{X}}_u=[\tilde{\mathbf{x}}\left(1\right),\dots ,\tilde{\mathbf{x}}\left(n\right)]$ with

$$\tilde{\mathbf{x}}\left(t\right)=\sqrt{\frac{m}{n}}\frac{\mathbf{x}\left(t\right)}{\sqrt{\mathbf{x}{\left(t\right)}^T{\mathbf{W}}_u^{-1}\mathbf{x}\left(t\right)}},t=1,\dots ,n.$$

(47)

For a fixed $\mathbf{r}$, the matrix $\mathbf{E}$ minimizing Equation (46) can be verified to be [15]

$$\begin{array}{c}\hfill \mathbf{E}=Diag\left([\mathbf{r};1]\right){[\mathbf{A},{\mathbf{I}}_m]}^H{\mathbf{W}}^{-1}{\mathbf{X}}_u,\end{array}$$

(48)

and for a fixed $\mathbf{E}$, the vector $\mathbf{r}$ minimizing Equation (46) can be readily given by

$$\begin{array}{c}\hfill \mathbf{r}=sqrt({(diag(\mathbf{E}{\mathbf{E}}^H))}_{[1:k+m-1]}\oslash {\mathbf{w}}_u).\end{array}$$

(49)

The sequences of $\mathbf{E}$ and $\mathbf{r}$, generated by alternately iterating Equations (48) and (49) from $\mathbf{r}>0$, converge to the globally optimal solution of the convex problem in Equation (46) [14,15].

Due to ${\mathbf{X}}_u{\mathbf{X}}_u^H={\mathbf{M}}_u$, iterating Equation (49) is just iterating Equation (45). Thus, the sequence of $\mathbf{r}$ generated by Equation (45) converges to the minimizer $\mathbf{r}$ of Equation (46). □

The solver in Proposition 3 is named the SPICE-like solver, and the MM1 algorithm in Proposition 1 with it is called the MM1-SPICE algorithm. Algorithm 2 summarily illustrates the MM1-SPICE algorithm for Equation (15).

Algorithm 2 The MM1-SPICE algorithm for spatially non-uniform white noise.

1: Give ${\mathbf{r}}_0$ and $ϵ$, 2: $u=0$, 3: repeat 4:     ${\mathbf{W}}_u=\mathbf{B}Diag\left({\mathbf{r}}_u\right){\mathbf{B}}^H+{\mathbf{J}}_m$, 5:     Calculate ${\mathbf{M}}_u$ and ${\mathbf{w}}_u$, 6:     $d=0$, 7:     ${\mathbf{r}}^0={\mathbf{r}}_u$, 8:     repeat 9:         ${\mathbf{W}}^d=\mathbf{B}Diag\left({\mathbf{r}}^d\right){\mathbf{B}}^H+{\mathbf{J}}_m$, 10:         ${\mathbf{r}}^{d+1}=sqrt(diag(Diag\left({\mathbf{r}}^d\right){\mathbf{B}}^H{\left({\mathbf{W}}^d\right)}^{-1}{\mathbf{M}}_u{\left({\mathbf{W}}^d\right)}^{-1}\mathbf{B}Diag\left({\mathbf{r}}^d\right))\oslash {\mathbf{w}}_u)$, 11:         $d=d+1$, 12:     until $\parallel {\mathbf{r}}^d-{\mathbf{r}}^{d-1}{\parallel }_2/{\parallel {\mathbf{r}}^{d-1}\parallel }_2<ϵ$, 13:     $u=u+1$, 14:     ${\mathbf{r}}_u={\mathbf{r}}^d$, 15: until$\parallel {\mathbf{r}}_u-{\mathbf{r}}_{u-1}{\parallel }_2/{\parallel {\mathbf{r}}_{u-1}\parallel }_2<ϵ$.

From Propositions 1–3, it is found that the proposed MM1 algorithm is by an inner–outer iteration loop. The MM1-CD and the MM1-SPICE algorithms are with the identical outer loop but different nested inner loops. In addition, relationship and difference between the MM1-CD and the MM1-SPICE are discussed in Section 4.3.

4.1.2. MM2 Algorithm

When ${\mathbf{r}}_u>0$, it is found from [36] that, for any $\mathbf{r}>0$,

$$\left[\begin{array}{cc}{\mathbf{C}}_u{(Diag\left([\mathbf{r};1]\right))}^{-1}{\mathbf{C}}_u^H& {\mathbf{I}}_m\\ {\mathbf{I}}_m& \mathbf{W}\end{array}\right]={\mathbf{F}}^H\mathbf{F}\ge 0,$$

(50)

where

$$\begin{array}{cc}\hfill {\mathbf{C}}_u& ={\mathbf{W}}_u^{-1}[\mathbf{A},{\mathbf{I}}_m]Diag\left([{\mathbf{r}}_u;1]\right),\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill \mathbf{F}& =[{(Diag\left([\mathbf{r};1]\right))}^{-\frac{1}{2}}{\mathbf{C}}_u^H,{(Diag\left([\mathbf{r};1]\right))}^{\frac{1}{2}}{[\mathbf{A},{\mathbf{I}}_m]}^H].\hfill \end{array}$$

(52)

From Equation (50), we have

$$\begin{array}{c}\hfill {\mathbf{W}}^{-1}\le {\mathbf{C}}_u{(Diag\left([\mathbf{r};1]\right))}^{-1}{\mathbf{C}}_u^H\end{array}$$

(53)

with the equality achieved at $\mathbf{r}={\mathbf{r}}_u$. The inequality in Equation (53) is also valid for any $\mathbf{r}\in \mathrm{\Omega }$ due to $\mathbf{W}>0$.

Denote

$$g_2(\mathbf{r}\mid {\mathbf{r}}_u,\lambda )=tr({\mathbf{C}}_u^H{\mathbf{M}}_u{\mathbf{C}}_u{(Diag\left([\mathbf{r};1]\right))}^{-1})+{\mathbf{w}}_u^T\mathbf{r}.$$

(54)

It is clear from Equation (53) that, when ${\mathbf{r}}_u>0$, for any $\mathbf{r}\in \mathrm{\Omega }$,

$$g_1(\mathbf{r}\mid {\mathbf{r}}_u,\lambda )\le g_2(\mathbf{r}\mid {\mathbf{r}}_u,\lambda ),$$

(55)

where the equality is achieved at $\mathbf{r}={\mathbf{r}}_u$. Therefore, at any ${\mathbf{r}}_u>0$, $g_2(\mathbf{r}\mid {\mathbf{r}}_u,\lambda )+c_2$ majorizes the objective function in Equation (15) for $\mathbf{r}\in \mathrm{\Omega }$.

Proposition 4 

(The MM2 algorithm for Equation (15)). The sequence $\left\{{\mathbf{r}}_u\right\}$ generated by

$${\mathbf{r}}_{u+1}=sqrt(diag(Diag({\mathbf{r}}_u){\mathbf{B}}^H{\mathbf{W}}_u^{-1}{\mathbf{M}}_u{\mathbf{W}}_u^{-1}\mathbf{B}Diag({\mathbf{r}}_u))\oslash {\mathbf{w}}_u)$$

(56)

with any initial value ${\mathbf{r}}_0>0$ converges to a locally optimal solution of the problem in Equation (15).

Proof. 

Through the convergence analysis in [36], we have that, although  Equation (55) is valid only when ${\mathbf{r}}_u>0$, the sequence $\left\{{\mathbf{r}}_u\right\}$ generated

$$\begin{array}{c}\hfill {\mathbf{r}}_{u+1}=\underset{\mathbf{r}\in \mathrm{\Omega }}{argmin}{g}_2(\mathbf{r}\mid {\mathbf{r}}_u,\lambda )\end{array}$$

(57)

with any initial value ${\mathbf{r}}_0>0$ converges to a locally optimal solution of the problem in Equation (15). It is worth mentioning that the elements of the coefficient vector $diag\left({\mathbf{C}}_u^H{\mathbf{M}}_u{\mathbf{C}}_u\right)$ in Equation (54) are positive. By solving the optimization problem in Equation (57) using the gradient method, the iteration procedure in Equation (57) is found to be exactly Equation (56).  □

The ${\mathbf{r}}_0$ involved in the iteration procedure in Equation (56) can be the one given in Equation (32). Summarily, Algorithm 3 gives the detailed steps of the MM2 algorithm for the problem in Equation (15).

Algorithm 3 The MM2 algorithm for spatially non-uniform white noise.

1: Give ${\mathbf{r}}_0$ and $ϵ$, 2: $u=0$, 3: repeat 4:     ${\mathbf{W}}_u=\mathbf{B}Diag\left({\mathbf{r}}_u\right){\mathbf{B}}^H+{\mathbf{J}}_m$, 5:     Calculate ${\mathbf{M}}_u$ and ${\mathbf{w}}_u$, 6:     ${\mathbf{r}}_{u+1}=sqrt(diag(Diag\left({\mathbf{r}}_u\right){\mathbf{B}}^H{\mathbf{W}}_u^{-1}{\mathbf{M}}_u{\mathbf{W}}_u^{-1}\mathbf{B}Diag\left({\mathbf{r}}_u\right))\oslash {\mathbf{w}}_u)$, 7:     $u=u+1$, 8: until$\parallel {\mathbf{r}}_u-{\mathbf{r}}_{u-1}{\parallel }_2/{\parallel {\mathbf{r}}_{u-1}\parallel }_2<ϵ$.

4.2. DOA Estimation for Spatially Uniform White Noise

The DOA estimation in the case of spatially uniform white noise amounts to solving the optimization problem in Equation (21). Problems in Equations (21) and (15) are with similar forms, but Equation (21) involves a smaller number of unknown parameters. By the same way as the problem in Equation (15) is analyzed in Section 4.1, we can naturally derive the algorithms to solve Equation (21). The algorithms for Equation (21) are still named MM1 (including MM1-CD and MM1-SPICE) and MM2.

4.2.1. MM1 Algorithm

Proposition 5 

(The MM1 algorithm for Equation (21)). Denote

$$\begin{array}{c}\hfill h_1(\mathbf{q}\mid {\mathbf{q}}_u,\lambda )=tr({\mathbf{Z}}_u{\mathbf{Q}}^{-1})+{\mathbf{e}}_u^T\mathbf{q}\end{array}$$

(58)

with

$$\begin{array}{cc}\hfill {\mathbf{Z}}_u& =\frac{m}{n}\sum _{t=1}^n\frac{\mathbf{x}\left(t\right)\mathbf{x}{\left(t\right)}^H}{\mathbf{x}{\left(t\right)}^H{\mathbf{Q}}_u^{-1}\mathbf{x}\left(t\right)},\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill {\mathbf{e}}_u& =diag({\mathbf{A}}^H{\mathbf{Q}}_u^{-1}\mathbf{A})+\lambda {\mathbf{1}}_k,\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill {\mathbf{Q}}_u& =\mathbf{A}Diag\left({\mathbf{q}}_u\right){\mathbf{A}}^H+{\mathbf{I}}_m,\hfill \end{array}$$

(61)

and then the sequence $\left\{{\mathbf{q}}_u\right\}$ generated by

$$\begin{array}{c}\hfill {\mathbf{q}}_{u+1}=\underset{\mathbf{q}\ge 0}{argmin}{h}_1(\mathbf{q}\mid {\mathbf{q}}_u,\lambda )\end{array}$$

(62)

with any initial value $\mathbf{q}\ge 0$ converges to a locally optimal solution of the problem in Equation (21).

To solve the optimization problem in Equation (62), using the same method as Equation (31) is solved, we offer two different iterative solvers for Equation (62).

Proposition 6 

(The CD solver for Equation (62)). The sequence of $\mathbf{q}$ generated by iterating

$$\begin{array}{c}\hfill q_j=\left(\frac{{\tilde{\beta }}_j}{{\tilde{\alpha }}_j}-1\right)\frac{1}{{\tilde{\gamma }}_j}\delta ({\tilde{\beta }}_j-{\tilde{\alpha }}_j),j=1,\dots ,k\end{array}$$

(63)

with any initial value $\mathbf{q}\ge 0$ converges to the globally optimal solution of the problem in Equation (62), where

$$\begin{array}{cc}\hfill {\tilde{\alpha }}_j& ={\left(e_j\right)}_u,\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill {\tilde{\beta }}_j& ={\mathbf{a}}_j^H{\left({\mathbf{H}}_j\right)}^{-1}{\mathbf{Z}}_u{\left({\mathbf{H}}_j\right)}^{-1}{\mathbf{a}}_j,\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill {\tilde{\gamma }}_j& ={\mathbf{a}}_j^H{\left({\mathbf{H}}_j\right)}^{-1}{\mathbf{a}}_j,\hfill \end{array}$$

(66)

with ${\left(e_j\right)}_u$ being the jth element of the vector ${\mathbf{e}}_u$, ${\mathbf{a}}_j$ being the jth column of the matrix $\mathbf{A}$ and

$$\begin{array}{c}\hfill {\mathbf{H}}_j=\sum _{l\ne j}{q}_l{\mathbf{a}}_l{\mathbf{a}}_l^H+{\mathbf{I}}_m.\end{array}$$

(67)

Proposition 6 introduces the nested CD inner loop of the MM1 algorithm for the problem in Equation (21), and can be proven similar to Proposition 2.

Proposition 7 

(The SPICE-like solver for Equation (62)). The sequence of $\mathbf{q}$ generated by iterating

$$\begin{array}{cc}\hfill \mathbf{q}& =sqrt(diag(Diag(\mathbf{q}){\mathbf{A}}^H{\mathbf{Q}}^{-1}{\mathbf{Z}}_u{\mathbf{Q}}^{-1}\mathbf{A}Diag(\mathbf{q}))\oslash {\mathbf{e}}_u)\hfill \end{array}$$

(68)

with any initial value $\mathbf{q}>0$ converges to the globally optimal solution of the problem in Equation (62).

The iterative procedure in Proposition 7 above is the nested SPICE-like inner loop of the MM1 algorithm for Equation (21), and can be proven similar to Proposition 3.

The MM procedure in Proposition 5 with the CD nested loop in Proposition 6 is the MM1-CD algorithm for Equation (21). The MM1-SPICE algorithm for Equation (21) is the MM procedure in Proposition 5 with the SPICE-like nested loop in Proposition 7.

4.2.2. MM2 Algorithm

Proposition 8 

(The MM2 algorithm for Equation (21)). The sequence $\left\{{\mathbf{q}}_u\right\}$ generated by

$${\mathbf{q}}_{u+1}=sqrt(diag(Diag({\mathbf{q}}_u){\mathbf{A}}^H{\mathbf{Q}}_u^{-1}{\mathbf{Z}}_u{\mathbf{Q}}_u^{-1}\mathbf{A}Diag({\mathbf{q}}_u))\oslash {\mathbf{e}}_u)$$

(69)

with any initial value $\mathbf{q}>0$ converges to a locally optimal solution of the problem in Equation (21).

Due to Equation (19), by means of inserting in $\mathbf{q}$ the initial estimates of $\mathbf{p}$ and ${\sigma }^2$ given in [13], ${\mathbf{q}}_0={[{\left(q_1\right)}_0,\dots ,{\left(q_k\right)}_0]}^T$ required by the iterations in Equations (62) and (69) can be

$${\left(q_j\right)}_0=\frac{{\mathbf{a}}_i^H\widehat{\mathbf{R}}{\mathbf{a}}_i}{\widehat{\sigma }{\parallel {\mathbf{a}}_i\parallel }_2^4}>0,j=1,\dots ,k,$$

(70)

where $\widehat{\mathbf{R}}$ is given by Equation (33) and $\widehat{\sigma }$ is the mean of the first m smallest values of the set

$$\left\{\frac{{\mathbf{a}}_1^H\widehat{\mathbf{R}}{\mathbf{a}}_1}{\parallel {\mathbf{a}}_1{\parallel }_2^2},\dots ,\frac{{\mathbf{a}}_k^H\widehat{\mathbf{R}}{\mathbf{a}}_k}{\parallel {\mathbf{a}}_k{\parallel }_2^2},{\widehat{\mathbf{R}}}_{1,1},\dots ,{\widehat{\mathbf{R}}}_{m,m}\right\}.$$

(71)

4.3. Discussions on the MM1 and MM2 Algorithms

The following arguments are focused on the case of spatially non-uniform noise, which are also applicable to the case of spatially uniform noise.

Both the MM1 and MM2 algorithms decrease the objective function in Equation (15) at each MM iteration and converge locally, in which the different majorizing functions are adopted. Compared to the MM2 algorithm, the MM1 algorithm has a majorizing function closer to the objective function in Equation (15) due to Equation (55). The computational burdens of the two algorithms are mainly caused by the matrix inversion operations.

Although the MM1-CD and MM1-SPICE algorithms have different nested inner iteration procedures, they converge to the same local solution theoretically because their outer MM iteration procedures are both Equation (31). Each nested inner iteration of the MM1-CD algorithm, detailed by Steps 9–17 in Algorithm 1, requires $k+m-1$ matrix inverse operations. In each nested inner iteration of the MM1-SPICE algorithm, presented by Steps 9–10 in Algorithm 2, only one matrix inverse operation is entailed.

It is somewhat interesting to find that the iteration of the MM2 algorithm (see Equation (56)) and the nested inner iteration of the MM1-SPICE algorithm (see Equation (45)) have similar forms. In each iteration of the MM2 algorithm, only one matrix inverse operation is needed.

4.4. Selection of Penalty Parameter

The penalty parameter $\lambda$ in Equations (15) and (21) affects the sparsity levels of the estimates of $\mathbf{r}$ and $\mathbf{q}$. A modified Bayesian information criterion (BIC) method [37], which is common and statistically sound, is provided here to choose an appropriate $\lambda$. Let ${\widehat{\mathbf{r}}}_{\lambda }$ and ${\widehat{\mathbf{q}}}_{\lambda }$ be the solutions of the problems in Equations (15) and (21) with a fixed $\lambda$, respectively, and denote ${\widehat{\mathbf{W}}}_{\lambda }=\mathbf{B}Diag\left({\widehat{\mathbf{r}}}_{\lambda }\right){\mathbf{B}}^H+{\mathbf{J}}_m$ and ${\widehat{\mathbf{Q}}}_{\lambda }=\mathbf{A}Diag\left({\widehat{\mathbf{q}}}_{\lambda }\right){\mathbf{A}}^H+{\mathbf{I}}_m$. The appropriate $\lambda$ for spatially non-uniform noise and uniform noise are

$$\underset{\lambda }{argmin}{L}_0\left({\widehat{\mathbf{W}}}_{\lambda }\right)+{\widehat{z}}_{\lambda ,1}log\left(n\right)/\left(2n\right)$$

(72)

and

$$\underset{\lambda }{argmin}{L}_0({\widehat{\mathbf{Q}}}_{\lambda })+{\widehat{z}}_{\lambda ,2}log\left(n\right)/\left(2n\right),$$

(73)

respectively, where ${\widehat{z}}_{\lambda ,1}$ and ${\widehat{z}}_{\lambda ,2}$ are the numbers of nonzero elements of the vectors ${\left({\widehat{\mathbf{r}}}_{\lambda }\right)}_{[1:k]}$ and ${\widehat{\mathbf{q}}}_{\lambda }$, respectively.

Note that the elements of ${\left({\widehat{\mathbf{r}}}_{\lambda }\right)}_{[1:k]}$ and ${\widehat{\mathbf{q}}}_{\lambda }$ can be treated as 0 if they are smaller than some certain values respectively, e.g., $\epsilon max\left({\left({\widehat{\mathbf{r}}}_{\lambda }\right)}_{[1:k]}\right)$ and $\epsilon max\left({\widehat{\mathbf{q}}}_{\lambda }\right)$ with a very small $\epsilon >0$.

Notice that the $L_0({\widehat{\mathbf{W}}}_{\lambda })$ and $L_0({\widehat{\mathbf{Q}}}_{\lambda })$ in Equations (72) and (73) are substitutes for the marginal likelihood also called as Bayesian evidence required in the BIC criterion [43]. By the way, when modeling in the Bayesian framework, the marginal likelihood usually cannot be analytically calculated, but can be approximated by several computational methods in [44,45,46].

5. Numerical Experiments

In this section, we numerically show the performance of the proposed methods. Consider a uniform linear array (ULA), and, for each $j=1,\dots ,k_0$, the steering vector corresponding to the DOA ${\xi }_j$ is

$$\begin{array}{c}\hfill \mathbf{a}\left({\xi }_j\right)={[1,exp(ı\pi cos\left({\xi }_j\right)),\dots ,exp(ı\pi (m-1)cos\left({\xi }_j\right))]}^T.\end{array}$$

(74)

The array output data $\mathbf{x}\left(t\right)$ are generated from the model

$$\begin{array}{c}\hfill \mathbf{x}\left(t\right)=\sum _{j=1}^{k_0}\mathbf{a}\left({\xi }_j\right){\overline{s}}_j\left(t\right)+\mathbf{v}\left(t\right),t=1,\dots ,n,\end{array}$$

(75)

and both the source signals ${\overline{s}}_1\left(t\right),\dots ,{\overline{s}}_{k_0}\left(t\right)$ and the observation noise $\mathbf{v}\left(t\right)$ are temporally independent. The SNR is defined as

$$\mathrm{SNR}\left(\mathrm{dB}\right)=10{log}_{10}\left(\frac{tr(\overline{\mathbf{A}}\overline{\mathbf{P}}{\overline{\mathbf{A}}}^H)}{tr\left(\mathbf{V}\right)}\right),$$

(76)

where $\overline{\mathbf{A}}=[\mathbf{a}\left({\xi }_1\right),\dots ,\mathbf{a}\left({\xi }_{k_0}\right)]$, $\overline{\mathbf{P}}$ is the covariance matrix of $\overline{\mathbf{s}}\left(t\right)={[{\overline{s}}_1\left(t\right),\dots ,{\overline{s}}_{k_0}\left(t\right)]}^T$, and $\mathbf{V}$ is the covariance matrix of $\mathbf{v}\left(t\right)$. The root mean-square error (RMSE) of the DOA estimate is employed to evaluate the estimation performance, which is approximated by $R=500$ Monte Carlo runs as

$$RMSE=\frac{1}{R}\sum _{i=1}^R\sqrt{\frac{1}{k_0}\sum _{j=1}^{k_0}{({\widehat{\xi }}_{j,i}-{\xi }_j)}^2},$$

(77)

where ${\widehat{\xi }}_{j,i}$ is the estimate of ${\xi }_j$ in the ith Monte Carlo run. In the following experiments, we applied the proposed methods and the SPICE [13,15] and LIKES methods [15] to estimate the DOAs. Set the grid points $\mathbf{\theta }=\{0^{\circ },1^{\circ },\dots ,{180}^{\circ }\}$, and the tolerance $ϵ={10}^{-8}$ for convergence. All methods were coded in MATLAB and executed on a workstation with two 2.10 GHz Intel Xeon CPUs.

5.1. Experimental Settings

5.1.1. Experiment 1

Consider a scenario with Gaussian source signals and noise:

• The number of sensors in the array is $m=15$, and the number of snapshots is $n=100$.
• There are $k_0=4$ sources locating at ${\xi }_1={60}^{\circ },{\xi }_2={80}^{\circ },{\xi }_3={100}^{\circ },{\xi }_4={120}^{\circ }$, respectively.
• The signal $\overline{\mathbf{s}}\left(t\right)$ is zero-mean circular complex Gaussian with covariance matrix $\overline{\mathbf{P}}=Diag\left([100,100,100,100]\right)$.
• The noise $\mathbf{v}\left(t\right)$ is zero-mean circular complex Gaussian with covariance matrix $\mathbf{V}={\sigma }^2{\mathbf{I}}_m$.

5.1.2. Experiment 2

Consider the scenario where both the source signals and the observation noise are non-Gaussian. Let $m=15$. We used the experimental settings:

• The $k_0=4$ source signals from ${\xi }_1={10}^{\circ }, {\xi }_2={40}^{\circ }, {\xi }_3={40}^{\circ }+\Delta \xi , {\xi }_4={70}^{\circ }$ are, respectively, as ${\overline{s}}_1\left(t\right)=5exp\left(ı{\phi }_1\left(t\right)\right), {\overline{s}}_2\left(t\right)=10exp\left(ı{\phi }_2\left(t\right)\right), {\overline{s}}_3\left(t\right)=10exp\left(ı{\phi }_3\left(t\right)\right), {\overline{s}}_4\left(t\right)=5exp\left(ı{\phi }_4\left(t\right)\right)$, where $\Delta \xi$ is the angle separation between ${\xi }_2$ and ${\xi }_3$, and ${\phi }_1\left(t\right),\dots ,{\phi }_4\left(t\right)$ are independent and respectively, distributed as

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\phi }_1\left(t\right)& \sim 2\pi \mathrm{Gamma}(1,1),\hfill \\ \hfill {\phi }_2\left(t\right)& \sim 2\pi \mathrm{Gamma}(1,1),\hfill \\ \hfill {\phi }_3\left(t\right)& \sim \mathrm{U}[0,2\pi ],\hfill \\ \hfill {\phi }_4\left(t\right)& \sim \mathrm{U}[0,2\pi ],\hfill \end{array}\end{array}$$

  (78)

  where $\mathrm{Gamma}(a,b)$ denotes the Gamma distribution with the shape parameter a and scale parameter b, and $\mathrm{U}[a,b]$ denotes the uniform distribution on the interval $[a,b]$. Note that the source signals ${\overline{s}}_1\left(t\right)$, ${\overline{s}}_2\left(t\right)$, ${\overline{s}}_3\left(t\right)$ and ${\overline{s}}_4\left(t\right)$ have constant modulus, which is common in communication applications [13].
• The observation noise $\mathbf{v}\left(t\right)$ is distributed as

  $$\begin{array}{c}\hfill \mathbf{v}\left(t\right)\sim \varrho {\mathbf{n}}_1,\end{array}$$

  (79)

  where $\varrho \sim \mathrm{Gamma}(1,6)$ and ${\mathbf{n}}_1$ is zero-mean Gaussian with covariance matrix $\mathbf{V}={\sigma }^2{\mathbf{I}}_m$.

5.1.3. Experiment 3

Consider a scenario with non-Gaussian source signals and non-Gaussian spatially non-uniform white noise. Let $m=8$ and $n=100$. The $k_0=3$ source signals locating at ${\xi }_1={40}^{\circ },{\xi }_2={50}^{\circ },{\xi }_3={80}^{\circ }$ are ${\overline{s}}_1\left(t\right)=10exp\left(ı{\phi }_2\left(t\right)\right)$, ${\overline{s}}_2\left(t\right)=30exp\left(ı{\phi }_3\left(t\right)\right)$ and ${\overline{s}}_3\left(t\right)=10exp\left(ı{\phi }_4\left(t\right)\right)$ with ${\phi }_2\left(t\right),{\phi }_3\left(t\right),{\phi }_4\left(t\right)$ being given by Equation (78). The observation noise is distributed as

$$\begin{array}{c}\hfill \mathbf{v}\left(t\right)\sim \varrho {\mathbf{n}}_2,\end{array}$$

(80)

where $\varrho \sim \mathrm{Gamma}(2,2)$ and ${\mathbf{n}}_2$ is zero-mean Gaussian with covariance matrix $\mathbf{V}={\sigma }^{2}Diag\left([m,m-1,\dots ,1]\right)$.

5.2. Experimental Results

In Experiment 1, the MM1-CD, MM1-SPICE and MM2 algorithms with $\lambda =0$ were firstly applied.

Table 1. Comparison of computational complexities and RMSEs.

SNR	Method	${\mathbf{n}}_1$	${\mathbf{n}}_2$	${\mathbf{n}}_3$	$\mathbf{\tau }$	RMSE
$-1\mathrm{dB}$	MM1-CD	$525.5600$	$16.0440$	$32.7550$	$27.9710$	$0.0877$
	MM1-SPICE	$317.6500$	$16.9280$	$18.7370$	$0.3805$	$0.0926$
	MM2	$129.5000$	−	−	$0.3266$	$0.1059$
$-3\mathrm{dB}$	MM1-CD	$539.2700$	$16.0000$	$33.7040$	$28.7700$	$0.2345$
	MM1-SPICE	$310.0300$	$16.6260$	$18.6070$	$0.3714$	$0.2434$
	MM2	$129.3800$	−	−	$0.3260$	$0.2568$
$-5\mathrm{dB}$	MM1-CD	$561.0900$	$16.4380$	$34.1500$	$29.3910$	$0.4292$
	MM1-SPICE	$303.1400$	$16.3340$	$18.5130$	$0.3590$	$0.4681$
	MM2	$129.1400$	−	−	$0.3195$	$0.4840$

reports the iteration numbers $n_1,n_2,n_3$, the computational time $\tau$ (in seconds), and the DOA estimation RMSEs (in degrees). Specifically, $n_1$ is the number of total iterations, $n_2$ is the number of iterations in the outer loop, amd $n_3$ is the average of the iterations in a nested inner loop. Besides the RMSEs, each value in Table 1 is the average of the results of 500 Monte Carlo runs. Note that, even though $n_1=n_2\times n_3$ in a Monte Carlo run, $n_1$ is not exactly equal to $n_2\times n_3$ in Table 1 because they are the average of 500 Monte Carlo runs.

As shown in Table 1, the MM1-CD and MM1-SPICE algorithms had similar RMSEs and the MM1-CD algorithm took much more time than the MM1-SPICE algorithm. Considering that the MM1-CD and MM1-SPICE algorithms theoretically converge to the identical local solution, we believe that they will always have similar RMSEs. Moreover, the MM2 algorithm had the smallest iteration numbers, the least computational time and satisfactory RMSEs. In the following, we present the DOA estimation performance only of the MM1-SPICE and MM2 algorithms. Note that the MM1-SPICE and MM2 algorithms may converge to different local solutions, and then they are referred to as two different estimation methods hereinafter. Specifically, we compared the following six methods:

• LIKES: The ML method proposed in [15] under the assumption that the array output $\mathbf{x}\left(t\right)$ is Gaussian distributed.
• SPICE: A sparse covariance-based estimation method proposed in [13] with no distribution assumption.
• MM1-SPICE-0: The MM1-SPICE method with the penalty parameter $\lambda =0$;
• MM1-SPICE-P: The MM1-SPICE method with the penalty parameter selected by the criterion in Section 4.4.
• MM2-0: The MM2 method with the penalty parameter $\lambda =0$.
• MM2-P: The MM2 method with the penalty parameter selected by the criterion in Section 4.4.

It is worth presenting that in, Experiments 1–3, the standard MUSIC algorithm was found to be almost ineffective, thus we do not illustrate it.

The RMSE comparisons of the above six methods are illustrated in Figure 1, Figure 2 and Figure 3 for Experiments 1–3 with different SNRs, respectively. Figure 4 and Figure 5 for Experiment 2 show the DOA estimation RMSEs of the above six methods versus the number of snapshots and the angle separation, respectively.

In Figure 1, we can see that the MM1-SPICE-0 and MM2-0 methods were comparable to the LIKES and SPICE methods in the Gaussian cases. As shown in Figure 2, Figure 3, Figure 4 and Figure 5, the scenarios of non-Gaussian random source signals and heavy-tailed random noise, the MM1-SPICE-0 and MM2-0 methods performed much better than the LIKES and SPICE methods. In other words, the ML methods designed for the CES outputs were effective in the simulation scenarios of Gaussian and non-Gaussian distributions.

As shown in Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5, we found that the penalized ML methods, i.e., the MM1-SPICE-P and MM2-P methods, had smaller RMSEs than the other four methods. As shown in Figure 4, with the increase of the number of snapshots, the performance of the MM1-SPICE-P and MM2-P methods improved but the performance of the MM1-SPICE-0 and MM2-0 methods remained virtually unchanged. Figure 5 shows that, as $\Delta \xi$ increased from $2^{\circ }$ to $6^{\circ }$, the performance of the MM1-SPICE-P and MM2-P methods became better, while, when $\Delta \xi$ was larger than $4^{\circ }$, the performance of the MM1-SPICE-0 and MM2-0 methods no longer improved. As shown in Figure 1, Figure 2 and Figure 3, unsurprisingly, as the SNR increased, the performance of all the six methods became better.

For illustrating the difference between the penalized and un-penalized ML methods, we evaluated the normalized spectrum (NS):

$$\begin{array}{c}\hfill \mathrm{NS}=\frac{{\widehat{\mathbf{r}}}_{[1:k]}}{max\left({\widehat{\mathbf{r}}}_{[1:k]}\right)}\mathrm{or}\frac{\widehat{\mathbf{q}}}{max\left(\widehat{\mathbf{q}}\right)},\end{array}$$

(81)

where $\widehat{\mathbf{r}}$ and $\widehat{\mathbf{q}}$ are, respectively, the estimates of $\mathbf{r}$ and $\mathbf{q}$.

Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 report the NSs of Experiments 1–3 by the MM2-0 and MM2-P methods (the curves of the NSs of the MM-SPICE-0 and MM2-SPICE-P methods are similar to those of the MM2-0 and MM2-P methods, respectively). Figure 6 and Figure 10 are for Experiments 1 and 3 with $SNR=0\mathrm{dB}$, respectively. Figure 7 and Figure 8 are both for Experiment 2 with $n=100$ and $SNR=0\mathrm{dB}$, but they involve two scenarios with different angle separations. Particularly, Figure 8 is for the very small $\Delta \xi =3^{\circ }$, as opposed to Figure 7 for $\Delta \xi ={15}^{\circ }$. Figure 9 is also for Experiment 2, but it is about the scenario with small snapshot size $n=20$, small $SNR=-10\mathrm{dB}$ and moderate angle separation $\Delta \xi ={10}^{\circ }$. Note that the red dashed lines mark the true DOAs and the NS curves in each figure are the results of a randomly selected realization. We can clearly see in Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 that the MM2-P method having the $l_1$-norm penalty added yielded a higher angular resolution.

6. Conclusions

This paper provides a penalized ML method for DOA estimation under the assumption that the array output has a CES distribution. Two MM algorithms, named MM1 and MM2, are developed for solving the non-convex penalized ML optimization problem for spatially uniform or non-uniform noise. They converge locally with different rates. The numerical experiments showed that the MM2 ran faster and performed as well as the MM1 when estimating the DOAs. Their rates of convergence will be further explored theoretically in future.

It is worth mentioning that, in numerical simulations, the proposed $l_1$-norm penalized likelihood method effectively estimated the DOAs, although the values of nonzero elements of $\mathbf{r}$ or $\mathbf{p}$ were not estimated very accurately. By replacing the $l_1$-norm penalty by other proper penalties (e.g., the smoothly clipped absolute deviation (SCAD) [47], adaptive LASSO [48], or $l_p$-norm penalty ($0\le p<1$)), more accurate estimates of the unknown parameter may be derived in the future.

If the $l_1$-norm penalty is replaced by the adaptive LASSO penalty, then the algorithms proposed in this paper can be applied almost without modification because the adaptive LASSO penalty is a weighted $l_1$-norm penalty. When the non-convex SCAD or the $l_p$-norm penalty ($0\le p<1$) is employed, the convex majorizing functions given in [40,49] can be exploited based on the MM framework, and then the algorithms in this paper with minor modifications are also applicable.

When we have some informative prior knowledge on the directions of source signals, the problem of estimating the DOAs can be formulated as maximizing a sparse posterior likelihood from the Bayesian perspective. The sparse Bayesian method of estimating DOAs from the CES array outputs is interesting and is worth studying, while how to formulate and solve a sparse posterior likelihood optimization and how to do model selection are challenges. The BIC criterion can be used for the model selection, in which the Bayesian evidence difficult to be analytically derived can be approximated by the methods in [44,45,46]. More research along this line will be done in the future.