# L0 Regularization parameter for sparse DOA estimation of coherent signals with modeling errors

Alice Delmer, Anne Ferréol, Pascal Larzabal

► To cite this version:

Alice Delmer, Anne Ferréol, Pascal Larzabal. L0 Regularization parameter for sparse DOA estimation of coherent signals with modeling errors. Signal Processing, 2023, 209, pp.1-7. 10.1016/j.sigpro.2023.109006. hal-04164606

HAL Id: hal-04164606

https://cnam.hal.science/hal-04164606v1

Submitted on 31 Mar 2025

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_1.png" /></figure>

Distributed under a Creative Commons CC0 - Public Domain Dedication 4.0 International License
# L0 Regularization Parameter for Sparse DOA Estimation of Coherent Signals with Modeling Errors

Alice Delmer<sup>a,b,*</sup>, Anne Ferréol<sup>a,b</sup>, Pascal Larzabal<sup>b</sup>

## Abstract

Sparse methods have been recently introduced in Direction-Of-Arrival (DOA) estimation as an alternative to subspace-based methods and Maximum Likelihood (ML) techniques. This paper proposes a low complexity sparse L0-regularized method for DOA estimation of coherent signals in multipath environment, a scenario in which subspace-based methods fail. The first eigenvector of the covariance matrix is used as an observation and its statistics in the presence of modeling errors are derived. Thanks to those statistics, we present a theoretical statistical analysis of an interval in which the regularization parameter, usually empirically tuned, should stand. For a regularization parameter in this interval, the global solution of the L0-regularized problem coincides with the solution of the deterministic ML. An off-line selection of the regularization parameter is therefore proposed. Simulations confirm the relevance of this interval and of the proposed method.

**Keywords:** Direction-of-arrival estimation, regularization parameter, sparse methods, performance analysis

## 1. Introduction

In the field of Direction-Of-Arrival (DOA) estimation, high-resolution methods such as MUSIC have brought many improvements compared to conventional beamforming. However, it still suffers from known limitations in the case of close sources and in environment with multipath propagation [1]. It even fails in scenario with coherent signals because the signal subspace is rank one and not the number of signals deficient. An alternative is to used a preprocessing technique as spatial smoothing [2] to decorrelate coherent signals. However, such process reduces the resolution and the degree of freedom of an array system. Maximum Likelihood (ML) methods can be used in this case in the presence of

<footnote><sup>a</sup>Thales, 4 avenue des Louvresses, 92230 Gennevilliers, France</footnote>

<footnote><sup>b</sup>SATIE, Université Paris-Saclay, CNRS, 4 avenue des Sciences, 91190 Gif-sur-Yvette, France</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Email address: alice.delmer@ens-paris-saclay.fr (Alice Delmer)</footnote>
multipath propagation at the cost of numerical high computational complexity. But recently, sparse DOA estimation has shown great potential to overcome some drawbacks of subspace-based methods [3, 4, 5] and with a low complexity. The purpose of this paper is thus to establish a sparse framework for the presence of coherent signals that allows to overcome ML methods drawbacks.

Sparsity naturally arises in DOA estimation when considering a calibration table $\mathbf{A} = [\mathbf{a}(\psi_1), \dots, \mathbf{a}(\psi_G)]$ containing a collection of array responses for $G$ angles belonging on a predefined grid. This calibration table acts as a dictionary for on-grid sparse DOA estimation. Considering a noisy observation vector $\tilde{\mathbf{u}} = \mathbf{A}\mathbf{z}_0 + \delta\mathbf{u}$, where $\mathbf{z}_0$ is the true discretized angular spectrum, incident angles can be retrieved from the resolution of the $\ell_0$-constrained problem

$$ \min_{\mathbf{z}} \|\mathbf{A}\mathbf{z} - \tilde{\mathbf{u}}\|_2^2, \quad \text{s.t.} \quad \|\mathbf{z}\|_0 \le M, \qquad (1) $$

where the $\ell_0$-pseudo norm is defined as $\|\mathbf{z}\|_0 = \text{Card}\{g \in \{1, \dots, G\} : z_g \neq 0\}$. The solution of this constrained problem coincides with the deterministic Maximum Likelihood (ML) in the presence of $M$ signals. However, this sparse $\ell_0$ minimization problem is known to be NP-hard: its optimal resolution requires an exhaustive search, which is prohibitive in terms of computational complexity in the presence of several sources or multipath propagation. The largely used $\ell_1$ approximation is not applicable for an operational scenario with a low number $N$ of antennas ($N \le 7$) as the mutual coherence of the dictionary $\mathbf{A}$ is thus close to 1. Suboptimal minimization methods have thus been proposed as Iterative Hard Thresholding (IHT), aiming to estimate the solution of the following $\ell_0$-regularized problem: Therefore, more efficient approaches are needed. Convex relaxation is largely used: replacing the $\ell_0$ norm by the $\ell_1$ norm in (1) leads to the well known Basis Pursuit Denoising (BPDN) and the $\ell_1$ optimization problem can be solved in a polynomial time [6]. But successful recovery of $\mathbf{z}_0$ can not be guaranteed in an operational multi-sources or multipath scenario with a low number $N$ of antennas ($N \le 7$) as the mutual coherence of the dictionary $\mathbf{A}$ is thus close to 1 [7]. Interesting suboptimal low complexity algorithms have also been proposed. For example, iterative hard thresholding (IHT) aims to solve the following $\ell_0$-regularized problem:

$$ \min_{\mathbf{z}} \left\{ J_{\ell_0}(\lambda, \mathbf{z}) = \frac{1}{2} \|\mathbf{A}\mathbf{z} - \tilde{\mathbf{u}}\|_2^2 + \lambda \|\mathbf{z}\|_0 \right\}. \qquad (2) $$

There still are several difficulties in solving this problem: the $\ell_0$-regularized error surface presents numerous local minima, and convergence of iterative descent algorithms is usually only proved to stationary points. Moreover, However, the solution highly depends on the regularization parameter $\lambda$, which is usually empirically tuned. Several heuristic methods for tuning this parameter were introduced in [8]: in order to gain in computation time, performing the minimization for several values of $\lambda$ is not desired. Homotopy-based methods are empirically evaluated in [9], but those methods are not applicable in the case of closely spaced signals with a highly coherent dictionary. There is a need for a theoretical analysis of the regularization parameter $\lambda$ to ensure that the solution of this problem is the same
as that of the initial $\ell_0$-constrained problem and, consequently, as that of the
deterministic ML.

Nikolova theoretically provided an admissible interval for the regularization parameter so that the minimizer of has a given sparsity. However, Nikolova's interval is conditional to a realization of $\tilde{u}$ and its computation is combinatorial, so it is useless in practical applications. In order to extend this idea for a practical implementation, there is a need to study the stochastic behavior of the bounds of this interval. Such investigation has been started in, where Nikolova's work has been generalized to a stochastic observation which was the vectorized covariance matrix in the presence of independent sources.

In this paper, the focus is the DOA estimation of coherent signals in multipath environment. This is why we The purpose of this paper is to establish a robust sparse framework for the presence of coherent signals that comes closer to ML methods performance while overcoming its drawbacks. Moreover, we consider the case of modeling errors: the calibration table A considered in the algorithms does not correspond exactly to the measurement system. Large degradation of performance are indeed induced by modeling errors for conventional methods [10, 11], and they particularly affect the resolution [12]: contributions of sparse methods in this case can thus be of great interest. Contributions are threefold. a) We originally use the eigenvector of the covariance matrix corresponding to the highest eigenvalue as an observation vector $\tilde{u}$ for sparse DOA estimation. This reduced observation allows to use the Single Measurement Vector (SMV) framework while exploiting several snapshots. Sparse DOA estimation is indeed less expensive with a SMV observation than with a Multiple Measurement Vector (MMV) observation. Moreover, we consider the case of modeling errors: the calibration table A considered in the algorithms does not correspond exactly to the measurement system. Large degradation of performance are indeed induced by modeling errors for conventional methods, and they particularly affect the resolution. b) We provide a statistical investigation for the determination of an admissible interval for the regularization parameter $\lambda$ of the $\ell_0$-regularized problem so that the global solution of this problem coincides with the ML solution. Nikolova [13] theoretically provided an admissible interval for the regularization parameter so that the minimizer of (2) has a given sparsity. However, Nikolova's interval is conditional to a realization of $\tilde{u}$ and its computation is combinatorial, so it is useless in practical applications. In order to extend this idea for a practical implementation, there is a need to study the stochastic behavior of the bounds of this interval. Such investigation has been started in [14, 15], where Nikolova's work has been generalized to a stochastic observation which was the vectorized covariance matrix in the presence of independent sources. A specific statistical investigation is provided in this paper for $\ell_0$-sparse bearing estimation with the eigenvector as an observation, and in the presence of modeling errors. Sparse bearing estimation with this new observation requires a specific statistical investigation for the determination of an admissible interval for the regularization parameter: this is the purpose of this paper. c) We propose a strategy for an off-line selection of the regularization parameter $\lambda$ in the admissible interval. This a priori selection additionnaly reduces the complexity compared to methods performing the resolution for several values of $\lambda$ and selecting a posteriori the best estimate. We show that the chosen value, used in
an existing low complexity algorithm, allows to come closer to the performance of ML.

We first provide the statistics of the eigenvector in the presence of modeling errors, and derive a theoretical interval for the regularization parameter $\lambda$ valid in spite of the statistical fluctuations of the observation. We study the influence of $\lambda$ on the behavior of a $l_0$-minimization algorithm and relate the obtained performance to the previous interval. Based on this study, we propose a strategy to perform an off-line selection of the parameter $\lambda$ to be used in the algorithms. This approach is validated by simulation. The paper is organized as follows: In Section 2, we present the sparse modeling of the eigenvector for coherent signals, and derive its statistics in the presence of modeling errors. In Section 3, we provide a theoretical interval for the regularization parameter $\lambda$ valid in spite of the statistical fluctuations of the observation. This theoretical interval is not dependent of the minimization algorithm used afterwards. In Section 4, we study the influence of $\lambda$ on the behavior of a $l_0$-minimization algorithm and relate the obtained performance to the previous interval. We perform the off-line selection of the parameter $\lambda$ and validate our approach by simulation.

**Notations.** Upper-case and lower-case boldface latin letters denote matrices and vectors, respectively. $^{\cdot}$ denotes the conjugate. $^{T}$ the transpose and $^{H}$ the conjugate transpose of a vector or matrix. $x_i$ is the $i^{th}$ component of vector $\mathbf{x}$.

## 2. SMV observation for coherent signals DOA estimation

### 2.1. Sparse representation

Let us consider an arbitrary array of $N$ antennas, and $M$ narrow-band signals impinging on this array from angles $\theta_m$, $m = 1...M$. Let us define $\Theta = \{\theta_1, ..., \theta_M\}$ as the set of DOAs. The output signal $\mathbf{x}(t) \in \mathbb{C}^N$ is a noisy linear combination of the array responses $\tilde{\mathbf{a}}(\theta_m)$ (called steering vectors):

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

where $\tilde{s}_m(t)$ is the complex envelope of the $m$-th incoming signal, $\mathbf{n}(t)$ a white Gaussian noise of zero mean and covariance $\mathbb{E}[\mathbf{n}(t)\mathbf{n}^H(t)] = \sigma_n^2 \mathbf{I}_N$, $\mathbf{I}_N$ is the $N \times N$ identity matrix, and $\tilde{\mathbf{A}}_{\Theta}$ is the mixing matrix $\tilde{\mathbf{A}}(\theta_1), ..., \tilde{\mathbf{A}}(\theta_M)$.

Considering multipath propagation from a single transmitter, the complex envelopes of the received signals verify $\tilde{s}_m(t) = q_m s(t - \tau_m)$ with $\mathbb{E}[|s(t)|^2] = \gamma$. In an urban environment, signals might be quasi-coherent, so that $\tau_1 \approx \tau_2 \approx \dots \approx \tau_M$, and the covariance matrix $\mathbf{R}_{ss} = \mathbb{E}[\tilde{s}(t)\tilde{s}^H(t)]$ can be computed as $\mathbf{R}_{ss} \approx \gamma \mathbf{q} \mathbf{q}^H$, for $\mathbf{q} = [q_1, ..., q_M]^T$.

The covariance matrix $\tilde{\mathbf{R}}_{xx} = \mathbb{E}[\mathbf{x}(t)\mathbf{x}^H(t)]$ of the output signal verifies

$$ \tilde{\mathbf{R}}_{xx} = \tilde{\mathbf{A}}_{\Theta} \mathbf{R}_{ss} \tilde{\mathbf{A}}_{\Theta}^H + \sigma_n^2 \mathbf{I}_N. \quad (4) $$
By noting $\tilde{\mathbf{R}} = \tilde{\mathbf{R}}_{xx} - \sigma_n^2 \mathbf{I}_N = \tilde{\mathbf{A}}_{\Theta} \mathbf{R}_{ss} \tilde{\mathbf{A}}_{\Theta}^H$, we obtain $\tilde{\mathbf{R}} = \mu_1 \tilde{\mathbf{u}} \tilde{\mathbf{u}}^H$ with $\mu_1 = \gamma \mathbf{q}^H \tilde{\mathbf{A}}_{\Theta}^H \tilde{\mathbf{A}}_{\Theta} \mathbf{q}$ and

$$ \tilde{\mathbf{u}} = \frac{1}{\sqrt{\mathbf{q}^H \tilde{\mathbf{A}}_{\Theta}^H \tilde{\mathbf{A}}_{\Theta} \mathbf{q}}} \tilde{\mathbf{A}}_{\Theta} \mathbf{q}. \qquad (5) $$

$\tilde{\mathbf{u}}$ is the eigenvector associated with the highest (and only nonzero) eigenvalue $\mu_1$ of $\tilde{\mathbf{R}}$.

In the framework of on-grid sparse estimation, it is assumed that the directions of arrival $\theta_m$ of the signals belong to a prescribed grid. Let us note $\Psi$ the set of the G candidate angles of arrival resulting from the discretization of the field-of-view: $\Psi = \{\psi_1, \dots, \psi_G\}$. It is thus assumed that $\Theta \subset \Psi$. Steering vectors of the candidate angles are concatenated in the matrix $\tilde{\mathbf{A}} = [\tilde{\mathbf{a}}(\psi_1), \dots, \tilde{\mathbf{a}}(\psi_G)]$. The mixing matrix $\tilde{\mathbf{A}}_{\Theta}$ is the submatrix of $\tilde{\mathbf{A}}$ containing only the columns corresponding to incident angles. The first eigenvector $\tilde{\mathbf{u}}$ is thus a sparse linear combination of the steering vectors $\tilde{\mathbf{a}}(\psi_g)$, for $\psi_g \in \Psi$:

$$ \tilde{\mathbf{u}} = \tilde{\mathbf{A}} \tilde{\mathbf{z}}, \qquad (6) $$

where $\tilde{\mathbf{z}}$ is the sparse vector whose non zero components are the components of $\tilde{\mathbf{z}}_{\Theta}$ defined as $\tilde{\mathbf{z}}_{\Theta} = (1/\sqrt{\mathbf{q}^H \tilde{\mathbf{A}}_{\Theta}^H \tilde{\mathbf{A}}_{\Theta} \mathbf{q}}) \mathbf{q}$. Consequently, we propose to use $\tilde{\mathbf{u}}$ as an observation for sparse DOA estimation. In an operational system, because of noise and modeling errors, $\tilde{\mathbf{u}}$ is subject to statistical fluctuations.

## 2.2. Statistics of the eigenvector in the presence of modeling errors

In an operational system, the steering vectors $\tilde{\mathbf{a}}(\psi_g)$ of the data collection system might not match exactly the alleged steering vectors $\mathbf{a}(\psi_g)$ used in the algorithms. A number of different modeling errors indeed coexist, such as errors in the sensors gains, sensors phases and sensors locations, but also mutual coupling. We model those errors as additive perturbations on the components of the nominal steering vectors:

$$ \tilde{\mathbf{A}} = \mathbf{A} + \mathbf{E} \qquad (7) $$

with $\mathbf{A} = [\mathbf{a}(\psi_1), \dots, \mathbf{a}(\psi_G)]$ the considered calibration matrix, with $\|\mathbf{a}(\psi_g)\|^2 = 1$, and $\mathbf{E} = [\mathbf{e}(\psi_1), \dots, \mathbf{e}(\psi_G)]$ the modeling error matrix. Vectors $\mathbf{e}(\psi_g)$ are assumed to be Gaussian, circular and independent on each sensor and each sources such that: $\mathbb{E}[\mathbf{e}(\psi_g)\mathbf{e}(\psi_j)^H] = \delta_{gj}\sigma_e^2/N\mathbf{I}_N$ and $\mathbb{E}[\mathbf{e}(\psi_g)\mathbf{e}(\psi_j)^T] = 0$. This perturbation model is often considered [16, 17, 18], and it is general enough to get a good picture of the performance of an approach.

In this article, we assume that the covariance matrix $\tilde{\mathbf{R}}_{xx}$ is perfectly known. This allows to concentrate on the degradation due to modeling uncertainties and avoids mixing them with discrepancies induced by the finite number of samples.

Let us note $\mathbb{E}_{\Theta}$ the submatrix of $\mathbb{E}$ containing only the columns corresponding to incident angles. Thus from (7), we have $\tilde{\mathbf{A}}_{\Theta} = \mathbf{A}_{\Theta} + \mathbb{E}_{\Theta}$. By straightforward calculations whose only difficulty is the length and boredom, the first
165 order Taylor expansion of $\tilde{u}$ as defined in (5) is expressed as:

$$ \tilde{u} \approx u + \delta u, \qquad (8) $$

with

$$ u = \frac{1}{\sqrt{q^H A_{\Theta}^H A_{\Theta} q}} A_{\Theta} q, \qquad (9) $$

$$ \delta u = \frac{1}{2 (q^H A_{\Theta}^H A_{\Theta} q)^{3/2}} \left( 2 q^H A_{\Theta}^H A_{\Theta} q E_{\Theta} q - q^H (E_{\Theta}^H A_{\Theta} + A_{\Theta}^H E_{\Theta}) q A_{\Theta} q \right). \qquad (10) $$

and $z_0$ is the sparse vector whose non-zero components are the components of $z_{\Theta} = (1/\sqrt{q^H A_{\Theta}^H A_{\Theta} q}) q$.

As $E_{\Theta}$ is Gaussian, thus $\tilde{u}$ first order approximation is complex Gaussian with mean $u$, and second order statistics $\Gamma = E[\delta u \delta u^H]$ and $C = E[\delta u \delta u^T]$ given by

170

$$ \Gamma = \frac{\sigma_e^2 q^H q}{2 q^H A_{\Theta}^H A_{\Theta} q} \left[ 2 I_N - \frac{A_{\Theta} q q^H A_{\Theta}^H}{q^H A_{\Theta}^H A_{\Theta} q} \right], \qquad (11) $$

$$ C = \frac{-\sigma_e^2 q^H q}{2 q^H A_{\Theta}^H A_{\Theta} q} \times \frac{A_{\Theta} q q^T A_{\Theta}^T}{q^H A_{\Theta}^H A_{\Theta} q}. \qquad (12) $$

Let us note that the covariance matrix $\Gamma$ is phase independent while the relation matrix $C$ is not.

Let us notice that $u = Az_0$ where $z_0$ is the sparse vector whose non-zero components are the components of $z_{\Theta} = (1/\sqrt{q^H A_{\Theta}^H A_{\Theta} q}) q$. By solving the sparse problem (2), it is thus possible to retrieve the DOAs of the impinging signals.

175

### 3. Admissible interval for the regularization parameter

The global minimum of the $l_0$-regularized problem (2) depends on the regularization parameter $\lambda$. The goal of this work is to define an admissible interval for $\lambda$ so that the global optimum corresponds to the ML solution and valid in spite of the statistical fluctuations of the proposed observation $\tilde{u}$ given by (8).

180

#### 3.1. Background

Nikolova [13] provides, conditionally to a realization, a theoretical expression of an interval in which the regularization parameter $\lambda$ should belong so that the $l_0$-regularized problem (2) admits the same global minimizer than the $l_0$-constrained problem (1). In this work, we denote $I(\tilde{u}) = (\lambda_M^-(\tilde{u}), \lambda_M^+(\tilde{u}))$ the interval computed for a given realization of $\tilde{u}$. The optimal sets of problems (1)
and (2) coincide for a given realization $\tilde{u}$ if and only if $\lambda_M^-(\tilde{u}) < \lambda < \lambda_M^+(\tilde{u})$ [13] with:

$$ \lambda_M^-(\tilde{u}) = \max_k \left\{ \frac{c_M(\tilde{u}) - c_k(\tilde{u})}{2(k-M)} \middle| M < k \le N \right\} \quad (13) $$

$$ \lambda_M^+(\tilde{u}) = \min_k \left\{ \frac{c_k(\tilde{u}) - c_M(\tilde{u})}{2(M-k)} \middle| 0 \le k < M \right\} \quad (14) $$

where $c_k$ is the optimal value of the least-square residual for a k-sparse vector:

$$ c_k(\tilde{u}) = \min_z \{ \| \mathbf{A} \mathbf{z} - \tilde{u} \|_2^2 \mid \mathbf{z} \in \mathbb{C}^G \text{ and } \| \mathbf{z} \|_0 \le k \} \quad (15) $$

$$ = \min_{\mathbf{\Phi}_k} \tilde{u}^H \mathbf{\Pi}_n(\mathbf{\Phi}_k) \tilde{u}. \quad (16) $$

$\mathbf{\Pi}_n(\mathbf{\Phi}_k) = \mathbf{I}_N - \mathbf{A}_{\mathbf{\Phi}_k} [\mathbf{A}_{\mathbf{\Phi}_k}^H \mathbf{A}_{\mathbf{\Phi}_k}]^{-1} \mathbf{A}_{\mathbf{\Phi}_k}^H$ is the orthogonal projector corresponding to the set $\mathbf{\Phi}_k$ containing $k$ angles from $\Psi$. Computing $\lambda_M^-(\tilde{u})$ and $\lambda_M^+(\tilde{u})$ following those formula requires the resolution of the $l_0$-constrained problem for all degrees of sparsity: it is combinatorial and totally unproductive. For those reasons, we propose a statistical analysis of the bounds $\lambda_M^-$ and $\lambda_M^+$ in the case $M=2$ signals. Those statistics necessarily depend on the statistics of $\tilde{u}$ previously provided.

## 3.2. Approximated distribution of $\lambda_2^-$ and $\lambda_2^+$ for stochastic observations

We consider two signals incoming from directions $\theta_1$ and $\theta_2$. We observe that the values $c_k(\mathbf{y})$ of (15) corresponds to the minimum of the ML criterion considering $k$ incoming signals. In the noiseless case, it is thus obvious that $c_2(\mathbf{y}) = \tilde{u}^H \mathbf{\Pi}_n(\{\theta_1, \theta_2\}) \tilde{u}$. In the noisy case, it can be considered that:

$$ c_2(\mathbf{y}) \approx \tilde{u}^H \mathbf{\Pi}_n(\{\theta_1, \theta_2\}) \tilde{u}, \quad (17) $$

$$ c_3(\mathbf{y}) \approx \min_{\psi_g \neq \theta_1, \theta_2} \tilde{u}^H \mathbf{\Pi}_n(\{\theta_1, \theta_2, \psi_g\}) \tilde{u}. \quad (18) $$

$\lambda_2^-(\tilde{u})$ and $\lambda_2^+(\tilde{u})$, defined in (13) and (14), can then be approximated by the maximum and minimum, respectively, over hermitian quadratic forms as follows:

$$ \lambda_2^-(\tilde{u}) \approx \frac{c_2(\tilde{u}) - c_3(\tilde{u})}{2} \approx \frac{1}{2} \max_{g \mid \psi_g \neq \theta_1, \theta_2} Q_g^-(\tilde{u}) \quad (19) $$

$$ \lambda_2^+(\tilde{u}) \approx \frac{c_1(\tilde{u}) - c_2(\tilde{u})}{2} \approx \frac{1}{2} \min_g Q_g^+(\tilde{u}) \quad (20) $$

with

$$ Q_g^-(\tilde{u}) = \tilde{u}^H (\mathbf{\Pi}_n(\{\theta_1, \theta_2\}) - \mathbf{\Pi}_n(\{\theta_1, \theta_2, \psi_g\})) \tilde{u} \quad (21) $$

$$ Q_g^+(\tilde{u}) = \tilde{u}^H (\mathbf{\Pi}_n(\psi_g) - \mathbf{\Pi}_n(\{\theta_1, \theta_2\})) \tilde{u}. \quad (22) $$

Let us note $\mathbf{Q}^- = [Q_1^-, \dots, Q_g^-, \dots, Q_G^-]^T$, for $g|\psi_g \neq \theta_1, \theta_2$. From (19), we deduce that $\lambda_2^-(\tilde{u})$ is inferior to a value $t$ if, for all $g$ such that $\psi_g \neq \theta_1, \theta_2$, we
have $\frac{1}{2}Q_g^-(\tilde{u}) < t$. Mathematically, we thus express the distribution function $F_{\lambda_2^-}$ of $\lambda_2^-$ as:

$$F_{\lambda_2^-}(t) = \Pr(\lambda_2^-(\tilde{u}) \le t) = \Pr(Q_g^-(\tilde{u}) \le 2t, \forall g | \psi_g \ne \theta_1, \theta_2) = F_{Q^-}(2t, \dots, 2t), \quad (23)$$

where $F_{Q^-}$ is the multivariate distribution of the $(G-2)$-dimensional random vector $Q^-$.

In order to derive this multivariate distribution, we first express the marginal distributions of $Q_g^-$ for $g = 1 \dots G$, $g | \psi_g \ne \theta_1, \theta_2$. The expressions $Q_g^-(\tilde{u})$ can be transformed in real quadratic forms [19] whose variable statistics are computed from the previously provided statistics u, $\Gamma$ and $C$ given by equations (9), (11) and (12). Details are given in Appendix A.1. Distributions of $Q_g^-$ can then be approximated as non-central chi-square distributions. We then approximate the distributions of $Q_g^-$ as non-central chi-square distributions following the method in [20]: readers are referred to Appendix A.2 for details on the exact computation of the different parameters of the distributions. Let us note $F_g^-$ and $f_g^-$ the obtained cumulative and probability density functions of $Q_g^-$, $g = 1 \dots G$. Unfortunately, those variables are not independent. A copula function $c$ can be used in order to derive the multivariate distribution $F_{Q^-}$ of the $(G-2)$-dimensional random vector $Q^-$ with the desired marginal distributions [21]:

$$F_{Q^{-}}(\varrho_{1}, \dots, \varrho_{G}) = c(F_{1}^{-}(\varrho_{1}), \dots, F_{G}^{-}(\varrho_{G})). \quad (24)$$

The Gaussian copula is particularly convenient in high dimension, and the dependence between variables is modeled by a symmetric and definite positive matrix $\Sigma$. The Gaussian copula is a function defined as follows:

$$c(\zeta_1, \dots, \zeta_G) = \frac{1}{\sqrt{\det \Sigma}} \exp \left( -\frac{1}{2} \beta^T (\Sigma^{-1} - I) \beta \right) \quad (25)$$

with $\beta^T = [\Phi^{-1}(\zeta_1), \dots, \Phi^{-1}(\zeta_G)]$, $\Phi^{-1}$ being the inverse cumulative distribution function of a standard normal distribution. The $ij^{th}$ component of the correlation matrix $\Sigma$ is given by:

$$\Sigma_{ij} = \frac{\mathrm{Cov}(Q_i^-, Q_j^-)}{\sigma_{Q_i^-} \sigma_{Q_j^-}}. \quad (26)$$

Let us write, for $g = 1 \dots G$, $P_g^-$ = $\Pi_n(\{\theta_1, \theta_2\}) - \Pi_n(\{\theta_1, \theta_2, \psi_g\})$ and thus $Q_g^- = \tilde{u}^H P_g^- \tilde{u}$. The covariance between $Q_i^-$ and $Q_j^-$ is computed as:

$$\begin{aligned} \mathrm{Cov}(Q_i^-, Q_j^-) ={}& \mathrm{tr}((\mathbf{P}_i^- \mathbf{C})(\mathbf{P}_j^-)^T \mathbf{C}^*) \\ & + \mathrm{tr}((\mathbf{P}_i^- \mathbf{\Gamma})(\mathbf{P}_j^- \mathbf{\Gamma})) + \mathbf{u}^H \mathbf{P}_i^- \mathbf{\Gamma} \mathbf{P}_j^- \mathbf{u} + \mathbf{u}^H \mathbf{P}_j^- \mathbf{\Gamma} \mathbf{P}_i^- \mathbf{u} \\ & + \mathbf{u}^H \mathbf{P}_i^- \mathbf{C} (\mathbf{P}_j^-)^T \mathbf{u}^* + \mathbf{u}^T (\mathbf{P}_j^-)^T \mathbf{C}^* \mathbf{P}_i^- \mathbf{u}. \end{aligned} \quad (27)$$
The standard deviation $\sigma_{Q_g^-}$ is computed using $\sigma_{Q_g^-}^2 = \text{Cov}(Q_g^-, Q_g^-)$ given from equation (27). Steps (24)-(27) are injected into (23) to derive the statistics of $\lambda_2^-$. Similar steps are used to compute the statistics of $\lambda_2^+$ from (20), with:

$$F_{\lambda_2^+}(t) = \Pr(\lambda_2^+(\tilde{\mathbf{u}}) \le t) = \Pr(\exists g | Q_g^+(\tilde{\mathbf{u}}) \le 2t) \quad (28)$$

$$= \sum_{g=1}^{G} F_{Q_g^+}(2t) + \dots + (-1)^{k-1} \sum_{1 \le g_1 < \dots < g_k \le G} F_{(Q_{g_1}, \dots, Q_{g_k})}(2t, \dots, 2t) \\ + \dots + (-1)^{G-1} F_{Q^+}(2t, \dots, 2t). \quad (29)$$

# 4. Numerical study

Let us consider a Uniform Circular Array (UCA) with $N = 7$ antennas and radius $r = 1/2\lambda_0$, where $\lambda_0$ is the wavelength. Although UCAs allow for a two dimensional estimation, we limit ourselves in the sequel of this paper to the one dimensional estimation problem (azimuth angle): we assume that the sources and sensors are coplanar. The field-of-view is the range $[0, 360]^{\circ}$ with a grid spacing of $1^{\circ}$. UCAs have the interesting property of being omni-directional and present $\theta$ invariant performance. For two given angles of arrival $\theta_1$ and $\theta_2$, we are thus interested in the spacing $|\theta_2 - \theta_1|$ rather than with the angles themselves. The null-to-null spacing of the main lobe (FNBW) of this array is $90^{\circ}$: a scenario in which directions of arrival are separated from less than half of this main lobe is therefore called a close sources scenario, as conventional methods can have difficulties to detect the two signals.

## 4.1. Behavior of CEL0 minimization algorithm

In part because of the numerous local minima, the minimization of the $\ell_0$-regularized criterion is particularly challenging. Results obtained by the proximal-gradient algorithm IHT are highly dependent on the noise and initialization, thus results can be inconsistent. Continuous relaxation and corresponding algorithms have been proposed to improve the robustness and performance. In this paper, we will consider the CEL0-regularized criterion, which is an exact relaxation [22] of the $\ell_0$-regularized cost function:

$$J_{\text{CEL0}}(\lambda, \mathbf{z}) = \frac{1}{2} \|\mathbf{A}\mathbf{z} - \tilde{\mathbf{u}}\|_2^2 + \Phi_{\text{CEL0}}(\lambda, \mathbf{z}). \quad (30)$$

The CEL0 penalty is defined, for normalized column in $\mathbf{A}$, by:

$$\Phi_{\text{CEL0}}(\lambda, \mathbf{z}) = G\lambda - \sum_{i \in \mathbb{I}_G} \frac{(|z_i| - \sqrt{2\lambda})^2}{2} \mathbf{1}_{\{|z_i| \le \sqrt{2\lambda}\}} \quad (31)$$

The CEL0 criterion (30) presents less local minima than the $\ell_0$ one, and preserves the global minimum [22]. We consider the adaptation of iterative reweighted $\ell_1$ algorithm proposed in [23] to implement the CEL0 minimization:
240 this algorithm is called IRL1-CEL0. It is indeed show that it can yield consistent performance, outperforming some other sparse methods [5]. We provide the steps of IRL1-CEL0 as a reminder in Algorithm 1. A post-processing step is then applied to detect two peaks in the estimated spectrum.

**Algorithm 1 IRL1-CEL0**

**Input:** dictionary **A** of unitary columns, observation **u**, $\beta = 1/L$, $L$ Lipschitz-constant of **A**, stopping criteria $\epsilon$, $n$, $\epsilon_{outer}$, $n_{outer}$

**Initialization:** $\mathbf{z}^{(0)} = \mathbf{0}$, $J^{(0)} = \tilde{\mathbf{u}}^H \tilde{\mathbf{u}}$, $i = 0$, $i_{outer} = 0$

while $\frac{|J^{(i)} - J^{(i-1)}|}{|J^{(i)}|} > \epsilon_{outer}$ and $i_{outer} < n_{outer}$ do

$i_{outer} = i_{outer} + 1$

• Compute **w** weighting vector by:

$$w_g = (\sqrt{2\lambda} - |z_g^{(i)}|) \mathbf{1}_{\{|z_g^{(i)}| \le \sqrt{2\lambda}\}}, g = 1..G$$

• $j = 0$, $T^{(0)} = 1$, $\mathbf{z}_{temp}^{(0)} = \mathbf{z}^{(i)}$

while $\frac{\|\mathbf{z}^{(i)} - \mathbf{z}^{(i-1)}\|_2^2}{\|\mathbf{z}^{(i-1)}\|_2^2} > \epsilon$ and $j < n$ do

$j = j + 1$, $i = i + 1$

• $\mathbf{z}^{(i)} = \text{prox}_{\|\cdot\|_1, \lambda\beta\mathbf{w}}(\mathbf{z}_{temp}^{(j-1)} - \beta\mathbf{A}^H(\mathbf{A}\mathbf{z}_{temp}^{(j-1)} - \tilde{\mathbf{u}})),$

with $[\text{prox}_{\|\cdot\|_1, \lambda\beta\mathbf{w}}]_g(\mathbf{v}) = \max\{0, 1 - \frac{\lambda\beta w_g}{|v_g|}\} \cdot v_g \mathbf{1}_{\{|v_g| \ne 0\}}$

• $T^{(j)} = \frac{1+\sqrt{1+(2T^{(j-1)})^2}}{2}$

• $\mathbf{z}_{temp}^{(j)} = \mathbf{z}^{(i)} + \frac{T^{(j-1)}-1}{T^{(j)}} (\mathbf{z}^{(i)} - \mathbf{z}^{(i-1)})$

• $J^{(i)} = \|\mathbf{A}\mathbf{z}^{(i)} - \tilde{\mathbf{u}}\|_2^2$

end while

end while

return $\mathbf{z}^{(i)}$

Performance of this algorithm as a function of $\lambda$ is represented on Fig. 1. There are compared to the histograms of $\lambda_2^-$ and $\lambda_2^+$. Performance is presented first in term of % of errors then in terms of RMSE. An error means that we fail to detect two peaks in the estimated pseudo-spectrum $\mathbf{z}$. For well separated histograms of $\lambda_2^-$ and $\lambda_2^+$ ($|\theta_2 - \theta_1| = 48^\circ$, which is just over half the beamwidth), the interval between the histograms corresponds to an interval for which the percentage of errors is 0 (represented by dotted lines on Fig. 1 right). We deduce that, in this interval, the algorithm effectively reaches a solution close to the global solution. There is therefore, for separated enough signals, a theoretical interval $I(|\theta_2 - \theta_1|)$ in which the regularization parameter $\lambda$ should stand to estimate $\mathbf{z}_0$ from problem (2) for those two given angles of arrival. However, for close directions of arrival (as for example $|\theta_2 - \theta_1| = 28^\circ$), histograms start to overlap, resulting in the non existence of an interval presenting exactly 0% of errors.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Figure 1: Performance of CEL0 minimization algorithm as a function of $\lambda$ for two coherent signals, for $\sigma_e = 0.15$. Left: $|\theta_2 - \theta_1| = 28^\circ$. Right: $|\theta_2 - \theta_1| = 48^\circ$. Top: histograms (pdf) of $\lambda_2^-$ and $\lambda_2^+$ for stochastic observations. Middle and bottom: percentage of errors and RMSE.</figcaption></figure>

## 4.2. Off-line determination of $\lambda$

The interval $I(|\theta_2 - \theta_1|)$ will be defined between the 95% of the cumulative distribution function of $\lambda_2^-$ and 5% of the one of $\lambda_2^+$. As $\theta_1$ and $\theta_2$ are unknown, there is a need to study this interval for all candidate angles. The bounds of the interval $I$ depending on $|\theta_2 - \theta_1|$ are represented on Fig. 2. For close directions of arrival, the interval $I$ doesn't exist anymore. Let us note that the upper bound of $\lambda_2^-$ is proportional to $\sigma_e^2$.

<figure><img src="image_4.png" /><figcaption>Figure 2: Interval $I$ as a function of $|\theta_2 - \theta_1|$, for $\sigma_e = 0.15$.</figcaption></figure>

From Fig. 2, it is possible to select a unique value for $\lambda$ that will allow to estimate $z$ whatever the spacing between the signals, as long as they are not too close. This off-line, theoretical and a priori selection allows to perform DOA estimation for an unique and appropriate value of $\lambda$ and not for many different values, as opposed to other parameters selection techniques [24]. We set $\lambda$ at the intersection between the upper bound of $\lambda_2^-$ and the lower bound of $\lambda_2^+$: for $\sigma_e = 0.15$, we thus have $\lambda = 0.007$ for an intersection at $26^\circ$. This defines a theoretical resolution limit (i.e. the minimum angle separation for which the two signals are detected and separated) that is also verified by our simulations. Let
us remember that for such angle separation with modeling errors and highly-correlated signals, performance of MUSIC algorithm is already largely degraded, and that the method fails for coherent signals. ML methods, on the other hands, are particularly computationally heavy.

### 4.3. Performance comparison

In the previous paragraphs, we proposed a) the use of the eigenvector corresponding to the highest eigenvalue of the covariance matrix as an observation vector; b) a theoretical analysis of the regularization parameter $\lambda$ leading to its off-line selection. Let us now compare the performance of IRL1-CEL0 using our framework with others existing methods. Performance is compared with BPDN - using the solver SPGL1 [25, 26] and UCA-RB-MUSIC [27], a beamspace version of MUSIC with forward-backward averaging and spatial-smoothing developed for DOA estimation of correlated or coherent signals with UCAs. ML performance is also presented as a benchmark: in order to gain in simulation time, performance of the ML estimator is computed from the performance prediction in the presence of modeling errors given in [28]. We consider 500 independent trials. Results are presented as a function of $\sigma_e$ on Fig. 3 and as a function of the angles of arrival separation $|\theta_2 - \theta_1|$ on Fig. 4. Simulations are performed with MATLAB Version R2019a by using a laptop with a 2.70 GHz processor and 16 GB of memory. Average simulation time of IRL1-CEL0, BPDN and ML has been computed on 500 trials for $\sigma_e = 0.15$ and $|\theta_2 - \theta_1| = 30^{\circ}$: results are presented on Table 1. The maximum number of iterations is fixed to the same value for IRL1-CEL0 and BPDN.

Performance of IRL1-CEL0 algorithm is compared with a performance prediction of the Maximum Likelihood (ML) estimator in the presence of modeling errors given in. Results are presented as a function of the angles of arrival separation on Fig. 4 and as a function of $\sigma_e$ on Fig. 4. Performance of CEL0 algorithm coupled with the theoretical determination of $\lambda$ is really close to the provided bound, which is the benchmark for such scenario.

<figure><img src="image_5.png" /><figcaption>Figure 3: RMSE as a function of $\sigma_e$ for two coherent signals with $|\theta_2 - \theta_1| = 30^{\circ}$, $\theta_1 = 36^{\circ}$.</figcaption></figure>

Comparatively to UCA-RB-MUSIC and BPDN, IRL1-CEL0 based sparse DOA estimation provides performance closer to the predicted one of the more
<figure><img src="image_6.png" /><figcaption>Figure 4: RMSE as a function of the angles separation for two coherent signals impinging from directions $\theta_1$ and $\theta_2$. $\theta_1$ is fixed to 36°, and $\sigma_e = 0.15$.</figcaption></figure>

<table><caption>Table 1: Average estimation time (s) on 500 trials for $\sigma_e = 0.15$ and $|\theta_2 - \theta_1| = 30^{\circ}$.</caption><thead><tr><td>Method</td><td>Time (s)</td></tr></thead><tbody><tr><td>IRL1-CEL0</td><td>0.062</td></tr><tr><td>Basis Pursuit Denoising (BPDN)</td><td>0.014</td></tr><tr><td>Maximum Likelihood (ML)</td><td>0.925</td></tr></tbody></table>

computationally expensive ML method. The performance must indeed be analyzed in relation to the computing time. The computing time of ML is nearly 15 times that of IRL1-CEL0; IRL1-CEL0 computing time is only 4 times that of BPDN (using an optimized solver). As ML requires a combinatorial search, it is obvious that the more multipath they are, the more interesting our framework will be: its complexity does not depend on the number of signals. This advantage might be investigated in a future work.

## 5. Conclusion

For coherent signals in multipath environment, the first covariance eigen-vector is used as an observation for a sparse DOA estimation. Performance of sparse ℓ₀-regularized DOA estimation algorithms is highly dependent on the choice of the regularization parameter $\lambda$. For the proposed observation in the presence of modeling errors, we statistically compute an interval $I$ in which the regularization parameter $\lambda$ must stand so that the theoretical solution of the ℓ₀-regularized problem corresponds to the ML solution. Simulations confirm the relevance of the selected interval. Indeed, for $\lambda$ belonging to $I$, CELO-regularized sparse DOA estimation provides performance very close to the predicted one of the more computationally expensive ML methods. We propose a strategy for an off-line determination of $\lambda$ in the admissible interval that defines a new theoretical resolution limit and we show the relevance of this choice. IRL1-CEL0 algorithm along with our framework comes closer to the predicted performance of the more computationally expensive ML methods. Moreover, unlike ML, the complexity of IRL1-CEL0 using our framework does not depend on the number of paths: our work therefore pave the way
Appendix A.

Appendix A.1. Transformation of a hermitian quadratic form into a real quadratic form

330

Let us consider a hermitian quadratic form $Q(\tilde{u}) = \tilde{u}^H P \tilde{u}$, with $\tilde{u} \sim CN(u, \Gamma, C)$ and $P$ a hermitian matrix. Let us define the real vector $\tilde{u}_R$ as:

$$ \tilde{u}_R = [\Re(\tilde{u})^T, \Im(\tilde{u})^T]^T. \qquad (A.1) $$

Basic calculations show that $\tilde{u}_R$ follows a multivariate normal distribution $N(u_R, \Gamma_R)$:

335

$$ u_R = [\Re(u)^T, \Im(u)^T]^T, \qquad (A.2) $$

$$ \Gamma_R = \frac{1}{2} \begin{pmatrix} \Re(\Gamma + C) & \Im(-\Gamma + C) \\ \Im(\Gamma + C) & \Re(\Gamma - C) \end{pmatrix}. \qquad (A.3) $$

Then the quadratic form $Q(\tilde{u})$ can be expressed as [19]:

$$ Q(\tilde{u}) = \tilde{u}_R^T P_R \tilde{u}_R \qquad (A.4) $$

with

$$ P_R = \begin{pmatrix} (P + P^T)/2 & i(P - P^T)/2 \\ i(P^T - P)/2 & (P + P^T)/2 \end{pmatrix}. \qquad (A.5) $$

$P_R$ is a symmetric matrix.

Appendix A.2. Non-central chi-square distribution approximation for a real quadratic form

340

Let us consider the quadratic form $Q(\tilde{u}) = \tilde{u}^H P \tilde{u} = \tilde{u}_R^T P_R \tilde{u}_R$ as given in the previous paragraph. The mean, standard deviation, skewness and kurtosis of $Q(\tilde{u})$ are respectively given by:

$$ \mu_Q = s_1, \quad \sigma_Q = \sqrt{2s_2}, \quad \beta_1 = \sqrt{8}b_1, \quad \beta_2 = 12b_2 \qquad (A.6) $$

with $s_k = \text{Tr}((P_R \Gamma_R)^k) + k u_R^T (P_R \Gamma_R)^{k-1} P_R u_R$, $b_1 = s_3/s_2^{3/2}$ and $b_2 = s_4/s_2^2$.

345

The distribution of $Q(\tilde{u})$ can be approximated by a non-central chi-square distribution $\chi_l^2(\delta)$ as:

$$ \Pr(Q(\tilde{u}) > t) \approx \Pr\left(\frac{\chi_l^2(\delta) - \mu_\chi}{\sigma_\chi} > \frac{t - \mu_Q}{\sigma_Q}\right), \qquad (A.7) $$

where

$$ \mu_\chi = \mathbb{E}\{\chi_l^2(\delta)\} = l + \delta \qquad (A.8) $$

$$ \sigma_\chi = \sqrt{\text{var}(\chi_l^2(\delta))} = \sqrt{2a}, \quad \text{with} \quad a = \sqrt{l + 2\delta}. \qquad (A.9) $$

350

The non-centrality parameter $\delta$ and the number of degrees of freedom $l$ are computed so that it matches the third-order moment and has the best match with the fourth-order moment [20]:

* if $b_1^2 > b_2$, $\delta = b_1 a^3 - a^2$ and $l = a^2 - 2\delta$ with $a = 1/(b_1 - \sqrt{b_1^2 - b_2})$,

* if $b_1^2 \le b_2$, $\delta = 0$ and $l = 1/b_1^2$.
References

[1] P. Stoica, A. Nehorai, Music, maximum likelihood, and cramer-rao bound, IEEE Transactions on Acoustics, Speech, and Signal Processing 37 (5) (1989) 720-741. doi:10.1109/29.17564.

[2] S. Pillai, B. Kwon, Forward/backward spatial smoothing techniques for coherent signal identification, IEEE Transactions on Acoustics, Speech, and Signal Processing 37 (1) (1989) 8-15. doi:10.1109/29.17496.

[3] Z. Yang, J. Li, P. Stoica, L. Xie, Chapter 11 - sparse methods for direction-of-arrival estimation, in: R. Chellappa, S. Theodoridis (Eds.), Academic Press Library in Signal Processing, Volume 7, Academic Press, 2018, pp. 509-581. doi:https://doi.org/10.1016/B978-0-12-811887-0.00011-0.

[4] X. Wu, W.-P. Zhu, J. Yan, A high-resolution doa estimation method with a family of nonconvex penalties, IEEE Transactions on Vehicular Technology 67 (6) (2018) 4925-4938. doi:10.1109/TVT.2018.2817638.

[5] E. Soubies, A. Chinatto, P. Larzabal, J. M. T. Romano, L. Blanc-Féraud, Direction-of-arrival estimation through exact continuous $ℓ_{2,0}$-norm relaxation, IEEE Signal Processing Letters 28 (2021) 16-20. doi:10.1109/LSP.2020.3042771.

[6] D. L. Donoho, M. Elad, Optimally sparse representation in general (nonorthogonal) dictionaries via 11 minimization, Proceedings of the National Academy of Sciences of the United States of America 100 (5) (2003) 2197-2202.

[7] J. A. Tropp, Greed is good: Algorithmic results for sparse approximation, IEEE Transactions on Information Theory 50 (10) (2004) 2231-2242.

[8] B. Rao, K. Engan, S. Cotter, J. Palmer, K. Kreutz-Delgado, Subset selection in noise based on diversity measure minimization, IEEE Transactions on Signal Processing 51 (3) (2003) 760-770. doi:10.1109/TSP.2002.808076.

[9] C. Soussen, J. Idier, J. Duan, D. Brie, Homotopy based algorithms for $ℓ_0$-regularized least-squares, IEEE Transactions on Signal Processing 63 (13) (2015) 3301-3316. doi:10.1109/TSP.2015.2421476.

[10] B. Friedlander, A sensitivity analysis of the music algorithm, in: International Conference on Acoustics, Speech, and Signal Processing, 1989, pp. 2811-2814 vol.4. doi:10.1109/ICASSP.1989.267053.

[11] B. Friedlander, Sensitivity analysis of the maximum likelihood direction-finding algorithm, IEEE Transactions on Aerospace and Electronic Systems 26 (6) (1990) 953-968. doi:10.1109/7.62247.

[12] A. Ferreol, P. Larzabal, M. Viberg, On the resolution probability of music in presence of modeling errors, IEEE Transactions on Signal Processing 56 (5) (2008) 1945-1953.
395

[13] M. Nikolova, Relationship between the optimal solutions of least squares regularized with L0-norm and constrained by k-sparsity, Applied and Computational Harmonic Analysis 41 (1) (2016) 237-265. doi:10.1016/j.acha.2015.10.010.

400

[14] A. Delmer, A. Ferréol, P. Larzabal, On regularization parameter for L0-sparse covariance fitting based doa estimation, in: 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2020, pp. 4552-4556. doi:10.1109/ICASSP40776.2020.9053963.

[15] A. Delmer, A. Ferréol, P. Larzabal, L0-sparse doa estimation of close sources with modeling errors, in: 2020 28th European Signal Processing Conference (EUSIPCO), 2021, pp. 1861-1865. doi:10.23919/Eusipco47968.2020.9287732.

405

[16] A. Ferreol, P. Larzabal, M. Viberg, On the asymptotic performance analysis of subspace doa estimation in the presence of modeling errors: case of music, IEEE Transactions on Signal Processing 54 (3) (2006) 907-920. doi:10.1109/TSP.2005.861798.

410

[17] A. Ferreol, P. Larzabal, M. Viberg, Modeling error sensitivity of the music algorithm conditioned on resolved sources, in: 16th European Signal Processing Conference (EUSIPCO 2008), EURASIP, 2008.

[18] J.-J. Fuchs, Doa estimation in the presence of modeling errors, the global matched filter approach, in: 17th European Signal Processing Conference (EUSIPCO 2009), EURASIP, 2009, pp. 1963-1967.

415

[19] A. A. Mohsenipour, S. B. Provost, On the Distribution of Quadratic Expressions in Various Types of Random Vectors, Ph.D. thesis (2012). URL https://ir.lib.uwo.ca/etd://ir.lib.uwo.ca/etd/955

420

[20] H. Liu, Y. Tang, H. H. Zhang, A new chi-square approximation to the distribution of non-negative definite quadratic forms in non-central normal variables, Computational Statistics and Data Analysis 53 (4) (2009) 853-856. doi:10.1016/j.csda.2008.11.025.

[21] R. Nelsen, An Introduction to Copulas, 2nd Edition, Springer Series in Statistics, Springer, 2006.

425

[22] E. Soubies, L. Blanc-Féraud, G. Aubert, A continuous exact L0 penalty (cel0) for least squares regularized problem, SIAM Journal on Imaging Sciences 8 (3) (2015) 1607-1639. doi:10.1137/151003714.

[23] S. Gazagnes, E. Soubies, L. Blanc-Féraud, High density molecule localization for super-resolution microscopy using cel0 based sparse approximation, in: 14th International Symposium on Biomedical Imaging (ISBI 2017), IEEE, 2017, pp. 28-31.

430
435

[24] A. Panahi, M. Viberg, Performance analysis of sparsity-based parameter estimation, IEEE Transactions on Signal Processing 65 (24) (2017) 6478–6488. doi:10.1109/TSP.2017.2755602.

[25] E. van den Berg, M. P. Friedlander, SPGL1: A solver for large-scale sparse reconstruction, https://friedlander.io/spgl1 (December 2019).

[26] E. van den Berg, M. P. Friedlander, Probing the pareto frontier for basis pursuit solutions, SIAM Journal on Scientific Computing 31 (2) (2008) 890–912. doi:10.1137/080714488. URL http://link.aip.org/link/?SCE/31/890

440

[27] C. Mathews, M. Zoltowski, Eigenstructure techniques for 2-d angle estimation with uniform circular arrays, IEEE Transactions on Signal Processing 42 (9) (1994) 2395–2407. doi:10.1109/78.317861.

[28] A. Ferreol, P. Larzabal, M. Viberg, Performance prediction of maximum-likelihood direction-of-arrival estimation in the presence of modeling errors, IEEE Transactions on Signal Processing 56 (10) (2008) 4785–4793.

445