<figure><img src="image_1.png" /></figure>

HAL
open science

# ML-DOA estimation using a sparse representation of array covariance

Thomas Aussaguès, Anne Ferréol, Alice Delmer, Pascal Larzabal

<figure><img src="image_2.png" /></figure>

## To cite this version:

Thomas Aussaguès, Anne Ferréol, Alice Delmer, Pascal Larzabal. ML-DOA estimation using a sparse representation of array covariance. 6th Junior Conference on Wireless and Optical Communications (JWOC), Université Paris-Saclay, Oct 2024, Gif-sur-Ivette, France. hal-04712924

HAL Id: hal-04712924

https://hal.science/hal-04712924v1

Submitted on 28 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ML-DOA estimation using a sparse representation of array covariance

Thomas Aussaguès *\*†*, Anne Ferréol*†*, Alice Delmer* and Pascal Larzabal †

*Thales, 4 avenue des Louvresses, 92230 Gennevilliers, France

†SATIE, Université Paris-Saclay, UMR CNRS 8029, 4 avenue des Sciences, 91190 Gif-sur-Yvette, France

**Abstract**—Sparse Direction-of-Arrival (DOA) estimators depend on the regularization parameter $\lambda$ which is often empirically tuned. In this work, conducted under the vectorized covariance matrix model, we are looking for theoretical equivalence between the Maximum Likelihood (ML) and sparse estimators. We show that under mild conditions, $\lambda$ can be chosen thanks to the distribution of the minimum of the ML criterion in the case of two impinging sources. The corresponding $\lambda$ choice is $\theta$-invariant, only requiring an upper bound on the number of sources. Furthermore, it guarantees the global minimum of the sparse $l_0$-regularized criterion to be the ML solution.

Numerical experiments confirm that, for the proposed $\lambda$, sparse and ML estimators yield the same statistical performance.

**Index Terms**—sparse estimation, regularization parameter, maximum likelihood

## I. INTRODUCTION

DOA estimation is a pivotal area in many critical applications such as radar or telecommunications. Throughout the last 60 years, a plethora of estimation techniques has been proposed [1], [2]. On the first hand conventional techniques such as Capon's beamformer [3] can not handle multiple sources. On the other hand, High Resolution (HR) estimators like MUSIC [4], [5] or ESPRIT [6] overcome this issue by exploiting the eigenstructure of the array covariance matrix. Nevertheless, these subspace-based approaches have limited performance in adversarial scenarios with few array snapshots or in presence of modeling errors [7]. Finally, the ML estimator [8] achieves the Cramér-Rao Lower Bound (CRLB) at high Signal-to-Noise Ratio (SNR) [9]. Although its statistical efficiency, the ML is rarely employed as it requires multidimensional non-convex optimization.

In recent years, there has been a growing interest within the signal processing processing community on sparse methods applied to DOA estimation [10] as they exhibit enhanced performances in tough scenarios. Out of the numerous modeling of the sparse DOA estimation problem, the sparse covariance matrix representation [11] emerged as a popular choice. Indeed, it benefits from the Virtual Array [12] concept thus having increased resolving power and number of identifiable sources.

For this model, the DOAs can be estimated through the minimization of a non-convex $l_0$-regularized objective $J_{l_0}$ which is parametrized by $\lambda$ the regularization parameter. Many works empirically tuned $\lambda$ which is unfeasible in practice. In [13], the authors proposed an interval for the regularization pa-

rameter. However, the interval bounds depends on the sources directions thus complexifying off-line selection.

In this work, we introduce a novel $\theta$-invariant regularization parameter choice $\lambda$ relying on the ML criterion distribution. This choice ensures the equivalence between ML and sparse DOA estimators as shown through numerical simulations.

## II. MAXIMUM LIKELIHOOD ON THE VECTORIZED COVARIANCE MATRIX MODEL

### A. Vectorized covariance matrix model

We consider a scenario of $M$ independent sources of directions $\Theta = \{\theta_1, \dots, \theta_M\}$ impinging on an array of $N$ sensors. Assuming narrowband hypothesis, the array observation is given by:

$$ \mathbf{x}(t) = \sum_{m=1}^{M} \mathbf{a}(\theta_m) s_m(t) + \mathbf{n}(t) = \mathbf{A}(\Theta) \mathbf{s}_{\Theta}(t) + \mathbf{n}(t) \quad (1) $$

with $\mathbf{A}(\Theta)$ the array manifold matrix formed by the steering vectors $\mathbf{a}(\theta_m) \in \mathbb{C}^N, m = 1 \dots M$, $\mathbf{s}_{\Theta}(t) = [s_1(t), \dots, s_M(t)]^T \in \mathbb{C}^M$ a vector containing the complex envelopes of the emitted signals and $\mathbf{n}(t) \in \mathbb{C}^N$ a complex Gaussian circular noise, independent of $\mathbf{s}(t)$, with covariance $\mathbb{E}[\mathbf{n}(t)\mathbf{n}^H(t)] = \sigma^2 \mathbf{I}_N$ where $\mathbb{E}[\cdot]$ is the temporal mean, $(\cdot)^H$ denotes the complex conjugate transpose and $\mathbf{I}_N$ the identity matrix of size $N$.

The array covariance matrix is then given by:

$$ \mathbf{R}_x = \mathbb{E}[\mathbf{x}(t)\mathbf{x}^H(t)] = \mathbf{A}(\Theta)\mathbf{R}_s\mathbf{A}^H(\Theta) + \sigma^2\mathbf{I}_N \quad (2) $$

where $\mathbf{R}_s = \mathbb{E}[\mathbf{s}_{\Theta}(t)\mathbf{s}_{\Theta}^H(t)]$ is the sources covariance matrix.

In practice, the covariance matrix is not available. Instead, $\mathbf{R}_x$ is replaced by its estimate $\hat{\mathbf{R}}_x$ computed using $K$ identically and independently distributed (i.i.d) array snapshots. Assuming temporally white noise ($\forall i \neq j, \mathbb{E}[\mathbf{n}^H(t_i)\mathbf{n}(t_j)] = 0$), $\hat{\mathbf{R}}_x$ can be written as:

$$ \hat{\mathbf{R}}_x = \frac{1}{K} \sum_{k=1}^{K} \mathbf{x}(t_k)\mathbf{x}^H(t_k) = \mathbf{R}_x + \Delta\mathbf{R}_x \quad (3) $$

where $K\Delta\mathbf{R}_x$ is a Wishart noise due to the finite number of snapshots [14].
Assuming that sources are uncorrelated, the vectorized covariance matrix model is obtained from (2, 3) by applying the column-wise vectorization operator $\text{vec}(\cdot)$ to $\hat{\mathbf{R}}_x - \sigma^2 \mathbf{I}_N$:

$$ \mathbf{r} = \text{vec}(\hat{\mathbf{R}}_x - \sigma^2 \mathbf{I}_N) = \sum_{m=1}^{M} \mathbf{b}(\theta_m)\gamma_m + \delta = \mathbf{B}(\Theta)\gamma_{\Theta} + \delta \quad (4) $$

where $\mathbf{b}(\theta_m) = \mathbf{a}^*(\theta_m) \otimes \mathbf{a}(\theta_m)$ is the virtual array steering vector with $\otimes$ the Kronecker product and $(\cdot)^*$ the complex conjugate. The corresponding mixing matrix is referred as $\mathbf{B}(\Theta) = [\mathbf{b}(\theta_1), \dots, \mathbf{b}(\theta_M)]$. $\gamma_m = \mathbb{E}[|s_m(t)|^2]$ denotes the power of the m-th source and $\gamma_{\Theta} = [\gamma_1, \dots, \gamma_M]^T$ the vector of the sources powers. Finally, $\delta = \text{vec}(\Delta\mathbf{R}_x)$ is the vectorized noise.

Asymptotically, the vectorized observation asymptotically follows a multivariate complex Gaussian law $\mathbb{CN}(\mu, \Gamma, C)$ where [15]:

$$ \Gamma = \mathbb{E} \left[ (\mathbf{r} - \mu) (\mathbf{r} - \mu)^H \right] = \frac{1}{K} \mathbf{R}_x^T \otimes \mathbf{R}_x \quad (5) $$

and C is defined in [15].

## B. Maximum Likelihood estimation

For additive Gaussian noise a ML estimator of $\Theta$, known to achieve the Cramér-Rao lower bound at high Signal-to-Noise Ratio (SNR), can be easily obtained from (4) [8]. However, for non-diagonal covariance matrices such as (5) this yields a multi-term log-likelihood function. This objective can be simplified for diagonal covariance matrices leading us to apply a pre-whitening step to (4):

$$ \mathbf{y} = \hat{\mathbf{W}}\mathbf{r} = \hat{\mathbf{W}}\mathbf{B}(\Theta)\gamma_{\Theta} + \delta_w \quad (6) $$

where $\hat{\mathbf{W}} = \sqrt{K}\hat{\mathbf{R}}_x^{-T/2} \otimes \hat{\mathbf{R}}_x^{-1/2}$ is the estimated whitening matrix and $\delta_w = \hat{\mathbf{W}}\delta$ a spatially white noise with $\mathbb{E}[\delta_w \delta_w^H] = \mathbf{I}_{N^2}$. Applying the ML to (6) yields the following estimator of $\Theta$:

$$ \begin{gathered} J_{ML}(\Theta) = \text{tr} (\Pi^{\perp}(\Theta)yy^H) = y^H \Pi^{\perp}(\Theta)y \quad (7) \\ \hat{\Theta} = \arg \min_{\Theta \in \mathbb{R}^M} J_{ML}(\Theta) \end{gathered} $$

where $\Pi^{\perp}(\Theta) = \mathbf{I}_{N^2} - (\hat{\mathbf{W}}\mathbf{B}(\Theta))(\hat{\mathbf{W}}\mathbf{B}(\Theta))^{\#}$ is the noise projector computed for directions $\Theta$ with $(\cdot)^{\#}$ the Moore-Penrose pseudo inverse and $\text{tr}(\cdot)$ the trace operator. In what follows, the optimal value of the ML criterion $J_{ML}$ is referred as:

$$ \epsilon = \min_{\Theta \in \mathbb{R}^M} J_{ML}(\Theta) = J_{ML}(\hat{\Theta}) = y^H \Pi^{\perp}(\hat{\Theta})y \quad (8) $$

## III. SPARSE ESTIMATION

### A. Sparse modeling

Under the hypothesis that DOAs lie within a predefined grid of directions $\mathbf{\Phi} = \{\varphi_1, \dots, \varphi_G\}$ ($\Theta \subset \mathbf{\Phi}$), sparsity can be introduced in (6) using an overcomplete dictionary $\mathbf{B}(\mathbf{\Phi}) = [\mathbf{b}(\varphi_1), \dots, \mathbf{b}(\varphi_G)] \in \mathbb{C}^{N^2 \times G}$ with $G \gg N^2$. Thus, we obtain the following sparse model for the observation $\mathbf{y}$:

$$ \mathbf{y} = \hat{\mathbf{W}}\mathbf{B}(\mathbf{\Phi})\gamma_0 + \delta_w \quad (9) $$

where $\gamma_0$ is a sparse vector which has exactly $M$ non-zeros entries corresponding to the sources powers. Thereafter, $\gamma_0$ is referred as the sparse spectrum.

### B. Sparse problem formulation

The aim of DOA estimation is to find the directions of the impinging signals from the observation $\mathbf{y}$ which corresponds, in the sparse model (9), to the non-zeros entries of $\gamma_0$. Thus, we need to estimate the sparse spectrum $\gamma_0$ in order to provide an estimate of $\Theta$. This problem is ill-posed since $G \gg N^2$. Hence, the sparsity of $\gamma_0$ is exploited to ensure the uniqueness of the solution. Following Delmer's work [13], an estimate of $\gamma$ can be obtained through the minimization of the $\ell_0$-regularized objective function given by:

$$ \min_{\gamma \in \mathbb{C}^G} \left\{ J_{\ell_0}(\lambda, \gamma) = \frac{1}{2} \|\mathbf{y} - \hat{\mathbf{W}}\mathbf{B}(\mathbf{\Phi})\gamma\|_2^2 + \lambda \|\gamma\|_0 \right\} \quad (10) $$

with a suitable regularization parameter $\lambda > 0$ which balances data fidelity towards the solution sparsity.

## IV. SELECTION OF THE REGULARIZATION PARAMETER

The choice of $\lambda$ in (10) is of utmost importance since it controls the trade-off between estimation error and sparsity. Many approaches exploiting the noise level have been proposed for $\ell_1$-regularization [16]–[18]. For the $\ell_0$ optimization framework $\lambda$ is generally empirically tuned given that only few results are available. Recently, Delmer [13] et al. proposed to choose $\lambda$ so that $J_{\ell_0}$ has a global minimum in $\hat{\gamma}$ ie.:

$$ \lambda > 0 \mid \forall \gamma \in \mathbb{C}^G, J_{\ell_0}(\lambda, \hat{\gamma}) \le J_{\ell_0}(\lambda, \gamma) \quad (11) $$

where $\hat{\gamma}$ is the sparse spectrum, in the neighboring of $\gamma_0$, obtained using the ML DOA estimate $\hat{\Theta}$ (7). For a given observation $\mathbf{y}$, the values of $\lambda$ that achieves condition (11) are within a stochastic admissible interval $I_M(\mathbf{y}) = [\lambda_M^-, \lambda_M^+]$. This interval ensures the global minimizer of (10) to be the ML solution which is an $M$-sparse vector.

For the case $M=2$, the corresponding interval is $I_2(\mathbf{y}) = [\lambda_2^-, \lambda_2^+]$ where:

$$ \lambda_2^- \approx \frac{c_2(\mathbf{y}) - c_3(\mathbf{y})}{2} \quad \text{and} \quad \lambda_2^+ \approx \frac{c_1(\mathbf{y}) - c_2(\mathbf{y})}{2} \quad (12) $$

with $c_k(\mathbf{y}) = \inf\{\|\mathbf{y} - \hat{\mathbf{W}}\mathbf{B}(\mathbf{\Phi})\gamma\|_2^2, \|\gamma\|_0 = k\}$ the ML criterion optimal value (8) for a prescribed sparsity level $k$. Note that we have $c_3(\mathbf{y}) \le c_2(\mathbf{y}) \le c_1(\mathbf{y})$. In the following, $I_2$ is assumed non-empty (ie. $c_1(\mathbf{y}) > c_3(\mathbf{y})$).

For arrays presenting robustness against second order ambiguities for $M=2$ sources with equal power, the relation $c_2(\mathbf{y}) \ll c_1(\mathbf{y})$ holds. Thus, the following inequality can be proven:

$$ \lambda_2^- \le \frac{1}{2} c_2(\mathbf{y}) \le \lambda_2^+ \quad (13) $$

Remarking that $c_2(\mathbf{y})$ corresponds to the minimum of the ML criterion (8) $\epsilon$, a suitable regularization parameter choice satisfying (13) can be formulated as:

$$ \lambda = \frac{1}{2} F_{\epsilon}(\eta) \quad (14) $$
## VI. CONCLUSION

In this paper, we proposed a novel regularization choice for $ℓ_0$-regularized DOA estimation by exploiting the distribution of the ML criterion thus leading to $\theta$-invariant choice. Through numerical experiments, we confirmed that the proposed $\lambda$ belongs to Delmer's interval hence ensuring the equivalence with the ML estimator. This approach must be extended to scenarios with $M > 2$ sources. This is an ongoing work.

## REFERENCES

[1] H. Krim and M. Viberg, "Two decades of array signal processing research: The parametric approach," *Signal Processing Magazine, IEEE*, vol. 13, pp. 67–94, 08 1996.

[2] M. Pesavento, M. Trinh-Hoang, and M. Viberg, "Three more decades in array signal processing research: An optimization and structure exploitation perspective," *IEEE Signal Processing Magazine*, vol. 40, no. 4, pp. 92–106, 2023.

[3] J. Capon, "High-resolution frequency-wavenumber spectrum analysis," *Proceedings of the IEEE*, vol. 57, no. 8, pp. 1408–1418, 1969.

[4] G. Bienvenu and L. Kopp, "Optimality of high resolution array processing using the eigensystem approach," *IEEE Transactions on Acoustics, Speech, and Signal Processing*, vol. 31, no. 5, pp. 1235–1248, 1983.

[5] R. Schmidt, "Multiple emitter location and signal parameter estimation," *IEEE Transactions on Antennas and Propagation*, vol. 34, no. 3, pp. 276–280, 1986.

[6] R. Roy and T. Kailath, "Esprit-estimation of signal parameters via rotational invariance techniques," *IEEE Transactions on Acoustics, Speech, and Signal Processing*, vol. 37, no. 7, pp. 984–995, 1989.

[7] B. Friedlander, "A sensitivity analysis of the music algorithm," in *International Conference on Acoustics, Speech, and Signal Processing*, 1989, pp. 2811–2814 vol.4.

[8] B. Ottersten, M. Viberg, P. Stoica, and A. Nehorai, *Exact and Large Sample Maximum Likelihood Techniques for Parameter Estimation and Detection in Array Processing*. Berlin, Heidelberg: Springer Berlin Heidelberg, 1993, pp. 99–151.

[9] P. Stoica and A. Nehorai, "Music, maximum likelihood, and cramer-rao bound," *IEEE Transactions on Acoustics, Speech, and Signal Processing*, vol. 37, no. 5, pp. 720–741, 1989.

[10] Z. Yang, J. Li, P. Stoica, and L. Xie, "Sparse methods for direction-of-arrival estimation," 2017.

[11] J. S. Picard and A. J. Weiss, "Direction finding of multiple emitters by spatial sparsity and linear programming," in *2009 9th International Symposium on Communications and Information Technology*, 2009, pp. 1258–1262.

[12] P. Chevalier, L. Albera, A. Ferreol, and P. Comon, "On the virtual array concept for higher order array processing," *IEEE Transactions on Signal Processing*, vol. 53, no. 4, pp. 1254–1271, 2005.

[13] A. Delmer, A. Ferréol, and P. Larzabal, "On regularization parameter for 10-sparse covariance fitting based doa estimation," in *ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, 2020, pp. 4552–4556.

[14] N. R. Goodman, "Statistical analysis based on a certain multivariate complex gaussian distribution (an introduction)," *The Annals of mathematical statistics*, vol. 34, no. 1, pp. 152–177, 1963.

[15] M. Mahot, F. Pascal, P. Forster, and J.-P. Ovarlez, "Asymptotic properties of robust complex covariance matrix estimates," *IEEE Transactions on Signal Processing*, vol. 61, no. 13, pp. 3348–3356, 2013.

[16] D. Malioutov, M. Cetin, and A. Willsky, "A sparse signal reconstruction perspective for source localization with sensor arrays," *IEEE Transactions on Signal Processing*, vol. 53, no. 8, pp. 3010–3022, 2005.

[17] Z. He, Q. Liu, L. Jin, and S. Ouyang, "Low complexity method for doa estimation using array covariance matrix sparse representation," *Electronics Letters*, vol. 49, no. 3, pp. 228–230, 2013.

[18] W. Cui, T. Qian, and J. Tian, "Enhanced covariances matrix sparse representation method for doa estimation," *Electronics Letters*, vol. 51, no. 16, pp. 1288–1290.

[19] T. Aussaguès, A. Ferréol, A. Delmer, and P. Larzabal, "Looking for equivalence between maximum likelihood and sparse doa estimators," in *2024 32th European Signal Processing Conference (EUSIPCO)*, 2024.

[20] P. L. Combettes and J.-C. Pesquet, "Proximal splitting methods in signal processing," 2010.

where $F_e(\eta)$ is the inverse Cumulative Distribution Function of $\epsilon$ for probability $\eta$. In [19], we derived the distribution of $\epsilon$ under complex Gaussian non-circular noise. We proved that $\epsilon$ asymptotically has a $\chi_2$ distribution with $N^2 - M$ degrees of freedom. The choice of $\eta$ is discussed in section V.

## V. NUMERICAL EXPERIMENTS

We consider an array of $N = 6$ antennas with 5 antennas distributed around a circle of radius $0.8\lambda_0$ where $\lambda_0$ denotes the wavelength and one central sensor. $M = 2$ sources of directions $\theta_1 = 180^\circ, \theta_2 = 200^\circ$ with $K = 200$ and SNR = 0 dB impinge on this array. Fig. 1 depicts the distribution of the interval bounds $\lambda_2^-, \lambda_2^+$ and $\lambda = \frac{1}{2}\epsilon$. For this scenario, our regularization parameter choice belongs to the interval thus ensuring that $J_{\ell_0}$ global minimizer corresponds to the ML DOA solution.

<figure><img src="image_4.png" /><figcaption>Fig. 1. Histograms of the distributions of $\lambda_2^-, \lambda_2^+$ and $\lambda = \frac{1}{2}\epsilon$ with Cumulative Distribution Function (CDF) normalization.</figcaption></figure>

According to [19], we propose to choose $\eta = 0.05$ leading to $\lambda = \frac{1}{2}F^{-1}_{\chi^2(34)}(0.05)$ where $F^{-1}_{\chi^2(34)}$ is the inverse CDF of the $\chi_2$ law with 34 degrees of freedom.

Finally, we verify that, for the proposed $\lambda$ choice, ML and sparse methods lead to the same statistical performance. To this end, we conducted 10000 independent Monte-Carlo experiments for SNR values between -20 dB and 20 dB. The sources are considered resolved if two peaks are detected and $\max\{|\hat{\theta}_1 - \theta_1|, |\hat{\theta}_2 - \theta_2|\} < 10^\circ$ where $10^\circ$ corresponds to the half beamwidth of the virtual array. Results, obtained with the Forward-Backward splitting algorithm [20], are represented on Fig.2. For SNR $\ge -4$ dB, ML and sparse estimators have the same statistical performance thus validating the proposed $\lambda$ choice.

<figure><img src="image_5.png" /><figcaption>Fig. 2. Probability of detection (top) and RMSE (bottom). Given that the two sources have the same power, results are only represented for direction $\theta_1$.</figcaption></figure>