PUMA criterion = MODE criterion

TL;DR

This paper demonstrates that the PUMA estimator for array processing is mathematically equivalent to the well-known MODE estimator, unifying two approaches under a common criterion.

Contribution

It reveals that PUMA and MODE estimators minimize the same criterion function, establishing a theoretical connection between the two methods.

Findings

PUMA and MODE estimators are mathematically equivalent.

Both estimators minimize the same criterion function.

The equivalence simplifies understanding and application of array processing techniques.

Abstract

We show that the recently proposed (enhanced) PUMA estimator for array processing minimizes the same criterion function as the well-established MODE estimator. (PUMA = principal-singular-vector utilization for modal analysis, MODE = method of direction estimation.)

Equations

$$\mathbf{y}(t)=\mathbf{A}(\mathbf{\varphi })\mathbf{s}(t)+\mathbf{n}(t)\in {\mathbb{C}}^m$$

$$\mathbf{R}=\mathbf{AP}{\mathbf{A}}^{\ast }+{\sigma }^2{\mathbf{I}}_m,$$

$$\hat{\mathbf{\varphi }}=\underset{\mathbf{\varphi }}{\mathrm{argmin}}\mathrm{tr}\left\{{\mathbf{\Pi }}_{\mathbf{A}}^{\perp }\hat{\mathbf{R}}\right\},$$

$$\hat{\mathbf{R}}=\frac{1}{T}\sum _{t=1}^T\mathbf{y}(t){\mathbf{y}}^{\ast }(t)$$

$$\mathbf{A}=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ e^{j{\varphi }_1}& e^{j{\varphi }_2}& \cdots & e^{j{\varphi }_r}\\ ⋮& ⋮& & ⋮\\ e^{j(m-1){\varphi }_1}& e^{j(m-1){\varphi }_2}& \cdots & e^{j(m-1){\varphi }_r}\end{array}\right].$$

$$\mathbf{TA}=0$$

$$\mathbf{T}=\left[\begin{array}{cccccc}{c}_0& c_1& \cdots & c_r& & \\ 0& \ddots & \ddots & & \ddots & 0\\ & & c_0& c_1& \cdots & c_r\end{array}\right]\in {\mathbb{C}}^{(m-r)\times m}$$

$$c_0+c_{1}z+\cdots +c_r{z}^r=c_0\prod _{k=1}^r(1-e^{-j{\varphi }_k}z),c_0\ne 0.$$

$${\mathbf{\Pi }}_{\mathbf{A}}^{\perp }={\mathbf{\Pi }}_{\mathbf{T}}={\mathbf{T}}^{\ast }(\mathbf{T}{\mathbf{T}}^{\ast }{)}^{-1}\mathbf{T}$$

$$\hat{\mathbf{c}}=\underset{\mathbf{c}}{\mathrm{argmin}}{V}_{\text{ML}}(\mathbf{c}),$$

$$V_{\text{ML}}(\mathbf{c})=\mathrm{tr}\left\{{\mathbf{\Pi }}_{\mathbf{T}}\hat{\mathbf{R}}\right\}=\mathrm{tr}\left\{(\mathbf{T}{\mathbf{T}}^{\ast }{)}^{-1}\mathbf{T}\hat{\mathbf{R}}{\mathbf{T}}^{\ast }\right\}.$$

$$\mathbf{R}={\mathbf{U}}_s\mathbf{\Lambda }{\mathbf{U}}_s^{\ast }+{\sigma }^2{\mathbf{U}}_n{\mathbf{U}}_n^{\ast }$$

$${\hat{\mathbf{U}}}_s\hat{\mathbf{\Gamma }}{\hat{\mathbf{U}}}_s^{\ast },$$

$$\hat{\mathbf{\Gamma }}\triangleq \text{diag}\left(\frac{({\hat{\lambda }}_1-{\hat{\sigma }}^2{)}^2}{{\hat{\lambda }}_1},\dots ,\frac{({\hat{\lambda }}_r-{\hat{\sigma }}^2{)}^2}{{\hat{\lambda }}_r}\right)$$

$$V_{\text{MODE}}(\mathbf{c})=\mathrm{tr}\left\{(\mathbf{T}{\mathbf{T}}^{\ast }{)}^{-1}\mathbf{T}{\hat{\mathbf{U}}}_s\hat{\mathbf{\Gamma }}{\hat{\mathbf{U}}}_s^{\ast }{\mathbf{T}}^{\ast }\right\}.$$

$$V_{\text{PUMA}}(\mathbf{c})={\mathbf{e}}^{\ast }\hat{\mathbf{W}}\mathbf{e},$$

$$\begin{array}{rl}\hat{\mathbf{W}}& \triangleq (\hat{\mathbf{\Gamma }}\otimes (\mathbf{T}{\mathbf{T}}^{\ast }{)}^{-1})\end{array}$$

$$\begin{array}{rl}{V}_{\text{PUMA}}(\mathbf{c})& ={\mathbf{e}}^{\ast }\hat{\mathbf{W}}\mathbf{e}\\ & =\mathrm{vec}(\mathbf{T}{\hat{\mathbf{U}}}_s{)}^{\ast }\left(\hat{\mathbf{\Gamma }}\otimes (\mathbf{T}{\mathbf{T}}^{\ast }{)}^{-1}\right)\mathrm{vec}(\mathbf{T}{\hat{\mathbf{U}}}_s)\\ & =\mathrm{vec}(\mathbf{T}{\hat{\mathbf{U}}}_s{)}^{\ast }\mathrm{vec}((\mathbf{T}{\mathbf{T}}^{\ast }{)}^{-1}\mathbf{T}{\hat{\mathbf{U}}}_s\hat{\mathbf{\Gamma }})\\ & =\mathrm{tr}\left\{{\hat{\mathbf{U}}}_s^{\ast }{\mathbf{T}}^{\ast }(\mathbf{T}{\mathbf{T}}^{\ast }{)}^{-1}\mathbf{T}{\hat{\mathbf{U}}}_s\hat{\mathbf{\Gamma }}\right\}\\ & =\mathrm{tr}\left\{(\mathbf{T}{\mathbf{T}}^{\ast }{)}^{-1}\mathbf{T}{\hat{\mathbf{U}}}_s\hat{\mathbf{\Gamma }}{\hat{\mathbf{U}}}_s^{\ast }{\mathbf{T}}^{\ast }\right\}\\ & =V_{\text{MODE}}(\mathbf{c}),\end{array}$$

$$\begin{array}{rl}\mathrm{vec}(\mathbf{XYZ})& =({\mathbf{Z}}^{\top }\otimes \mathbf{X})\mathrm{vec}(\mathbf{Y})\\ \text{tr}\{{\mathbf{X}}^{\ast }\mathbf{Y}\}& =\mathrm{vec}(\mathbf{X}{)}^{\ast }\mathrm{vec}(\mathbf{Y}).\end{array}$$

Full text

I Problem formulation

The standard signal model in array processing is

[TABLE]

where $\boldsymbol{\phi}=[\phi_{1}\>\cdots\>\phi_{r}]^{\top}$ parameterizes the unknown directions of arrival from $r<m$ far-field sources, $\mathbf{s}(t)$ is a vector of unknown source signals, $\mathbf{n}(t)$ is a noise term, and $\mathbf{A}(\cdot)$ is a known function describing the array response [1, 2]. The covariance matrix of the received signals is

[TABLE]

where $\mathbf{P}$ and $\sigma^{2}\mathbf{I}_{m}$ are the signal and noise covariances, respectively. The data is assumed to be circular Gaussian.

Given $T$ independent snapshots $\{\mathbf{y}(t)\}^{T}_{t=1}$, the maximum likelihood (ML) estimate of $\boldsymbol{\phi}$ is given by

[TABLE]

where

[TABLE]

denotes the sample covariance matrix and $\boldsymbol{\Pi}^{\perp}_{\mathbf{A}}$ is the orthogonal projector onto $\mathcal{R}(\mathbf{A})^{\perp}$ and is a nonlinear function of $\boldsymbol{\phi}$. The nonconvex problem in (3) can be viewed as fitting the signal subspace spanned by $\mathbf{A}$ to the data, and it can be tackled using numerical search techniques.

When considering uniform linear arrays, the columns of $\mathbf{A}$ have a Vandermonde structure:

[TABLE]

In this case we have the following orthogonal relation

[TABLE]

where

[TABLE]

is a Toeplitz matrix with coefficients $\mathbf{c}=[c_{0}\;c_{1}\cdots\;c_{r}]^{\top}$. These coefficients also define a polynomial with roots that lie on the unit circle,

[TABLE]

Therefore there is a direct correspondence between $\boldsymbol{\phi}$ and $\mathbf{c}$ [1, 2]. As a consequence of (4) the orthogonal projector can be written as

[TABLE]

which yields an equivalent problem to (3) in terms of $\mathbf{c}$:

[TABLE]

where

[TABLE]

Using this alternative parameterization, tractable minimization algorithms can be formulated. Next, we consider two alternative estimation criteria and prove that they are equivalent.

II PUMA criterion equals MODE criterion

Using the eigendecomposition, the covariance matrix can be written as

[TABLE]

where $\mathcal{R}(\mathbf{U}_{s})=\mathcal{R}(\mathbf{A})$ and $\boldsymbol{\Lambda}=\text{diag}(\lambda_{1},\dots,\lambda_{r})\succ\mathbf{0}$ is the matrix of eigenvalues that are larger than $\sigma^{2}$. Instead of fitting the subspace to the sample covariance $\widehat{\mathbf{R}}$, as in (6), consider fitting to a weighted estimate of the signal subspace [3, 4]:

[TABLE]

where

[TABLE]

and where $\{\hat{\lambda}_{i}\}$ and $\hat{\sigma}^{2}$ are obtained from the eigendecomposition of $\widehat{\mathbf{R}}$. Then the cost function in (5) is replaced by

[TABLE]

This leads to the asymptotically efficient ‘method of direction estimation’ (Mode) [3][2, ch. 8.5]. A simple two-step algorithm was proposed in [3] to approximate the minimum of the above estimation criterion.

Another approach for array processing, called ‘principal-singular-vector utilization for modal analysis’ (Puma), has been recently proposed in [5] (see also references therein for predecessors of that approach). It is motivated by properties of a related linear prediction problem and based on the following fitting criterion

[TABLE]

where

[TABLE]

is a weighting matrix and $\mathbf{e}$ is a function of $\mathbf{c}$ and the eigenvectors in $\widehat{\mathbf{U}}_{s}$. As shown in [5], $\mathbf{e}$ can be written as $\mathbf{e}=\text{vec}(\mathbf{T}\widehat{\mathbf{U}}_{s})$. It follows immediately that

[TABLE]

where we made use of the following results

[TABLE]

Therefore the Puma criterion is exactly equivalent to the Mode criterion. The algorithm proposed in [5] is thus an alternative technique for minimizing $V_{\text{MODE}}(\mathbf{c})$.

III Other variants

A fitting criterion on a similar form as $V_{\text{PUMA}}(\mathbf{c})$ was proposed in [6] and shown to reduce to $V_{\text{MODE}}(\mathbf{c})$ in a special case. Alternative minimization techniques are also discussed therein, see also [2, ch. 8]. See e.g. [7, 8] for additional variations of $V_{\text{MODE}}(\mathbf{c})$.

In scenarios with low signal-to-noise ratio or small sample size $T$, subspace-fitting methods such as Mode may suffer from a threshold breakdown effect due to ‘subspace swaps’ [9, 10]. To reduce the risk that the signal subspace is fitted to noise in these cases, a modification was proposed in [11] consisting of using $p<m-r$ extra coefficients in $\mathbf{c}$. Then after computing the corresponding directions of arrival, all possible subsets of $r$ directions are compared using the maximum likelihood criterion and the best subset is chosen as the estimate. This method is called the ModeX in [11] and its principle is exactly what is used in [5] to propose the Enhanced-Puma.

Interestingly, while both papers [3] and [11] are referenced in [5], the equivalence (as shown above) of the Puma estimation criterion proposed there to Mode [3] and ModeX estimation criteria [11] was missed in [5].